关于断点的问题，希望有大佬能看下这个代码，稍微有点长

wangzb · 发表于 2020-6-25 18:42:59

[Fortran] 纯文本查看 复制代码

001!-----------------------------------------------------------------------------------------------------
002!------------------------------J.Cai added filenames start--------------------------------------------
003!-----------------------------------------------------------------------------------------------------
004   character(256),parameter :: kvsgFile='11kvsg-4.9-5.0-1-H2O-0.0p-CO2-1p-300K-11.dat' !
005   character(256),parameter :: kvsgqFile='11kvsgq-4.9-5.0-1-H2O-0.0p-CO2-1p-300K.dat' !
006   character(256),parameter :: abscFile='11absc-4.9-5.0-1-H2O-0.0p-CO2-1p-300K.dat' !
007   character(256),parameter :: CO2File='02-3.7-3.8.par'
008   character(256),parameter :: CH4File='F:\JavaHAWKS\HITRAN04\02_hit2010.par'
009   character(256),parameter :: H2OFile='01_0-30000_HITEMP-35.par'
010   !character(256),parameter :: H2OFile='/opt/HITRAN/hitran96/by_molec/01_hit96.par' 
011! J.Cai added filenames end
012! J.Cai comment out starts
013!-----------------------------------------------------------------------------------------------------
014!! Selection of g-values for numerical quadrature, using a Numerical Recipes routine
015!! If Numerical Recipes is not available, set nq=12, comment out the following 8 lines of code,
016!! and uncomment the 5-line REAL declaration following it
017!   REAL              :: gqs(nq),wqs(nq),kq(nq),aq(numt,nq),gq(nq),wq(nq),gaujac,alf=3.,bet=-.9,wsum
018!! Get quadrature coefficients from Numerical recipies
019!    wsum=0.
020!    CALL GAUJAC(gqs,wqs,nq,alf,bet)
021!          do iq=1,nq
022!            gq(iq)=0.5*(1.-gqs(iq))
023!            wq(iq)=wqs(iq)/(2.**(alf+bet+1)*gq(iq)**alf*(1.-gq(iq))**bet)
024!            wsum=wsum+wq(iq)
025!          enddo
026!! Correction to make sum(wq)=1
027!         wq=wq/wsum
028!! End quadrature coefficients from Numerical recipies
029! J.Cai comment out ends
030! Selection of precalculated g-values for numerical quadrature, for nq=12,alf=3.,bet=0.
031! J.Cai uncomment out starts
032    REAL              :: kq(nq),aq(numt,nq), &
033           gq(nq)=(/ 5.120075E-02,1.170678E-01,2.015873E-01,3.007074E-01,4.095012E-01,5.225285E-01,  &
034                     6.341280E-01,7.387071E-01,8.310236E-01,9.064499E-01,9.612060E-01,9.925594E-01/),&
035           wq(nq)=(/ 5.556622E-02,7.576839E-02,9.258290E-02,1.048306E-01,1.118451E-01,1.132605E-01,  &
036                     1.090012E-01,9.927844E-02,8.457905E-02,6.563999E-02,4.341329E-02,1.904792E-02/)
037! J.Cai uncomment out starts
038! PRESSURE IN atm
039   patm=P/1.01325
040! IF P-BASED absc is requested make sure it is not a mixture
041    IF(ipl==1) THEN
042        cnchk=xmfr(1)*xmfr(2)+xmfr(1)*xmfr(3)+xmfr(3)*xmfr(2)
043        IF(cnchk > 1.d-6) PAUSE 'ERROR: you cannot use pressure-based absorption coefficient for mixture!'
044    ENDIF
045! Open output files
046! Output file for k,g(T,k),a(T,k) for all numT temperatures and all n_pwrk k-values
047   OPEN(7,FILE=kvsgFile) ! J.Cai centralize filename
048! Header formatted for TECPLOT, for a numT of 4
049   write(7,6)
0506 format('VARIABLES = k,g0,g1,g2,g3,g4,a1,a2,a3,a4')
051! Output file for wq,gq(Tref),kq(gq),aq(T,kq) for all numT temperatures and all nq g-values
052   IF(ipr==2) THEN
053       OPEN(8,FILE=kvsgqFile,STATUS='unknown') ! J.Cai centralize filename
054! Header formatted for readability, for a numT of 4
055        write(8,8)
056   ENDIF
0578 format('    wq',9x,'gq',9x,'kq',8x,'aq1',8x,'aq2',8x,'aq3',8x,'aq4')
058! File containing absorption coefficient
059   !IF(iwr>0) OPEN(9,FILE='c:\absco\absch2o-0p-2000K.dat',STATUS='unknown')
060   IF(iwr>0) OPEN(9,FILE=abscFile,STATUS='unknown') ! J.Cai centralize filename
061   hck=hhh*ccc/kkk
062   hckt0=hck/T0
063   IF(iwr==2) THEN
064! READ (GAS-ONLY) ABSORPTION COEFFICIENT FROM FILE
065      read(9,97) dummy
066      read(9,98) number
067      read(9,92) dummy,wvnm_b2,wvnm_e2,wvnmst2
068      number2=(wvnm_e2-wvnm_b2)/wvnmst2+1
069      IF(number2 /= number) PAUSE 'bad data file'
070      read(9,*) (absc(i),i=1,number)
071      CLOSE(9)
072      DO i=1,number
073         wvnm(i)=wvnm_b2+(i-1)*wvnmst2
074      ENDDO
075   ELSE
076! CALCULATE NECESSARY WAVENUMBERS AND ABSORPTION COEFFICIENTS
077! WVNM(I): WAVENUMBER
078! ABSC(I): ABSORPTION COEFFICIENT; INITIALIZE FOR SOOT
079   number=(wvnm_e-wvnm_b)/wvnmst+1
080   IF(number>nabs) PAUSE 'increase nabs'
081   DO i=1,number
082    wvnm(i)=wvnm_b+(i-1)*wvnmst
083    absc(i)=0.d0
084   ENDDO
085  kplancksum=0.d0
086! Scan over gases
087  DO ifg=1,3
088  IF(xmfr(ifg)<.9d-3) CYCLE
089! Open HITEMP database
090  IF(ifg==1) THEN
091     lu=11
092     write(*,*) 'Reading CO2 data'
093     OPEN(lu,FILE=CO2File, action='READ')  ! J.Cai change path
094  ELSEIF(ifg==2) THEN
095     lu=12
096     write(*,*) 'Reading H2O data'
097     OPEN(lu,FILE=H2OFile, action='READ')  ! J.Cai change path
098  ELSE
099     lu=13
100     write(*,*) 'Reading CH4 data'
101     OPEN(lu,FILE=CH4File, action='READ') ! J.Cai change path
102  ENDIF
103! GET LINE INFORMATION FROM HITEMP (all lines wvnm_b.le.lmbda.le.wvnm_e)
104   CALL Getdata(lines)
105   CLOSE(lu)
106   write(*,*) 'Gas ',ifg,',  lines read: ', lines
107   mass=1.d3*molemass(ifg)/nnn
108   hckt=hck/Tref
109   hcktt0=hckt0-hckt
110! MULTIPLIER FOR PRESSURE-BASED ABSORPTION COEFFICIENT
111   cr=nnn/8.314d1/Tref
112! MULTIPLIER FOR LINEAR ABSORPTION COEFFICIENT
113   IF(ipl==0) cr=cr*xmfr(ifg)*P
114   IF(ifg==1) THEN
115    write(*,*) 'start of absc loop for co2'
116! VIBRATIONAL PARTITION FUNCTION FOR CO2
117    V0=(1d0-DEXP(-666d0*hckt0))**2*(1d0-DEXP(-2396d0*hckt0))*   &
118            (1d0-DEXP(-1351d0*hckt0))
119    V=(1d0-DEXP(-666d0*hckt))**2*(1d0-DEXP(-2396d0*hckt))*   &
120            (1d0-DEXP(-1351d0*hckt))
121    VV0=V/V0
122! ROTATIONAL PARTITION FUNCTION
123    TT0=T0/Tref
124    RR0=TT0
125  ELSEIF(ifg==2) THEN
126   write(*,*) 'start of absc loop for h2o'
127! VIBRATIONAL PARTITION FUNCTION FOR H2O
128    V0=(1d0-DEXP(-3652d0*hckt0))*(1d0-DEXP(-1595d0*hckt0))*   &
129            (1d0-DEXP(-3756d0*hckt0))
130    V=(1d0-DEXP(-3652d0*hckt))*(1d0-DEXP(-1595d0*hckt))*   &
131            (1d0-DEXP(-3756d0*hckt))
132    VV0=V/V0
133! ROTATIONAL PARTITION FUNCTION
134    TT0=T0/Tref
135    RR0=TT0**1.5
136  ELSE
137   write(*,*) 'start of absc loop for ch4'
138! VIBRATIONAL PARTITION FUNCTION FOR CH4
139    V0=(1d0-DEXP(-1306d0*hckt0))**3*(1d0-DEXP(-1526d0*hckt0))**2*   &
140            (1d0-DEXP(-2914d0*hckt0))*(1d0-DEXP(-3020d0*hckt0))**3
141    V=(1d0-DEXP(-1306d0*hckt))**3*(1d0-DEXP(-1526d0*hckt))**2*   &
142            (1d0-DEXP(-2914d0*hckt))*(1d0-DEXP(-3020d0*hckt))**3
143    VV0=V/V0
144! ROTATIONAL PARTITION FUNCTION
145    TT0=T0/Tref
146    RR0=TT0**1.5
147  ENDIF
148   VR0=VV0*RR0
149! Scan over lines
150    c1sigt4(0)=c1/(sigma*Tref**4)
151    DO l=1,lines
152! print out every 10000 lines so that user sees machine is working
153    IF(l-(l/10000)*10000 == 0) write(*,*) 'reading line #',l
154      line_width=(Data(l,4)*xmfr(ifg)+Data(l,3)*(1.d0-xmfr(ifg)))*patm*TT0**Data(l,6)
155! molecule-based intensity
156      intensity=Data(l,2)*VR0*DEXP(hcktt0*Data(l,5))
157! pressure-based or linear intensity
158      intensity=intensity*cr*(1d0-DEXP(-data(l,1)*hckt))/(1d0-DEXP(-data(l,1)*hckt0))
159! calculate Planck-mean absorption coefficient to ensure adequacy of spectral resolution
160      kplancksum=kplancksum+intensity*c1sigt4(0)*data(l,1)**3/(EXP(hck/Tref*data(l,1))-1.d0)
161! absorption coefficient at line center
162      klmax=intensity/(PI*line_width)
163      IF(klmax<klmin) CYCLE
164! Find wavenumber (subscript) closest to line center
165      icl=(data(l,1)-wvnm_b)/wvnmst+1
166! Scan over adjacent wavenumbers to see whether line makes contribution to absco
167        DO i=icl,number
168          deta=(data(l,1)-wvnm(i))/line_width
169          dk=klmax/(1.d0+deta*deta)
170          absc(i)=absc(i)+dk
171          IF(dk<klmin) EXIT
172        ENDDO
173        DO i=icl-1,1,-1
174          deta=(data(l,1)-wvnm(i))/line_width
175          dk=klmax/(1.d0+deta*deta)
176          absc(i)=absc(i)+dk
177          IF(dk<klmin) EXIT
178        ENDDO
179    ENDDO ! l
180  ENDDO ! ifg
181   write(*,*) 'end of absc loop'
182! WRITE ABSC TO DATA FILE IF DESIRED
183  IF(iwr==1) THEN
184      write(9,96)
185      write(9,98) number
186      write(9,95) wvnm_b,wvnm_e,wvnmst
187      write(9,99) (absc(i),i=1,number)
188      CLOSE(9)
189      ENDIF
19092 FORMAT(a1,3f12.5)
19195 FORMAT('#',3f12.5)
19296 FORMAT('variables = "absco"')
19397 FORMAT(a25)
19498 FORMAT('zone i=',i8)
19599 FORMAT(e14.5)
196   ENDIF ! iwr
197! ADD SOOT CONTRIBUTION TO ABSCO (AFTER STORING/RETRIEVING GAS-ONLY VALUES)
198  sootfc=sootf(fvsoot,nsoot,ksoot)
199  DO i=1,number
200    absc(i)=absc(i)+sootfc*wvnm(i)
201  ENDDO
202! CALCULATE K-DISTRIBUTION
203! pwrK_MAX AND pwrK_MIN DEFINES THE MAXIMUM AND MINIMUM pwr(K) VALUES.
204! pwrK_STEP IS THE pwr(K) INTERVAL WE USE TO SCAN THE SPECTRUM.
205! pwrK(I): pwr(K) VALUES
206! FF(I): F*DELTAK
207! GG(I): G FUNCTION (CUMULATIVE K-DISTRIBUTION FUNCTION)
208! Find maximum kappa value
209    DO i=1,number
210      IF(absc(i)>kmax) then
211        kmax=absc(i)
212        imax=i
213      endif
214    ENDDO
215    IF(ipl==0) THEN
216        write(*,93) kmax,wvnm(imax)    !////在这个地方出现了断点的问题。
217    ELSE
218        write(*,94) kmax,wvnm(imax)
219    ENDIF
22093 FORMAT(' kmax =',f12.5,'cm-1 at ',f10.4,'cm-1')
22194 FORMAT(' kmax =',f12.5,'bar-1cm-1 at ',f10.4,'cm-1')
222   IF(kmax>kdmax .and. kdmax>0.d0) &
223        write(*,*) 'WARNING!!! kdmax is less than maximum absorption coefficient!!'
224    kmin=max(1.d-9,kdmin)
225    kmax=max(kdmax,kmax)
226    pwrk_min=kmin**pwr
227    pwrk_max=kmax**pwr
228    pwrk_step=(pwrk_max-pwrk_min)/REAL(n_pwrk-1)
229!*****************************************************************************