|  | 
 
| [Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode C     check for water balance error
      BALERD = STORD - (INTR + INTS + INTRMS + SNOW + SWAT + GWAT) + 
     *   PRECD - EVAPD - FLOWD - SEEPD
      IF (ABS(BALERD) .GT. .003) THEN
         PRINT*, 'Serious water balance error (BALERD .GT. 0.003 mm)'
         PRINT*, 'Must be a programming error'
         write(*,*)'intrms=',INTRMS
         WRITE(*,*)'EVAPD=',EVAPD
         WRITE(*,*)'IRVPMS=',IRVPMSD,IRVPMS
         CALL INTERMS (RNETMS, GIVP, MAI, FRINTLMS, CINTRLMS, 
     *   DTP, INTRMS, RINTMS, IRVPMS)
         STOP
      END IF我改了程序后这个水量平衡的等式就不成立了,然后逐个输出值找错误最后确定就是IRVPMS值出错,而他的出错主要是子程序EVP我新加的CATCHMS初始值不对,可是我后来给他附了初值可是也不太管用
 附件里B90v3是原版程序,B90v4是我改过的
 哪位高手大神老师在北京啊?我想现场找你们问问,正好你们省的看程序了。
 以下是我改的EVP  子程序
 [Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode C     *************************************************************************
      SUBROUTINE INTER (RFAL, PINT, LAI, SAI, FRINTL, FRINTS, CINTRL, 
     *   CINTRS, DTP, INTR, RINT, IRVP)
C     rain interception, used when NPINT > 1
C     same routine is used for snow interception, with different calling 
C        variables
      IMPLICIT NONE
C     input
        REAL   RFAL     ! rainfall rate, mm/d
        REAL   PINT     ! potential interception rate, mm/d
        REAL   LAI      ! projected leaf area index, m2/m2
        REAL   SAI      ! projected stem area index, m2/m2
        REAL   FRINTL   ! intercepted fraction of RFAL per unit LAI
        REAL   FRINTS   ! intercepted fraction of RFAL per unit SAI
        REAL   CINTRL   ! maximum interception storage of rain per unit LAI, mm
        REAL   CINTRS   ! maximum interception storage of rain per unit SAI, mm
        REAL   DTP      ! precipitation interval time step, d
        REAL   INTR     ! intercepted rain, mm
C     output
        REAL   RINT     ! rain catch rate, mm/d
        REAL   IRVP     ! evaporation rate of intercepted rain, mm/d
C     local
        REAL   INTRMX   ! maximum canopy storage for rain, mm
        REAL   CATCH    ! maximum RINT, mm/d
        REAL   NEWINT   ! first approximation to new canopy storage (INTR)
C
      CATCH = (FRINTL * LAI + FRINTS * SAI) * RFAL
      INTRMX = CINTRL * LAI + CINTRS * SAI
      NEWINT = INTR + (CATCH - PINT) * DTP
      IF (NEWINT .GT. 0.) THEN
C        canopy is wet throughout DTP
         IRVP = PINT
         IF (NEWINT .GT. INTRMX) THEN
C           canopy capacity is reached
            RINT = PINT + (INTRMX - INTR) / DTP
C           RINT can be negative if INTR exists and LAI or SAI is 
C              decreasing over time
         ELSE
C           canopy capacity is not reached
            RINT = CATCH
         END IF
      ELSE
C        canopy dries during interval or stays dry
         RINT = CATCH
         IRVP = (INTR / DTP) + CATCH
C        IRVP is < PINT
      END IF
      END
C     *************************************************************************
      SUBROUTINE INTER24 (RFAL, PINT, LAI, SAI, FRINTL, FRINTS, CINTRL, 
     *   CINTRS, DURATN, INTR, RINT, IRVP)
C     rain interception with duration in hours, used when NPINT = 1
C     same routine is used for snow interception, with different calling 
C        variables
      IMPLICIT NONE
C     input
        REAL   RFAL     ! 24-hour average rainfall rate, mm/d
        REAL   PINT     ! potential interception rate, mm/d
        REAL   LAI      ! projected leaf area index, m2/m2
        REAL   SAI      ! projected stem area index, m2/m2
        REAL   FRINTL   ! intercepted fraction of RFAL per unit LAI
        REAL   FRINTS   ! intercepted fraction of RFAL per unit SAI
        REAL   CINTRL   ! maximum interception storage of rain per unit LAI, mm
        REAL   CINTRS   ! maximum interception storage of rain per unit SAI, mm
        REAL   DURATN   ! average storm duration, hr
        REAL   INTR     ! intercepted rain storage, mm,
C     output
        REAL   RINT     ! rain catch rate, mm/d
        REAL   IRVP     ! evaporation rate of intercepted rain, mm/d
C     local
        REAL   INTRMX   ! maximum canopy storage for rain, mm
        REAL   INTRNU   ! canopy storage at end of hour, mm
        REAL   NEWINT   ! first approximation to INTRNU, mm
        REAL   RINTHR   ! rain catch rate for hour, mm/hr
        REAL   CATCH    ! maximum RINTHR, mm/hr
        REAL   IRVPHR   ! evaporation rate for hour, mm/hr
        REAL   SMINT    ! daily accumulated actual catch, mm
        REAL   SMVP     ! daily accumulated actual evaporation, mm
        INTEGER IHD     ! half DURATN in truncated integer hours
        INTEGER I       ! hour, 0 to 23
        REAL   DTH      ! time step, = 1 hr
C     intrinsic
C       REAL, INT
C
      IHD = INT((DURATN + .1) / 2)
      INTRMX = CINTRL * LAI + CINTRS * SAI
      INTRNU = INTR
      SMINT = 0.
      SMVP = 0.
      DTH = 1.
      DO 100 I = 0, 23
         IF (I .LT. (12 - IHD) .OR. I .GE. (12 + IHD)) THEN
C           before or after rain
            CATCH = 0.
         ELSE
C           during rain, mm/hr is rate in mm/d divided by hr of rain/d
            CATCH = (FRINTL * LAI + FRINTS * SAI) * RFAL / REAL(2 * IHD)
         ENDIF
         NEWINT = INTRNU + (CATCH - PINT / 24.) * DTH
         IF (NEWINT .GT. .0001) THEN
C           canopy is wet throughout hour, evap rate is PINT
            IRVPHR = PINT / 24.
            IF (NEWINT .GT. INTRMX) THEN
C              canopy capacity is reached
               RINTHR = IRVPHR + (INTRMX - INTRNU) / DTH
C              INTRMX - INTRNU can be negative if LAI or SAI is decreasing 
C                 over time
            ELSE
C              canopy capacity is not reached
               RINTHR = CATCH
            END IF
         ELSE
C           canopy dries during hour or stays dry
            RINTHR = CATCH
            IRVPHR = INTRNU / DTH + CATCH
C           IRVPHR for hour is < PI/24
         END IF
         INTRNU = INTRNU + (RINTHR - IRVPHR) * DTH
         SMVP = SMVP + IRVPHR * DTH
         SMINT = SMINT + RINTHR * DTH
  100 CONTINUE
      IRVP = SMVP
C                 / 1 d
      RINT = SMINT
C                 / 1 d
      END
    
C     *************************************************************************
      SUBROUTINE INTERMS (RNETMS, GIVP, MAI , FRINTLMS, CINTRLMS, 
     *   DTP, INTRMS, RINTMS, IRVPMS)
C     moss rain interception, used when NPINT > 1
C     same routine is used for snow interception, with different calling 
C        variables
      IMPLICIT NONE
C     input
        REAL   RNETMS    ! throughfall  minus  RSNO, mm/d
        REAL   GIVP      ! ground evaporation with interception, mm/d
        REAL   MAI       ! projected moss area index, m2/m2
        REAL   FRINTLMS  ! intercepted fraction of RNETMS per unit MAI
        REAL   CINTRLMS  ! maximum interception storage of rain per unit MAI, mm
        REAL   DTP       ! precipitation interval time step, d
        REAL   INTRMS    ! moss intercepted rain, mm
C     output
        REAL   RINTMS    ! moss catch rain rate, mm/d        
        REAL   IRVPMS    ! evaporation rate of  moss  intercepted rain, mm/d
C     local
        REAL   INTRMSMX  ! maximum moss storage for rain, mm
        REAL   CATCHMS   ! maximum RINTMS, mm/d
        REAL   NEWINTMS  ! first approximation to new moss storage (INTR)
C
         WRITE(*,*)"1CATCHMS=",CATCHMS
      CATCHMS = FRINTLMS* MAI* RNETMS
         WRITE(*,*)"2CATCHMS=",CATCHMS
      INTRMSMX = CINTRLMS* MAI
      NEWINTMS = INTRMS + (CATCHMS - GIVP) * DTP
      IF (NEWINTMS .GT. 0.) THEN
C        moss is wet throughout DTP
         IRVPMS = GIVP
         IF (NEWINTMS .GT. INTRMSMX) THEN
C           moss capacity is reached
            RINTMS = GIVP + (INTRMSMX - INTRMS) / DTP
C           RINTMS can be negative if INTRMS exists and MAI is 
C              decreasing over time
            ELSE
C           moss capacity is not reached
            WRITE(*,*)"3CATCHMS=",CATCHMS
            RINTMS = CATCHMS
            WRITE(*,*)"1RINTMS=",RINTMS
         END IF
      ELSE
C        canopy dries during interval or stays dry
          WRITE(*,*)"4CATCHMS=",CATCHMS
         RINTMS = CATCHMS
         WRITE(*,*)"2RINTMS=",RINTMS
         IRVPMS = (INTRMS / DTP) + CATCHMS
         WRITE(*,*)"IRVPMS=",IRVPMS
C        IRVPMSis < GIVP
      END IF
      END
 | 
 |