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

wangzb · 发表于 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)
!*****************************************************************************