SUBROUTINE XAJMX(   FILELEN    ,&
                    N		,&		! 单元出流数组大小								//输入变量
					M		,&		! 无因次单位数组大小							//输入变量
					WUM		,&		! 上层张力水容量wum							//输入变量
                    W1		,&		! 下层张力水容量wl							//输入变量
                    WDM		,&		! 深层张力水容量wdm							//输入变量
                    KC	    ,&		!蒸发能力折算系数KC							//输入变量
                    C	    ,&		!深层蒸发系数c							    //输入变量
                    B	    ,&		!张力水蓄水容量系数b						    //输入变量
                    IMP1    ,&		!不透水面积比率imp1							//输入变量
                    SM	    ,&		!自由水蓄水容量sm								//输入变量
                    EX	    ,&		!自由水蓄水容量指数ex							//输入变量
                    KG	    ,&		!地下水出流系数kg								//输入变量
                    KSS	    ,&		!壤中流出流系数kss								//输入变量
                    KKG	    ,&		!地下水出流系数kkg								//输入变量
                    KKSS    ,&		!壤中流出流系数kkss								//输入变量
					AREA	,&		! 单元面积										//输入变量
					UH		,&		! 无因次单位线									//输入变量
					DT		,&		! 时段步长										//输入变量
					P		,&		! 降雨系列										//输入变量
					EP		,&		! 蒸发皿蒸发能力								//输入变量
					W		,&		! 土壤含水层	1.上层	2.下层	3.深层			//输入变量
					FR		,&		! 初始产流面积									//输入变量
					S		,&		! 初始自由水深									//输入变量
					QRSS0	,&		! 初始壤中流流量								//输入变量
					QRG0	,&		! 初始地下水径流量								//输入变量
					QR		)		! 单元出流										//输出变量

	IMPLICIT NONE

!///////////////////////////////////////变量声明//////////////////////////////////////////////
    INTEGER::N					! 
	INTEGER::M					!
	REAL::PAR(13)				! 1.上层张力水容量wum	2.下层张力水容量wl	3.深层张力水容量wdm
								! 4.蒸发能力折算系数KC.深层蒸发系数c		6.张力水蓄水容量系数b
								! 7.不透水面积比率imp1	8.自由水蓄水容量sm	9.自由水蓄水容量指数ex
								!10.地下水出流系数kg 11.壤中流出流系数kss	12.地下水出流系数kkg
								!13.壤中流出流系数kkss
    
    REAL::WUM                    !上层张力水容量wum	
 	REAL::W1                     !下层张力水容量wl	
 	REAL::WDM                   ! 深层张力水容量wdm
	REAL::KC					! 蒸发能力折算系数
	REAL::C						! 深层蒸发系数	
	REAL::B						! 张力水蓄水容量系数
 	REAL::IMP1					! 不透水面积比率		
	REAL::SM					! 自由水蓄水容量
	REAL::EX					! 自由水蓄水容量指数
	REAL::KG					! 地下水出流系数
	REAL::KSS					! 壤中流出流系数
	REAL::KKG					! 地下水出流系数	
	REAL::KKSS					! 壤中流出流系数		
    
	REAL::AREA					! 单元面积
	REAL::UH(M)					! 无因次单位线
	REAL::DT					! 时段步长
	REAL::P(N)					! 降雨系列
	REAL::EP(N)					! 蒸发皿蒸发能力
	REAL::QR(N)					! 单元出流
	REAL::W(3)					! 土壤含水层	1.上层	2.下层	3.深层
	REAL::FR					! 初始产流面积
	REAL::S						! 初始自由水深
	REAL::QRSS0					! 初始壤中流流量
	REAL::QRG0					! 初始地下水径流量

	INTEGER::D
    REAL::KSSD
    REAL::KGD 
 	REAL::E(3) 
 	REAL::WM(3)	

!以下变量是原vb程序中未声明的变量
	INTEGER::I					! 计数器					//临时变量
	INTEGER::J					! 计数器					//临时变量
	INTEGER::ICHECK				! 判断计算时段长度是否合适的识别码
	INTEGER::NN

	REAL(KIND=8),PARAMETER::C5=5.000000000

	REAL::U
	REAL::CI
	REAL::CG
	REAL::WM0
	REAL::W0
	REAL::PE
	REAL::R
	REAL::RIMP
	REAL::WMM
	REAL::A
	REAL::X
	REAL::RS
	REAL::RSS
	REAL::RGD
	REAL::RG
	REAL::SS
	REAL::Q
	REAL::KSSDD
	REAL::KGDD
	REAL::SMM
	REAL::SMMF
	REAL::SMF
	REAL::AU
	REAL::RSD
	REAL::RSSD
	REAL::QRS
	REAL::QRSS
	REAL::QRG
	REAL::QTR
    
    INTEGER:: FILELEN
    
    REAL::WTEMP(3)
    REAL::FRTEMP	
    REAL::STEMP	
    REAL::QRSS0TEMP		
    REAL::QRG0TEMP		
    
    
!///////////////////////////////////////变量声明//////////////////////////////////////////////
    DO I = 1,3
        WTEMP(I) = W(I)
    END DO
    FRTEMP = FR
    STEMP =S
    QRSS0TEMP=QRSS0
    QRG0TEMP=QRG0
    
!///////////////////////////////////////计算区域//////////////////////////////////////////////

    D = 0
    KSSD = 0.0
    KGD = 0.0 
    do i = 1,3
        E(I)=0.0
    END DO

    U        = 0.0
    CI       = 0.0
    CG       = 0.0
    WM0      = 0.0
    W0       = 0.0
    PE       = 0.0
    R        = 0.0
    RIMP     = 0.0
    WMM      = 0.0
    A        = 0.0
    X        = 0.0  
    RS       = 0.0
    RSS      = 0.0
    RGD      = 0.0
    RG       = 0.0
    SS       = 0.0
    Q        = 0.0
    KSSDD    = 0.0 
    KGDD     = 0.0
    SMM      = 0.0
    SMMF     = 0.0   
    SMF      = 0.0
    AU       = 0.0
    RSD      = 0.0
    RSSD     = 0.0
    QRS      = 0.0
    QRSS     = 0.0
    QRG      = 0.0
    QTR      = 0.0

    
	! 赋值
	ICHECK = 1

	DO I = 1, N 

		QR(I) = 0.0	

	END DO

	U = AREA / 3.6 / DT

	IF( DT.LE.24.0 )THEN

        D = 24 / DT
        CI = KKSS ** (1.0 / D)
        CG = KKG ** (1.0 / D)
        KSSD = (1.0 - (1.0 - (KG + KSS)) ** (1.0 / D)) / (1.0 + &
            KG / KSS)
        KGD = KSSD * KG / KSS
    
	ELSE
		ICHECK = 0
    END IF

   DO I = 1 , N

        IF (EP(I) .LT. 0.0)THEN
			EP(I) = 0.0
		END IF

        IF (P(I) .LT. 0.0)THEN
			P(I) = 0.0
		END IF

        EP(I) = EP(I) * KC
        WM0 = WUM + W1 + WDM
        W0 = W(1) + W(2) + W(3)
        PE = P(I) - EP(I)
 
		! 赋初值
        R = 0.0
        RIMP = 0.0

        IF(PE.GT.0.0)THEN

            WMM = (1.0 + B) * WM0 / (1.0 - IMP1)

            IF ((WM0 - W0).LE.0.0001)THEN

                A = WMM

            ELSE

                A = WMM * (1.0 - (1.0 - W0 / WM0) ** (1.0 / (1.0 + B)))

            END IF

            IF ((PE + A) .LT. WMM) THEN

                R = PE - WM0 + W0 + WM0 * ((1.0 - (PE + A) / WMM) &
                    ** (1.0 + B))

            ELSE

                R = PE - (WM0 - W0)

            END IF

            RIMP = PE * IMP1

        END IF

        IF ((W(1) + P(I)) .GT. EP(I)) THEN

            E(1) = EP(I)
            E(2) = 0.0
            E(3) = 0.0

        ELSE

            E(1) = W(1) + P(I)
            E(2) = (EP(I) - E(1)) * W(2) / W1

            IF (W(2) .LE.( C * WM(2))) THEN

                E(2) = C * (EP(I) - E(1))
                E(3) = 0.0

                IF (W(2) .GE. (C * (EP(I) - E(1)))) THEN

                    E(2) = C * (EP(I) - E(1))
                    E(3) = 0.0

                ELSE
                    E(2) = W(2)
                    E(3) = C * (EP(I) - E(1) - E(2))

                END IF

            END IF

        END IF
        
        W(1) = W(1) + P(I) - R - E(1)
        W(2) = W(2) - E(2)
        W(3) = W(3) - E(3)

        IF (W(1) .GT. WUM) THEN

            W(2) = W(1) - WUM + W(2)
            W(1) = WUM

            IF (W(2) .GT. W1) THEN

                W(3) = W(3) + W(2) - W1
                W(2) = W1

            END IF

        END IF

        X = FR

        IF (PE .LE.0.0) THEN

            RS = 0.0
            RSS = S * KSSD * FR
            RGD = S * FR * KGD
            S = S - (RSS + RG) / FR

        ELSE

            FR = R / PE
            S = X * S / FR
            SS = S
            Q = R / FR
            NN = INT(Q / C5) + 1		! 在vb程序中C5是双精度常数5
            Q = Q / NN
            KSSDD = (1.0 - (1.0 - (KGD + KSSD)) ** (1.0 / NN)) /&
            (1.0 + KGD / KSSD)
            KGDD = KSSDD * KGD / KSSD
            RS = 0.0
            RSS = 0.0
            RG = 0.0
            SMM = (1 + EX) * SM

            IF (EX .LT. 0.001) THEN

                SMMF = SMM

            ELSE

                SMMF = SMM * (1.0 - (1.0 - FR) ** (1.0 / EX))

            END IF

            SMF = SMMF / (1.0 + EX)

            DO J = 1 , NN

                IF (S .GT. SMF) THEN

					S = SMF

				END IF

                AU = SMMF * (1.0 - (1.0 - S / SMF) ** (1.0 /&
                    (1.0 + EX)))

                IF ((Q + AU) .LE. 0.0) THEN

                    RSD = 0.0
                    RSSD = 0.0
                    RGD = 0.0
                    S = 0.0

                ELSE IF ((Q + AU) .GE. SMMF) THEN

                    RSD = (Q + S - SMF) * FR
                    RSSD = SMF * KSSDD * FR
                    RGD = SMF * FR * KGDD
                    S = SMF - (KSSD + KGD) / FR

                ELSE IF ((Q + AU) .LT. SMMF) THEN

                    RSD = (Q - SMF + S + SMF * (1.0 - (Q + AU)&
                        / SMMF) ** (1.0 + EX)) * FR
                    RSSD = (S + Q - RSD / FR) * KSSDD * FR
                    RGD = (S + Q - RSD / FR) * KGDD * FR
                    S = S + Q - (RSD + RSSD + RGD) / FR

                END IF

				RS = RS + RSD
				RSS = RSS + RSSD
                RG = RG + RGD

            END DO

        END IF

        RS = RS * (1.0 - IMP1)
        RSS = RSS * (1.0 - IMP1)
        RG = RG * (1.0 - IMP1)
        QRS = (RS + RIMP) * U
        QRSS = QRSS0 * CI + RSS * (1.0 - CI) * U
        QRG = QRG0 * CG + RG * (1.0 - CG) * U
        QTR = QRS + QRSS + QRG

        DO J = 1 , M

            IF ((I + J - 1) .LE. N) THEN

                QR(I + J - 1) = QR(I + J - 1) + QTR * UH(J)

            END IF

        END DO

        QRSS0 = QRSS
        QRG0 = QRG

   END DO
   

    DO I = 1,3
        W(I) = WTEMP(I)
    END DO
    FR = FRTEMP
    S =STEMP
    QRSS0=QRSS0TEMP
    QRG0=QRG0TEMP
    
    
END SUBROUTINE XAJMX