Fortran Coder

查看: 66|回复: 3

[求助] 关于断点的问题,希望有大佬能看下这个代码,稍微有点长

[复制链接]

2

帖子

1

主题

0

精华

新人

F 币
33 元
贡献
12 点
发表于 2020-6-25 18:42:59 | 显示全部楼层 |阅读模式
[Fortran] 纯文本查看 复制代码
!-----------------------------------------------------------------------------------------------------
!------------------------------J.Cai added filenames start--------------------------------------------
!-----------------------------------------------------------------------------------------------------
   character(256),parameter :: kvsgFile='11kvsg-4.9-5.0-1-H2O-0.0p-CO2-1p-300K-11.dat' !
   character(256),parameter :: kvsgqFile='11kvsgq-4.9-5.0-1-H2O-0.0p-CO2-1p-300K.dat' !
   character(256),parameter :: abscFile='11absc-4.9-5.0-1-H2O-0.0p-CO2-1p-300K.dat' !
   character(256),parameter :: CO2File='02-3.7-3.8.par' 
   character(256),parameter :: CH4File='F:\JavaHAWKS\HITRAN04\02_hit2010.par' 
   character(256),parameter :: H2OFile='01_0-30000_HITEMP-35.par' 
   !character(256),parameter :: H2OFile='/opt/HITRAN/hitran96/by_molec/01_hit96.par' 
! J.Cai added filenames end
! J.Cai comment out starts
!-----------------------------------------------------------------------------------------------------
!! Selection of g-values for numerical quadrature, using a Numerical Recipes routine
!! If Numerical Recipes is not available, set nq=12, comment out the following 8 lines of code,
!! and uncomment the 5-line REAL declaration following it
!   REAL              :: gqs(nq),wqs(nq),kq(nq),aq(numt,nq),gq(nq),wq(nq),gaujac,alf=3.,bet=-.9,wsum
!! Get quadrature coefficients from Numerical recipies
!    wsum=0.
!    CALL GAUJAC(gqs,wqs,nq,alf,bet)
!          do iq=1,nq
!            gq(iq)=0.5*(1.-gqs(iq))
!            wq(iq)=wqs(iq)/(2.**(alf+bet+1)*gq(iq)**alf*(1.-gq(iq))**bet)
!            wsum=wsum+wq(iq)
!          enddo
!! Correction to make sum(wq)=1
!         wq=wq/wsum
!! End quadrature coefficients from Numerical recipies
! J.Cai comment out ends
! Selection of precalculated g-values for numerical quadrature, for nq=12,alf=3.,bet=0.
! J.Cai uncomment out starts
    REAL              :: kq(nq),aq(numt,nq), &
           gq(nq)=(/ 5.120075E-02,1.170678E-01,2.015873E-01,3.007074E-01,4.095012E-01,5.225285E-01,  &
                     6.341280E-01,7.387071E-01,8.310236E-01,9.064499E-01,9.612060E-01,9.925594E-01/),&
           wq(nq)=(/ 5.556622E-02,7.576839E-02,9.258290E-02,1.048306E-01,1.118451E-01,1.132605E-01,  &
                     1.090012E-01,9.927844E-02,8.457905E-02,6.563999E-02,4.341329E-02,1.904792E-02/)
! J.Cai uncomment out starts
! PRESSURE IN atm
   patm=P/1.01325
! IF P-BASED absc is requested make sure it is not a mixture
    IF(ipl==1) THEN
        cnchk=xmfr(1)*xmfr(2)+xmfr(1)*xmfr(3)+xmfr(3)*xmfr(2)
        IF(cnchk > 1.d-6) PAUSE 'ERROR: you cannot use pressure-based absorption coefficient for mixture!'
    ENDIF
! Open output files
! Output file for k,g(T,k),a(T,k) for all numT temperatures and all n_pwrk k-values
   OPEN(7,FILE=kvsgFile) ! J.Cai centralize filename
! Header formatted for TECPLOT, for a numT of 4
   write(7,6)
6 format('VARIABLES = k,g0,g1,g2,g3,g4,a1,a2,a3,a4')
! Output file for wq,gq(Tref),kq(gq),aq(T,kq) for all numT temperatures and all nq g-values
   IF(ipr==2) THEN
       OPEN(8,FILE=kvsgqFile,STATUS='unknown') ! J.Cai centralize filename
! Header formatted for readability, for a numT of 4
        write(8,8)
   ENDIF
8 format('    wq',9x,'gq',9x,'kq',8x,'aq1',8x,'aq2',8x,'aq3',8x,'aq4')
! File containing absorption coefficient
   !IF(iwr>0) OPEN(9,FILE='c:\absco\absch2o-0p-2000K.dat',STATUS='unknown')
   IF(iwr>0) OPEN(9,FILE=abscFile,STATUS='unknown') ! J.Cai centralize filename
   hck=hhh*ccc/kkk
   hckt0=hck/T0
   IF(iwr==2) THEN
! READ (GAS-ONLY) ABSORPTION COEFFICIENT FROM FILE
      read(9,97) dummy
      read(9,98) number
      read(9,92) dummy,wvnm_b2,wvnm_e2,wvnmst2
      number2=(wvnm_e2-wvnm_b2)/wvnmst2+1
      IF(number2 /= number) PAUSE 'bad data file'
      read(9,*) (absc(i),i=1,number)
      CLOSE(9)
      DO i=1,number
         wvnm(i)=wvnm_b2+(i-1)*wvnmst2
      ENDDO
   ELSE
! CALCULATE NECESSARY WAVENUMBERS AND ABSORPTION COEFFICIENTS
! WVNM(I): WAVENUMBER
! ABSC(I): ABSORPTION COEFFICIENT; INITIALIZE FOR SOOT
   number=(wvnm_e-wvnm_b)/wvnmst+1
   IF(number>nabs) PAUSE 'increase nabs'
   DO i=1,number
    wvnm(i)=wvnm_b+(i-1)*wvnmst
    absc(i)=0.d0
   ENDDO
  kplancksum=0.d0
! Scan over gases
  DO ifg=1,3
  IF(xmfr(ifg)<.9d-3) CYCLE
! Open HITEMP database
  IF(ifg==1) THEN
     lu=11
     write(*,*) 'Reading CO2 data'
     OPEN(lu,FILE=CO2File, action='READ')  ! J.Cai change path
  ELSEIF(ifg==2) THEN
     lu=12
     write(*,*) 'Reading H2O data'
     OPEN(lu,FILE=H2OFile, action='READ')  ! J.Cai change path
  ELSE
     lu=13
     write(*,*) 'Reading CH4 data'
     OPEN(lu,FILE=CH4File, action='READ') ! J.Cai change path
  ENDIF
! GET LINE INFORMATION FROM HITEMP (all lines wvnm_b.le.lmbda.le.wvnm_e)
   CALL Getdata(lines)
   CLOSE(lu)
   write(*,*) 'Gas ',ifg,',  lines read: ', lines
   mass=1.d3*molemass(ifg)/nnn
   hckt=hck/Tref
   hcktt0=hckt0-hckt
! MULTIPLIER FOR PRESSURE-BASED ABSORPTION COEFFICIENT
   cr=nnn/8.314d1/Tref
! MULTIPLIER FOR LINEAR ABSORPTION COEFFICIENT
   IF(ipl==0) cr=cr*xmfr(ifg)*P
   IF(ifg==1) THEN
    write(*,*) 'start of absc loop for co2'
! VIBRATIONAL PARTITION FUNCTION FOR CO2
    V0=(1d0-DEXP(-666d0*hckt0))**2*(1d0-DEXP(-2396d0*hckt0))*   &
            (1d0-DEXP(-1351d0*hckt0))
    V=(1d0-DEXP(-666d0*hckt))**2*(1d0-DEXP(-2396d0*hckt))*   &
            (1d0-DEXP(-1351d0*hckt))
    VV0=V/V0
! ROTATIONAL PARTITION FUNCTION
    TT0=T0/Tref
    RR0=TT0
  ELSEIF(ifg==2) THEN
   write(*,*) 'start of absc loop for h2o'
! VIBRATIONAL PARTITION FUNCTION FOR H2O
    V0=(1d0-DEXP(-3652d0*hckt0))*(1d0-DEXP(-1595d0*hckt0))*   &
            (1d0-DEXP(-3756d0*hckt0))
    V=(1d0-DEXP(-3652d0*hckt))*(1d0-DEXP(-1595d0*hckt))*   &
            (1d0-DEXP(-3756d0*hckt))
    VV0=V/V0
! ROTATIONAL PARTITION FUNCTION
    TT0=T0/Tref
    RR0=TT0**1.5
  ELSE
   write(*,*) 'start of absc loop for ch4'
! VIBRATIONAL PARTITION FUNCTION FOR CH4
    V0=(1d0-DEXP(-1306d0*hckt0))**3*(1d0-DEXP(-1526d0*hckt0))**2*   &
            (1d0-DEXP(-2914d0*hckt0))*(1d0-DEXP(-3020d0*hckt0))**3
    V=(1d0-DEXP(-1306d0*hckt))**3*(1d0-DEXP(-1526d0*hckt))**2*   &
            (1d0-DEXP(-2914d0*hckt))*(1d0-DEXP(-3020d0*hckt))**3
    VV0=V/V0
! ROTATIONAL PARTITION FUNCTION
    TT0=T0/Tref
    RR0=TT0**1.5
  ENDIF
   VR0=VV0*RR0
! Scan over lines
    c1sigt4(0)=c1/(sigma*Tref**4)
    DO l=1,lines
! print out every 10000 lines so that user sees machine is working
    IF(l-(l/10000)*10000 == 0) write(*,*) 'reading line #',l
      line_width=(Data(l,4)*xmfr(ifg)+Data(l,3)*(1.d0-xmfr(ifg)))*patm*TT0**Data(l,6)
! molecule-based intensity
      intensity=Data(l,2)*VR0*DEXP(hcktt0*Data(l,5))
! pressure-based or linear intensity
      intensity=intensity*cr*(1d0-DEXP(-data(l,1)*hckt))/(1d0-DEXP(-data(l,1)*hckt0))
! calculate Planck-mean absorption coefficient to ensure adequacy of spectral resolution
      kplancksum=kplancksum+intensity*c1sigt4(0)*data(l,1)**3/(EXP(hck/Tref*data(l,1))-1.d0)
! absorption coefficient at line center
      klmax=intensity/(PI*line_width)
      IF(klmax<klmin) CYCLE
! Find wavenumber (subscript) closest to line center
      icl=(data(l,1)-wvnm_b)/wvnmst+1
! Scan over adjacent wavenumbers to see whether line makes contribution to absco
        DO i=icl,number
          deta=(data(l,1)-wvnm(i))/line_width
          dk=klmax/(1.d0+deta*deta)
          absc(i)=absc(i)+dk
          IF(dk<klmin) EXIT
        ENDDO
        DO i=icl-1,1,-1
          deta=(data(l,1)-wvnm(i))/line_width
          dk=klmax/(1.d0+deta*deta)
          absc(i)=absc(i)+dk
          IF(dk<klmin) EXIT
        ENDDO
    ENDDO ! l
  ENDDO ! ifg
   write(*,*) 'end of absc loop'
! WRITE ABSC TO DATA FILE IF DESIRED
  IF(iwr==1) THEN
      write(9,96)
      write(9,98) number
      write(9,95) wvnm_b,wvnm_e,wvnmst
      write(9,99) (absc(i),i=1,number)
      CLOSE(9)
      ENDIF
92 FORMAT(a1,3f12.5)
95 FORMAT('#',3f12.5)
96 FORMAT('variables = "absco"')
97 FORMAT(a25)
98 FORMAT('zone i=',i8)
99 FORMAT(e14.5)
   ENDIF ! iwr
! ADD SOOT CONTRIBUTION TO ABSCO (AFTER STORING/RETRIEVING GAS-ONLY VALUES)
  sootfc=sootf(fvsoot,nsoot,ksoot)
  DO i=1,number
    absc(i)=absc(i)+sootfc*wvnm(i)
  ENDDO
! CALCULATE K-DISTRIBUTION
! pwrK_MAX AND pwrK_MIN DEFINES THE MAXIMUM AND MINIMUM pwr(K) VALUES.
! pwrK_STEP IS THE pwr(K) INTERVAL WE USE TO SCAN THE SPECTRUM.
! pwrK(I): pwr(K) VALUES
! FF(I): F*DELTAK
! GG(I): G FUNCTION (CUMULATIVE K-DISTRIBUTION FUNCTION)
! Find maximum kappa value
    DO i=1,number
      IF(absc(i)>kmax) then
        kmax=absc(i)
        imax=i
      endif
    ENDDO
    IF(ipl==0) THEN
        write(*,93) kmax,wvnm(imax)    !////在这个地方出现了断点的问题。
    ELSE
        write(*,94) kmax,wvnm(imax)
    ENDIF
93 FORMAT(' kmax =',f12.5,'cm-1 at ',f10.4,'cm-1')
94 FORMAT(' kmax =',f12.5,'bar-1cm-1 at ',f10.4,'cm-1')
   IF(kmax>kdmax .and. kdmax>0.d0) &
        write(*,*) 'WARNING!!! kdmax is less than maximum absorption coefficient!!'
    kmin=max(1.d-9,kdmin)
    kmax=max(kdmax,kmax)
    pwrk_min=kmin**pwr
    pwrk_max=kmax**pwr
    pwrk_step=(pwrk_max-pwrk_min)/REAL(n_pwrk-1)
!*****************************************************************************
回复

使用道具 举报

128

帖子

2

主题

0

精华

大师

F 币
769 元
贡献
400 点

规矩勋章

发表于 2020-6-28 08:36:52 | 显示全部楼层
什么问题?
回复

使用道具 举报

2

帖子

1

主题

0

精华

新人

F 币
33 元
贡献
12 点
 楼主| 发表于 2020-6-28 09:37:04 | 显示全部楼层

出现了断点

1532

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1365 元
贡献
907 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2020-6-28 11:06:17 | 显示全部楼层
无法根据代码片段(不完整),来推测运行时错误。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-7-8 17:45

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表