SUBROUTINE LINSET(ITR,IUNIT,IFRQ0,IFRQ1,XMAX,DOP,AGAM) C ====================================================== C C Set up frequency points and weights for a line C Auxiliary procedure for START C C Input: C ITR - index of the transition C INMOD - mode of evaluating frequency points in the line C = 0 - means that frequency points and weights have C already been read; C ne 0 - frequency points and weights will be evaluated: C the meaning of the individual values: C = 1 - equidistant frequencies, trapezoidal integration C = 2 - equidistant frequencies, Simpson integration C = 3 - modified Simpson integration, i.e. a series of C 3-point Simpson integrations with each subsequent C integration interval doubled, until the whole C integartion area is covered C = 4 - frequencies (in units of standard x) and weights C (for integration over x) are read; C C XMAX > 0 - means that the line is assumed symmetric around the C center, frequency points are set up between x=0 and C x=XMAX, where x is frequency difference from the C line center in units of the standard Doppler width C < 0 - frequency points are set between x=XMAX and x=-XMAX C DOP - Doppler width C AGAM - damping parameter (for lines with Voigt or non-standard C profile only) C C Output (to COMMON/FRQEXP) C FREQ - array of frequencies C Note that LINSET calculates values of frequencies only C for frequency points between IFR0(ITR) and IFR1(ITR). C W - corresponding integration weights C PROF - corresponding values of absorption profile C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' PARAMETER (BOL2=2.76108D-16, CIN=UN/2.997925D10, OS0=0.02654, * F0C1=1.25D-9, TTW=2./3., PISQ1=UN/1.77245385090551D0, * C18IN=UN/CAS,F0C2=3.906D-11) DIMENSION X(MFREQL),W0(MFREQL) C IF(ITR.EQ.0) GO TO 200 INMOD=INTMOD(ITR) INMOD0=MOD(IABS(INMOD),10) IJ0=IFR0(ITR) IJ1=IFR1(ITR) N=IJ1-IJ0+1 IF(INMOD.EQ.0) THEN S=OSC0(ITR)*0.02654D0 IP0=IPROF(ITR) IP=IABS(IP0) DO I=1,N PROF(I+IJ0-1)=PROFIL(FREQ(I+IJ0-1),AGAM,DOP,ITR,IP,0)* * S/DOP IJLIN(I+IJ0-1)=ITR END DO IF(IP0.GE.0) THEN IF(XMAX.LT.0.) THEN PROF(IJ0)=0. PROF(IJ1)=0. ELSE PROF(IJ1)=0. END IF END IF IF(IPROF(ITR).GE.0) THEN IF(XMAX.LT.0.) THEN INTMOD(ITR)=-2 ELSE INTMOD(ITR)=-1 END IF ELSE IF(XMAX.LT.0.) THEN INTMOD(ITR)=2 ELSE INTMOD(ITR)=1 END IF END IF RETURN END IF X0=0. IF(XMAX.LT.0) X0=XMAX M=(N-1)/2 X(1)=0. W0(1)=UN IF(N.LE.1) GO TO 100 C IF(INMOD0.LE.2) THEN C C Trapezoidal integration C HH=ABS(X0+XMAX)/(N-1) DO I=1,N X(I)=X0+(I-1)*HH END DO IF(INMOD0.EQ.2) GO TO 40 DO I=1,N W0(I)=HH END DO W0(1)=HALF *HH W0(N)=HALF *HH GO TO 100 C C Ordinary Simpson integration C 40 HH=HH/3.D0 IF(MOD(N,2).NE.1) * CALL QUIT('even number of points in Simpson - LINSET',n,n) DO I=1,M W0(2*I)=4.D0*HH W0(2*I+1)=2.D0*HH END DO W0(1)=HH W0(N)=HH ELSE IF(INMOD0.EQ.3) THEN C C Modified Simpson integration - a set of 3-point Simpson C integrations with continuosly increasing distance between the C integration points (schematically, distances between points are C h,h,2h,2h,4h,4h, etc.) C TWI=UN MM=M IF(MOD(N,2).NE.1) * CALL QUIT('even number of points in MSimpson - LINSET',n,n) IF(XMAX.LT.0) MM=M/2 DO I=1,MM TWI=TWI*2.D0 X(2*I+1)= TWI-UN X(2*I)=TWI-UN-TWI/4.D0 W0(2*I)=2.D0*TWI W0(2*I+1)=1.5D0*TWI END DO TWN=TWI-UN TWA=ABS(XMAX)/TWN HH=TWA/6.D0 DO I=1,MM X(2*I+1)=X(2*I+1)*TWA X(2*I)=X(2*I)*TWA W0(2*I)=W0(2*I)*HH W0(2*I+1)=W0(2*I+1)*HH END DO W0(1)=HH W0(N)=TWI*HH/2.D0 X(1)=0. IF(M.EQ.MM) GO TO 100 IF(MOD(N,4).NE.1) * CALL QUIT('conflict in MSimpson - LINSET',n,n) DO I=1,M X(M+1+I)=X(I+1) W0(M+1+I)=W0(I+1) END DO M2=2*(M+1) DO I=1,M X(I)=-X(M2-I) W0(I)=W0(M2-I) END DO X(M+1)=0. W0(M+1)=2.D0*HH C C frequencies (in units of standard x) and weights C (for integration over x) are read; C ELSE IF(INMOD0.EQ.4) THEN READ(IUNIT,*) (X(I),I=1,N),(W0(I),I=1,N) END IF C C For all types of integration: C set up arrays FR,WW,PRF C 100 S=OSC0(ITR)*0.02654D0 IP0=IPROF(ITR) IP=IABS(IP0) DO I=1,N FREQ(I+IJ0-1)=FR0(ITR)-DOP*X(I) W(I+IJ0-1)=DOP*W0(I) PROF(I+IJ0-1)=PROFIL(FREQ(I+IJ0-1),AGAM,DOP,ITR,IP,0)*S/DOP END DO C C for IPROF(ITR) ge 0 - endpoint(s) of the line profile are forced to C have zero cross-section C IF(IP0.GE.0) THEN IF(XMAX.LT.0.) THEN PROF(IJ0)=0. PROF(IJ1)=0. ELSE PROF(IJ1)=0. END IF END IF C C Recalculation of quadrature weights in order to enforce exact C normalization of the integral (absorption profile * weights) C SUM=0. DO I=1,N SUM=SUM+PROF(I+IJ0-1)*W(I+IJ0-1) END DO SUM=S/SUM DO I=1,N W(I+IJ0-1)=W(I+IJ0-1)*SUM END DO C C reset the switch INTMOD C IF(IPROF(ITR).GE.0) THEN IF(XMAX.LT.0.) THEN INTMOD(ITR)=-2 ELSE INTMOD(ITR)=-1 END IF ELSE IF(XMAX.LT.0.) THEN INTMOD(ITR)=2 ELSE INTMOD(ITR)=1 END IF END IF IF(ABS(INMOD).GE.10) INTMOD(ITR)=0 C IF(INDEXP(ITR).NE.0) THEN CALL IJALIS(ITR,IFRQ0,IFRQ1) END IF RETURN C 200 CONTINUE IF(ISPODF.GT.0) RETURN DO 220 IT=1,NTRANS IF(.NOT.LINE(IT)) GO TO 220 IP=IABS(IPROF(IT)) IF(IP.NE.2) GO TO 220 IAT=IATM(ILOW(IT)) AM=BOL2/AMASS(IAT)*TEFF DOPP=FR0(IT)*CIN*SQRT(AM+VTB*VTB) DOP1=UN/DOPP ANE=ELSTD IF(ANE.LE.0.) ANE=1.E14 F000=EXP(TTW*LOG(ANE)) II=NQUANT(ILOW(IT)) JJ=NQUANT(IUP(IT)) IZZ=IZ(IEL(ILOW(IT))) FAC=TWO F00=F0C1*F000 IF(IZZ.EQ.2) THEN FAC=UN F00=F0C2*F000 END IF CALL STARK0(II,JJ,IZZ,XKIJ,WL0,FIJ) FXK=F00*XKIJ DBETA=WL0*WL0*C18IN/FXK BETAD=DOPP*DBETA FID=OS0*FIJ*DBETA FID0=OS0*(OSC0(IT)-FIJ)*DOP1*PISQ1 CALL DIVSTR(IZZ) DO IJ=IFR0(IT),IFR1(IT) BETA=DBETA*ABS(FREQ(IJ)-FR0(IT)) SG=STARKA(BETA,fac)*FID SG0=0. V=(FREQ(IJ)-FR0(IT))*DOP1 IF(ABS(V).LE.13.) SG0=EXP(-V*V)*FID0 PROF(IJ)=SG+SG0 END DO 220 CONTINUE RETURN END