aboutsummaryrefslogtreecommitdiff
path: root/synthe/rh2ofast.for
blob: ae3f5ee8aefa5315d7dc3df754b3513385f4fd0b (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
PROGRAM RH2OFAST
C     FAST VERSION,  NO ENERGY LEVEL INFORMATION
C     READS PACKED BINARY VERSION OF PARTRIDGE AND SCHWENKE'S H2O LINELIST
      PARAMETER (kw=99)
      COMMON /LINDAT/WL,E,EP,LABEL(2),LABELP(2),OTHER1(2),OTHER2(2),
     1        WLVAC,CENTER,CONCEN, NELION,GAMMAR,GAMMAS,GAMMAW,REF,
     2      NBLO,NBUP,ISO1,X1,ISO2,X2,GFLOG,XJ,XJP,CODE,ELO,GF,GS,GR,GW,
     3        DWL,DGFLOG,DGAMMAR,DGAMMAS,DGAMMAW,DWLISO,ISOSHIFT,EXTRA3
      REAL*8 LINDAT8(14)
      REAL*4 LINDAT4(28)
      EQUIVALENCE (LINDAT8(1),WL),(LINDAT4(1),NELION)
      REAL*8 WL,E,EP,WLVAC,CENTER,CONCEN,UNPACKWL,WLVAC1
      REAL*8 LABEL,LABELP,OTHER1,OTHER2,LABELISO(4)
      CHARACTER*10 COTHER1,COTHER2,CLABEL,CLABELP
      EQUIVALENCE (COTHER1,OTHER1(1)),(COTHER2,OTHER2(1))
      EQUIVALENCE (CLABEL,LABEL(1)),(CLABELP,LABELP(1))
      INTEGER TYPE
      EQUIVALENCE (GF,G,CGF),(TYPE,NLAST)
      REAL*8 RESOLU,RATIO,RATIOLG,WLBEG,WLEND,RATIOLOG,WLBEG1,WLEND1
      REAL*8 RATIORATIO,VACAIR
      REAL*4 DECKJ(7,kw),XISO(4),X2ISO(4)
      REAL*4 TABLOG(32768),AIRSHIFT(100000)
      INTEGER*2 IELO,IGFLOG
      EQUIVALENCE (IWL,IWLBYTES(1))
      BYTE IWLBYTES(4),IELOBYTES(2),IGFLOGBYTES(2),ONEBYTE
      EQUIVALENCE (IELOBYTES(1),IELO),(IGFLOGBYTES(1),IGFLOG)
C               1H1H16O 1H1H17O 1H1H18O 1H2H16O
      DATA XISO/  .9976,  .0004,  .0020, .00001/
      DATA X2ISO/-0.001, -3.398, -2.690, -5.000/
      DATA LABELISO/2H16,2H17,2H18,2H26/
      data alpha/0./
C
      DO 1 I=1,32768
    1 TABLOG(I)=10.**((I-16384)*.001)
      IF(IFPRED.NE.1)CALL TABVACAIR(AIRSHIFT)
      OPEN(UNIT=11,STATUS='OLD',READONLY,SHARED,FORM='UNFORMATTED',
     1RECORDTYPE='FIXED',RECORDSIZE=2,ACCESS='DIRECT')
      OPEN(UNIT=12,TYPE='OLD',FORM='UNFORMATTED',ACCESS='APPEND')
      OPEN(UNIT=14,TYPE='OLD',FORM='UNFORMATTED',ACCESS='APPEND')
      READ(93)NLINES,LENGTH,IFVAC,IFNLTE,N19,TURBV,DECKJ,IFPRED,
     1WLBEG,WLEND,RESOLU,RATIO,RATIOLG,CUTOFF,LINOUT
      RATIOLOG=LOG(1.D0+1.D0/2000000.D0)
      RATIORATIO=RATIOLOG/RATIOLG
      IXWLBEG=DLOG(WLBEG)/RATIOLG
      IF(DEXP(IXWLBEG*RATIOLG).LT.WLBEG)IXWLBEG=IXWLBEG+1
C
      WLBEG1=WLBEG-1.
      WLEND1=WLEND+1.
      N=0
      READ(11,REC=1)IWL
c     on some computers need byte rotation
c     onebyte=iwlbytes(1)
c     iwlbytes(1)=iwlbytes(4)
c     iwlbytes(4)=onebyte
c     onebyte=iwlbytes(2)
c     iwlbytes(2)=iwlbytes(3)
c     iwlbytes(3)=onebyte
      WLVAC=EXP(IWL*RATIOLOG)
      IF(IFVAC.NE.1)THEN
      KWL=WLVAC*10.+.5
      WLVAC=WLVAC+AIRSHIFT(KWL)
      ENDIF
      PRINT 3334,WLVAC
 3334 FORMAT(' FIRST LINE IS        1','   WL',F12.4)
      IF(WLVAC.GT.WLEND1)GO TO 21
      LENGTHFILE=65912356
      READ(11,REC=LENGTHFILE)IWL
c     on some computers need byte rotation
c     onebyte=iwlbytes(1)
c     iwlbytes(1)=iwlbytes(4)
c     iwlbytes(4)=onebyte
c     onebyte=iwlbytes(2)
c     iwlbytes(2)=iwlbytes(3)
c     iwlbytes(3)=onebyte
      WLVAC=EXP(IWL*RATIOLOG)
      IF(IFVAC.NE.1)WLVAC=VACAIR(WLVAC)
      PRINT 3335,LENGTHFILE,WLVAC
 3335 FORMAT(' LAST LINE IS ',I9,'   WL',F12.4)
      IF(WLBEG1.GT.WLVAC)GO TO 21
C     FIND THE FIRST LINE AFTER ISTART
      LIMITBLUE=1
      LIMITRED=LENGTHFILE
   12 NEWLIMIT=(LIMITRED+LIMITBLUE)/2
      PRINT 3333,LIMITBLUE,NEWLIMIT,LIMITRED
 3333 FORMAT(3I10)
      READ(11,REC=NEWLIMIT)IWL
c     on some computers need byte rotation
c     onebyte=iwlbytes(1)
c     iwlbytes(1)=iwlbytes(4)
c     iwlbytes(4)=onebyte
c     onebyte=iwlbytes(2)
c     iwlbytes(2)=iwlbytes(3)
c     iwlbytes(3)=onebyte
      WLVAC=EXP(IWL*RATIOLOG)
      IF(IFVAC.NE.1)WLVAC=VACAIR(WLVAC)
      IF(WLVAC.LT.WLBEG1)GO TO 13
      LIMITRED=NEWLIMIT
      IF(LIMITRED-LIMITBLUE.LE.1)GO TO 14
      GO TO 12
   13 LIMITBLUE=NEWLIMIT
      IF(LIMITRED-LIMITBLUE.LE.1)GO TO 14
      GO TO 12
   14 ISTART=NEWLIMIT
      PRINT 3333,LIMITBLUE,LIMITRED,NEWLIMIT
      WRITE(6,6)ISTART
    6 FORMAT(I10,14H IS FIRST LINE)
C
      DO 20 ILINE=ISTART,LENGTHFILE
      READ(11,REC=ILINE)IWL,IELO,IGFLOG
c     on some computers need byte rotation
c     onebyte=iwlbytes(1)
c     iwlbytes(1)=iwlbytes(4)
c     iwlbytes(4)=onebyte
c     onebyte=iwlbytes(2)
c     iwlbytes(2)=iwlbytes(3)
c     iwlbytes(3)=onebyte
c     onebyte=ielobytes(1)
c     ielobytes(1)=ielobytes(2)
c     ielobytes(2)=onebyte
c     onebyte=igflogbytes(1)
c     igflogbytes(1)=igflogbytes(2)
c     igflogbytes(2)=onebyte
c
      WLVAC=EXP(IWL*RATIOLOG)
      FREQ=2.99792458E17/WLVAC
      IF(IFVAC.NE.1)THEN
      KWL=WLVAC*10.+.5
      WLVAC=WLVAC+AIRSHIFT(KWL)
      ENDIF
      IF(WLVAC.GT.WLEND1)GO TO 21
      IXWL=DLOG(WLVAC)/RATIOLG+.5D0
      ISO=1
      IF(IELO.GT.0.AND.IGFLOG.GT.0)GO TO 19
      ISO=2
      IF(IELO.GT.0)GO TO 19
      ISO=3
      IF(IGFLOG.GT.0)GO TO 19
      ISO=4
   19 NELION=534
      ELO=ABS(IELO)
      IGFLOG=ABS(IGFLOG)
      NBUFF=IXWL-IXWLBEG+1
      CONGF=.01502*TABLOG(IGFLOG)/FREQ*XISO(ISO)
      FRQ4PI=FREQ*12.5664
C     GAMMAS=0
C     LOG GAMMAW=-7
C     IGR=
      IGS=1
      IGW=9384
      GAMMAR=2.223E13/WLVAC**2*.001
      GAMRF=GAMMAR/FRQ4PI
C     GAMRF=TABLOG(IGR)/FRQ4PI
      GAMSF=TABLOG(IGS)/FRQ4PI
      GAMWF=TABLOG(IGW)/FRQ4PI
      WRITE(12)NBUFF,CONGF,NELION,ELO,GAMRF,GAMSF,GAMWF,alpha
      NLINES=NLINES+1
      IF(N.EQ.0)WRITE(6,1919)WLVAC
 1919 FORMAT(F12.4)
      N=N+1
   20 CONTINUE
   21 LINOUT=-ABS(LINOUT)
      WRITE(6,22)N
   22 FORMAT(I10,13H IS LAST LINE)
      WRITE(6,1919)WLVAC
   25 WRITE(6,26)NLINES
   26 FORMAT(I10,25H LINES WRITTEN ON TAPE 12)
      REWIND 93
      WRITE(93)NLINES,LENGTH,IFVAC,IFNLTE,N19,TURBV,DECKJ,IFPRED,
     1WLBEG,WLEND,RESOLU,RATIO,RATIOLG,CUTOFF,LINOUT
      CALL EXIT
      END
      SUBROUTINE TABVACAIR(AIRSHIFT)
      REAL*4 AIRSHIFT(100000)
      REAL*8 WLVAC,VACAIR
      DO 1 IWL=1,1999
    1 AIRSHIFT(IWL)=0.
      DO 2 IWL=2000,100000
      WLVAC=IWL*.1
    2 AIRSHIFT(IWL)=VACAIR(WLVAC)-WLVAC
      RETURN
      END
      FUNCTION VACAIR(W)
      IMPLICIT REAL*8 (A-H,O-Z)
C     W IS VACUUM WAVELENGTH IN NM
      WAVEN=1.D7/W
      VACAIR=W/(1.0000834213D0+
     1 2406030.D0/(1.30D10-WAVEN**2)+15997.D0/(3.89D9-WAVEN**2))
C    1(1.000064328+2949810./(1.46E10-WAVEN**2)+25540./(4.1E9-WAVEN**2))
      RETURN
      END