在Fortran模型中加入一个自己写的子程序出现读入问题

kellsi

胡文刚 发表于 2021-3-12 08:43
你给出的部分代码看不出任何问题。

搜索你的代码，是否存在 fopeng 这样的函数？或者类似

这个通道的问题出现在从soil文件中读入数据那里，请教如何更改unit的值
您说的FOPENG在TTUTIL.LIB中封装了，无法更改这个逻辑设置

胡文刚

你给出的部分代码看不出任何问题。

搜索你的代码，是否存在 fopeng 这样的函数？或者类似
ERROR in FOPENG: Unit number is < 10 or > 99
这样的字符串？

截图时，请让我们看到你的开发环境，确定是CVF？

kellsi

胡文刚 发表于 2021-3-11 22:18
我怀疑这个错误可能不是编译器抛出的，而是你的代码逻辑，比如类似这样的语句

if( unit < 1 .or. unit > 9 ...

!  soiltemperature.f90
!
!  FUNCTIONS:
!  Console3 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM:  Soil Temperature
!
!  PURPOSE:  FOR COMPUTING SOIL TEMPERATURE AND SOIL  HEAT FLUX
!  TI:  一天中的第几个小时，h
!  DT:  时间步长，second
!  TA:  Soil annual average temperature of the first layers,表土年均温℃
!  Tamp:  Amplitude of annual temperature wave at soil surface,表土年温度波幅 ℃
!  TB:  深层的基础温度，℃
!  TN0: 下一时间步长的初始温度, ℃
!  TN:  下一时间步长的温度，Tnew, ℃
!  DA:  日序，第几天
!  T(NSL): 逐层的上边界温度，℃
!  K0:  边界层温度传导率，W/(M^2  K)
!
!  CALL SUBROUTINE: WEATHR
!****************************************************************************

      Subroutine Tsoil(DOY,TI,WCLQT,TMDA,TKL)

      USE CHART
      Implicit none


      INTEGER   I,NSL
          PARAMETER (NSL=10)
          INTEGER   ITASK,IUNITD,IUNITL
      CHARACTER (100) FILEI2
      REAL      TA,Tamp,TB,TI,DT,DA,DT1
      REAL      C1(NSL),C2(NSL),C3(NSL),C4(NSL),T0,TN0
      REAL      BD(NSL),F,G,CLAYX(NSL),K0,SOILT(NSL)
  
      REAL      W(NSL+1), T(NSL+1), TN(NSL+1), K(NSL), CP(NSL+1), A(NSL), B(NSL)
      REAL      C(NSL), D(NSL), Z(NSL+1), WCLQT(NSL), KM(NSL+1)
   
      CHARACTER (1)  DUMMY
          REAL                DOY,TKL(NSL), TMDA
          REAL      IDOY, ISTAT2, RDD, TMMN, TMMX, VP, WN, RAIN

!从soil文件中读入所需参数值
      CALL RDINIT(IUNITD,IUNITL,FILEI2)
          CALL RDFREA('BD'  ,BD ,NSL,NSL)
          CALL RDFREA('CLAYX', CLAYX ,NSL, NSL)
          CALL RDFREA('TKL', TKL ,NSL, NSL)
          CALL RDFREA('SOILT', SOILT ,NSL, NSL)
!          CALL RDSREA('TA', TA ,NSL, NSL)
!          CALL RDSREA('Tamp', Tamp ,NSL, NSL)
!------- Reading of soil data completed
      CLOSE (IUNITD)

!参数初始值赋值
      TA=20;TAMP=15;TB=20   !   BD=1.3;
      K0=20 !BOUNDARY LAYER CONDUCTANCE IN W/(M^2  K)
          !从weather文件中读取逐日温度
          CALL WEATHR (IDOY  , ISTAT2, RDD, TMMN, TMMX, VP, WN, RAIN)
          TA=(TMMX+TMMN)/2
          TAMP=TMMX-TMMN
!剖面深度及温度的初始值预设
      Z(1)=0.0
      DO I=1,NSL
         Z(I+1)=Z(I)+TKL(I) !0.005*1.5**(I-1)  !进行剖面深度的预设
         T(I)=SOILT(I)   !定义剖面T的初始值
      END DO
!          WRITE(*,*)TKL
!          PAUSE
      T(NSL+1)=TB
      TN(NSL+1)=T(NSL+1)
      T0=TB  
      !TI IS TIME OF DAY;DT IS TIME STEP (SEC);DA IS DAY NUMBER
!      TI=0;DA=0;DT1=3600
          TI=0;DA=DOY;DT1=3600      
      F=0.6
      G=1-F
      !MC=0.12  !  CLAY FRACTION

      C1(I)=0.65-0.78*BD(I)+0.6*BD(I)**2
      C2(I)=1.06*BD(I)
      C3(I)=1+2.6/(CLAYX(I)**0.5)
      C4(I)=0.3+0.1*BD(I)**2
      !WCLQT=0.3  !传递实际土壤含水量

      DT=DT1

230   TI=TI+DT/3600   


      CP(1)= (2400000*BD(1)/2.65+4180000*WCLQT(1))*(Z(2)-0.0)/(2*DT)
      K(1)=(C1(1)+C2(1)*WCLQT(1)-(C1(1)-C4(1))*EXP(-(C3(1)*WCLQT(1))**4))/(Z(2)-0)
      DO I = 2,NSL
        CP(I)=(2400000*BD(I)/2.65+4180000*WCLQT(I))*(Z(I+1)-Z(I-1))/(2*DT)
        K(I)=(C1(I)+C2(I)*WCLQT(I)-(C1(I)-C4(I))*EXP(-(C3(I)*WCLQT(I))**4))/(Z(I+1)-Z(I))
      END DO


      !计算Ki-1/2和Ki+1/2
      !KM(1)=0.5*(K(1)+K0)
      !do I = 2,NSL
      !   KM(I)=0.5*(K(I-1)+K(I))
      !enddo
      !KM(NSL+1)=K(NSL)

      do I = 1,NSL
        KM(i)=K(i)
      end do
      KM(NSL+1)=K(NSL)



!C2---Matrix coefficience
      !TN0=TA+TAMP*SIN(0.261799*(TI-6))
          TN0=TA+TAMP*SIN(0.261799*(TI-6))
      !node=1
      A(1)= 0.0
      B(1)= F*(KM(1)+K0)+CP(1)
      C(1)= -KM(1)*F
      D(1)= G*K0*T0 + (CP(1)-G*(KM(1)+K0))*T(1)+ G*KM(1)*T(2)+ F*K0*TN0

      DO I=2,NSL-1
        A(I)=-KM(I-1)*F
        C(I)=-KM(I)*F
        B(I)=F*(KM(I)+KM(I-1))+CP(I)
        D(I)=G*KM(I-1)*T(I-1)+(CP(I)-G*(KM(I)+KM(I-1)))*T(I)+G*KM(I)*T(I+1)
      END DO

      I=NSL
      A(I)=-KM(I-1)*F
      C(I)= 0
      B(I)=F*(KM(I)+KM(I-1))+CP(I)
      D(I)=G*KM(I-1)*T(I-1)+(CP(I)-G*(KM(I)+KM(I-1)))*T(I)+G*KM(I)*T(I+1) + F*KM(I)*TN(I+1)

      !WRITE(*,*)T(NSL+1),TN(NSL+1)
     
      
!C3--- 求解三对角方程
      DO I=1,NSL-1
        C(I)=C(I)/B(I)
        D(I)=D(I)/B(I)
        B(I+1)=B(I+1)-A(I+1)*C(I)
        D(I+1)=D(I+1)-A(I+1)*D(I)
      ENDDO
      TN(NSL)=D(NSL)/B(NSL)
      DO I =NSL-1,1,-1
        TN(I)=D(I)-C(I)*TN(I+1)
      ENDDO
      
      
      WRITE (*,*)'DAY =',DA
      WRITE (*,*)'HOUR =',TI
      WRITE (*,*)'DT =',DT
      WRITE (*,*)'HEAT FLUX =',K0*(G*(T0-T(1))+F*(TN0-TN(1))) !"W/M2"
      !WRITE (*,*)'DEPTH','TEMPERATURE','K(I)'
      
      T0=TN0  !存储该一时刻的上边界温度
      WRITE (*,*)'T0=',T0
      DO I=1,NSL        
        WRITE (*,*) Z(I),TN(I),KM(I)
        T(I)=TN(I)
      ENDDO


      IF(TI+DT/3600.GT.24)THEN
         DT=(24-TI)*3600
      ENDIF

      IF(ABS(TI-24).LT.10E-4) THEN
        TI=0
        DT=DT1
        DA=DA+1
      ENDIF

          RETURN
      END

kellsi

胡文刚 发表于 2021-3-11 22:18
我怀疑这个错误可能不是编译器抛出的，而是你的代码逻辑，比如类似这样的语句

if( unit < 1 .or. unit > 9 ...

如果像您说的，是代码逻辑里面存在的问题，要怎么做才能不出现这个问题啊

kellsi

胡文刚 发表于 2021-3-11 22:18
我怀疑这个错误可能不是编译器抛出的，而是你的代码逻辑，比如类似这样的语句

if( unit < 1 .or. unit > 9 ...

kellsi

胡文刚 发表于 2021-3-11 22:03
CVF 是支持 1 做为通道号的。你确定你用的CVF？能否截图看看

确实是用的CVF...

胡文刚

我怀疑这个错误可能不是编译器抛出的，而是你的代码逻辑，比如类似这样的语句

if( unit < 1 .or. unit > 99 ) then
  write(*,*) 'ERROR in FOPENG: Unit number is < 10 or > 99 !'
end if

胡文刚

CVF 是支持 1 做为通道号的。你确定你用的CVF？能否截图看看

kellsi

胡文刚 发表于 2021-3-11 18:54
你的编译器过于陈旧。只支持 10 到 99 的通道号，建议更新你的编译器。

请问这个要怎么更新啊，现在用的模型还是在CVF下运行的，有很多调用封装到lib里了，在IVF中无法运行

胡文刚

你的编译器过于陈旧。只支持 10 到 99 的通道号，建议更新你的编译器。