|
[Fortran] 纯文本查看 复制代码 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] 纯文本查看 复制代码 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
|
|