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

kellsi

做作物生长模型的模型小白，在原有的模型基础上添加温度模块，调用原有lib读取数据文件时出错，但是不知道为什么会这样，求教！

kellsi

   FSE 2.1: Initialize model
SUCCESS

ERROR in FOPENG: Unit number is < 10 or > 99 !

Arguments of the CALL to FOPENG leading to this error:
   Unit             =     1
   File name        =                                                                                                   
   File status      =    RDO
   File type        =    FS
   Record length    =    0
   Delete privilege =
Fatal execution error, press <Enter>

胡文刚

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

kellsi

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

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

胡文刚

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

胡文刚

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

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

kellsi

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

确实是用的CVF...

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: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