text:       SUBROUTINE CALPUV9 (ISTL_)  
  C  
  C CHANGE RECORD  
  C ** SUBROUTINE CALPUV9 CALCULATES THE EXTERNAL SOLUTION FOR P, UHDYE,  
  C ** AND VHDXE, FOR FREE SURFACE FLOWS WITH PROVISIONS FOR WETTING  
  C ** AND DRYING OF CELLS  
  C  
        USE GLOBAL  
  
        REAL,SAVE,ALLOCATABLE,DIMENSION(:)::QSUMTMP  
  
        IF(.NOT.ALLOCATED(QSUMTMP))THEN
          ALLOCATE(QSUMTMP(LCM))  
          QSUMTMP=0.
        ENDIF
  C  
        IF(DEBUG)THEN
          IF(N.EQ.1.AND.ISDSOLV.EQ.1)THEN  
            OPEN(1,FILE='FUV1.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQCOEF1.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQTERM1.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='FP1.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
          ENDIF  
          IF(N.EQ.2.AND.ISDSOLV.EQ.1)THEN  
            OPEN(1,FILE='FUV2.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQCOEF2.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQTERM2.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='FP2.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
          ENDIF  
          IF(ISDSOLV.EQ.1)THEN  
            OPEN(1,FILE='FUV.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQCOEF.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQTERM.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='FP.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
          ENDIF
        ENDIF
  C  
        DELT=DT2  
        DELTD2=DT  
        IF(ISTL_.EQ.2)THEN  
          DELT=DT  
          DELTD2=0.5*DT  
        ENDIF  
        DELTI=1./DELT  
  C  
  C **  CALCULATE EXTERNAL BUOYANCY INTEGRALS AT TIME LEVEL (N)  
  C  
        IF(BSC.GT.1.E-6)THEN
          ! *** CALCULATE EXTERNAL BUOYANCY INTEGRALS AT TIME LEVEL (N)  
          CALL CALEBI  
  C
          ! *** CALCULATE EXPLICIT EXTERNAL PRESSURE GRADIENTS  
          DO L=2,LA  
            !SBX(L)=0.5*SUB(L)*DYU(L)  
            FPGXE(L)=ROLD*FPGXE(L)+RNEW*(  
       &        -SBX(L)*HU(L)*GP*((BI2(L)+BI2(LWEST(L)))*(HP(L)-HP(LWEST(L)))  
       &        +2.0*HU(L)*(BI1(L)-BI1(LWEST(L)))    
       &        +(BE(L)+BE(LWEST(L)))*(BELV(L)-BELV(LWEST(L)))))  
          ENDDO  
          DO L=2,LA  
            LS=LSC(L)  
            !SBY(L)=0.5*SVB(L)*DXV(L)  
            FPGYE(L)=ROLD*FPGYE(L)+RNEW*(  
       &         -SBY(L)*HV(L)*GP*((BI2(L)+BI2(LS))*(HP(L)-HP(LS))  
       &        +2.0*HV(L)*(BI1(L)-BI1(LS)) 
       &        +(BE(L)+BE(LS))*(BELV(L)-BELV(LS))))  
          ENDDO  
        ENDIF  
  C  
  C **  CALCULATE EXPLICIT EXTERNAL UHDYE AND VHDXE EQUATION TERMS  
  C **  HRU=SUB*HMU*DYU/DXU & HRV=SVB*HMV*DXV/DYV  
  C  
        IF(ISTL_.EQ.2)THEN  
          DO L=2,LA  
            HUTMP(L)=0.5*(HU(L)+H1U(L))  
            HVTMP(L)=0.5*(HV(L)+H1V(L))  
          ENDDO  
        ELSE  
          DO L=2,LA  
            HUTMP(L)=HU(L)  
            HVTMP(L)=HV(L)  
          ENDDO  
        ENDIF
  C  
        DO L=2,LA  
          LS=LSC(L)  
          LN=LNC(L)  
          !DXYU(L)=DXU(L)*DYU(L)  
          !DXIU(L)=1./DXU(L)  
          !HRUO(L)=SUBO(L)*DYU(L)*DXIU(L)  
          FUHDYE(L)=UHDY1E(L)  
       &      -DELTD2*SUB(L)*HRUO(L)*HUTMP(L)*(P1(L)-P1(LWEST(L)))  
       &      +SUB(L)*DELT*DXIU(L)*(DXYU(L)*(TSX1(L)-RITB1*TBX1(L))  
       &      +FCAXE(L)+FPGXE(L)-SNLT*FXE(L))  
  C
          !DXYV(L)=DXV(L)*DYV(L)  
          !DYIV(L)=1./DYV(L)  
          !HRVO(L)=SVBO(L)*DXV(L)*DYIV(L)
          FVHDXE(L)=VHDX1E(L)  
       &      -DELTD2*SVB(L)*HRVO(L)*HVTMP(L)*(P1(L)-P1(LS ))  
       &      +SVB(L)*DELT*DYIV(L)*(DXYV(L)*(TSY1(L)-RITB1*TBY1(L))  
       &      -FCAYE(L)+FPGYE(L)-SNLT*FYE(L))
        ENDDO  
  C
        IF(ISDSOLV.GE.1.AND.DEBUG)THEN  
          OPEN(1,FILE='FUV.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
          WRITE(1,1001)N,ISTL_  
          DO L=2,LA  
            WRITE(1,1001)IL(L),JL(L),UHDY1E(L),HRUO(L),HUTMP(L),P1(L),  
       &        P1(LWEST(L)),TSX1(L),TBX1(L),FCAXE(L),FPGXE(L),FXE(L)  
          ENDDO  
          CLOSE(1)  
          IF(N.EQ.1)THEN  
            OPEN(1,FILE='FUV1.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1001)IL(L),JL(L),UHDY1E(L),HRUO(L),HUTMP(L),P1(L),  
       &          P1(LWEST(L)),TSX1(L),TBX1(L),FCAXE(L),FPGXE(L),FXE(L)  
            ENDDO  
            CLOSE(1)  
          ENDIF  
          IF(N.EQ.2)THEN  
            OPEN(1,FILE='FUV2.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1001)IL(L),JL(L),UHDY1E(L),HRUO(L),HUTMP(L),P1(L),  
       &          P1(LWEST(L)),TSX1(L),TBX1(L),FCAXE(L),FPGXE(L),FXE(L)  
            ENDDO  
            CLOSE(1)  
          ENDIF  
        ENDIF  
  C  
  C **  SET IMPLICIT BOTTOM AND VEGETATION DRAG AS APPROPRIATE  
  C  
        RCX(2:LA)=1./(1.+FBODYFXI(2:LA)) !Multiple layers no vegetation
        RCY(2:LA)=1./(1.+FBODYFYI(2:LA)) 
  
        RCX(1)=0.  
        RCY(1)=0.  
        RCX(LC)=0.  
        RCY(LC)=0.  
  C  
  C * SINGLE LAYER NO VEGETATION  
  C  
        IF(KC.EQ.1.AND.ISVEG.EQ.0.AND.RITB.GT.0.)THEN  
          DO L=2,LA  
            RCX(L)=1./( 1.+FBODYFXI(L) 
       &        +RITB*DELT*HUI(L)*STBX(L)*SQRT(VU(L)*VU(L)+U(L,1)*U(L,1)))  
            RCY(L)=1./( 1.+FBODYFYI(L)
       &        +RITB*DELT*HVI(L)*STBY(L)*SQRT(UV(L)*UV(L)+V(L,1)*V(L,1)))  
            FUHDYE(L)=FUHDYE(L)*RCX(L)  
            FVHDXE(L)=FVHDXE(L)*RCY(L)  
          ENDDO  
        ENDIF  
  C  
  C * SINGLE LAYER WITH VEGETATION  
  C  
        IF(KC.EQ.1.AND.ISVEG.GE.1)THEN  
          DO L=2,LA  
            RCX(L)=1./( 1.+FBODYFXI(L)
       &        +RITB*DELT*HUI(L)*STBX(L)*SQRT(VU(L)*VU(L)+U(L,1)*U(L,1))  
       &        +DELT*FXVEGE(L) )  
            RCY(L)=1./( 1.+FBODYFYI(L)
       &        +RITB*DELT*HVI(L)*STBY(L)*SQRT(UV(L)*UV(L)+V(L,1)*V(L,1))  
       &        +DELT*FYVEGE(L) )  
            FUHDYE(L)=FUHDYE(L)*RCX(L)  
            FVHDXE(L)=FVHDXE(L)*RCY(L)  
          ENDDO  
        ENDIF  
  C  
  C * MULTIPLE LAYERS WITH VEGETATION  
  C  
        IF(KC.GT.1.AND.ISVEG.GE.1)THEN  
          DO L=2,LA  
            RCX(L)=1./( 1.+FBODYFXI(L)+DELT*FXVEGE(L) )  
            RCY(L)=1./( 1.+FBODYFYI(L)+DELT*FYVEGE(L) )  
            FUHDYE(L)=FUHDYE(L)*RCX(L)  
            FVHDXE(L)=FVHDXE(L)*RCY(L)  
          ENDDO  
        ENDIF  
  C
  C **  ADJUST VOLUME SOURCE AND SINKS  
  C  
        IF(ISGWIE.EQ.0)THEN  
          DO L=2,LA  
            IF(QSUME(L).LE.0.)THEN  
              IF(H1P(L).LE.HDRY)THEN  
                QSUMTMP(L)=0.  
              ELSE  
                QSUMTMP(L)=-(H1P(L)-HDRY)*DXYP(L)*DELTI  
                QSUMTMP(L)=MAX(QSUMTMP(L),QSUME(L))  
              ENDIF  
            ELSE  
              QSUMTMP(L)=QSUME(L)  
            ENDIF  
          ENDDO  
          DO L=2,LA  
            DIFQVOL=QSUME(L)-QSUMTMP(L)  
            DO K=1,KC  
              QSUM(L,K)=QSUM(L,K)-DIFQVOL*DZC(K)  
            ENDDO  
            QSUME(L)=QSUMTMP(L)  
          ENDDO  
        ENDIF  
  C  
  C **  ADJUST SOURCES AND SINKS ESTIMATING SURFACE AND GROUNDWATER  
  C **  AVAILABLE FOR EVAPOTRANSPIRATON AND INFILTRATION  
  C  
        IF(ISGWIE.GE.1)THEN  
          DO L=2,LA  
            RIFTR(L)=0.  
            EVAPSW(L)=0.  
            EVAPGW(L)=0.  
            IF(H1P(L).GT.HDRY)THEN  
  C  
  C       APPLY MAXIMUM ET  
  C  
              IF(EVAPCVT.LT.0.)THEN  
                SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/  
       &            (1.+0.00412*TEM(L,KC))))  
                EVAPT(L)=CLEVAP(L)*0.7464E-3*WINDST(L)*(SVPW-VPA(L))
       &            /PATMT(L)  
              ENDIF  
              EVAPSW(L)=EVAPT(L)*DXYP(L)  
              RIFTR(L)=0.  
  C  
  C       CALCULATE DEPTH OF ACTIVE GROUNDWATER ELEV BELOW SURFACE  
  C  
              DTAGW=BELV(L)-AGWELV(L)  
              IF(DTAGW.GT.0.0)THEN  
  C  
  C         INFLITRATION CAN OCCUR, CALCULATE LIMITING RATE TO BRING  
  C         GW ELEV TO SOIL SURFACE  
  C  
                RIFTRL=RNPOR*DTAGW*DELTI  
  C  
  C         SET RIFTRL TO MIN OF LIMITING RATE OR ACTUAL RATE  
  C  
                RIFTRL=MIN(RIFTRM,RIFTRL)  
  C  
  C         ESTIMATE RATE BASED ON AVAILABLE SURFACE WATER  
  C  
                RAVAIL=(H1P(L)-HDRY)*DELTI-EVAPT(L)  
  C  
  C         SET RIFTRL TO MIN OF AVAILABLE RATE OR LIMITING RATE  
  C  
                RIFTRL=MIN(RAVAIL,RIFTRL)  
  C  
  C         CONVERT TO VOLUME FLOW UNITS  
  C  
                RIFTR(L)=RIFTRL*DXYP(L)  
              ENDIF  
  C  
  C       ADJUST VOLUME OUTFLOWS OF WET CELLS  
  C  
              IF(QSUME(L).LT.0.0)THEN  
                QSUMIET=RIFTR(L)+EVAPSW(L)  
                QEAVAIL=DXYP(L)*(H1P(L)-HDRY)*DELTI-QSUMIET  
                QEAVAIL=MAX(QEAVAIL,0.0)  
                QEAVAIL=-QEAVAIL  
                QSUMTMP(L)=MAX(QSUME(L),QEAVAIL)  
              ELSE  
                QSUMTMP(L)=QSUME(L)  
              ENDIF  
            ELSE  
              RIFTR(L)=0.  
              EVAPSW(L)=0.  
              QSUMTMP(L)=MAX(QSUME(L),0.0)  
            ENDIF  
          ENDDO  
          DO L=2,LA  
            DIFQVOL=QSUME(L)-QSUMTMP(L)  
            DO K=1,KC  
              QSUM(L,K)=QSUM(L,K)-DIFQVOL*DZC(K)  
            ENDDO  
            QSUME(L)=QSUMTMP(L)  
          ENDDO  
        ENDIF  
  C  
  C **  SET OLD TIME LEVEL TERMS IN CONTINUITY EQUATION FOR  
  C **  NON BOUNDARY POINTS  
  C **  HRU=HMU*DYU/DXU & HRV=HMV*DXV/DYV  
  C **  DXYIP=1/(DXP*DYP)  
  C  
        DO L=2,LA  
          LN=LNC(L)  
          FP(L)=DELTI*DXYP(L)*P1(L)-0.5*G*(UHDY1E(LEAST(L))-UHDY1E(L)  
       &                                  +VHDX1E(LN )-VHDX1E(L))  
        ENDDO  
  C  
  C **  ADVANCE EXTERNAL VARIABLES FOR THREE TIME LEVEL STEP  
  C  
        IF(ISTL_.EQ.3)THEN  
          DO L=2,LA  
            UHDY2E(L)=UHDY1E(L)  
            VHDX2E(L)=VHDX1E(L)  
            UHDY1E(L)=UHDYE(L)  
            VHDX1E(L)=VHDXE(L)  
            U1V(L)=UV(L)  
            V1U(L)=VU(L)  
            P1(L)=P(L)  
            H1U(L)=HU(L)  
            H1V(L)=HV(L)  
            H1UI(L)=HUI(L)  
            H1VI(L)=HVI(L)  
            H2P(L)=H1P(L)  
            H1P(L)=HP(L)  
          ENDDO
          IF(ISGWIE.GE.1)THEN
  	      DO L=2,LA  
  	        AGWELV2(L)=AGWELV1(L)  
  	        AGWELV1(L)=AGWELV(L)  
  	      ENDDO  
          ENDIF  
        ENDIF
  C  
  C **  SET NEW TIME LEVEL TERMS IN CONTINUITY EQUATION INCLUDING  
  C **  HOST-GUEST CHANNAL INTERACTION FOR NON BOUNDARY POINTS  
  C  
        C1=0.5*G
        DO L=2,LA  
          LN=LNC(L)  
          ! ***  THE SUB & SVB SWITCHES ALREADY ACCOUNTED FOR
          FP(L)=FP(L)-C1*(FUHDYE(LEAST(L))-FUHDYE(L)  
       &                 +FVHDXE(LN )-FVHDXE(L)  
       &                -2.0*QSUME(L) )  
        ENDDO  
  C
        IF(ISGWIE.GE.1)THEN  
          DO L=2,LA  
            FP(L)=FP(L)-G*SPB(L)*(RIFTR(L)+EVAPSW(L))  
          ENDDO  
        ENDIF  
  C
        IF(ISDSOLV.GE.1.AND.DEBUG)THEN  
          OPEN(1,FILE='FP.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
          WRITE(1,1001)N,ISTL_  
          DO L=2,LA  
            WRITE(1,1001)IL(L),JL(L),FP1(L),FUHDYE(L),FUHDYE(LEAST(L)),  
       &        FVHDXE(L),FVHDXE(LNC(L)),QSUME(L),RIFTR(L),EVAPSW(L)  
          ENDDO  
          CLOSE(1)  
          IF(N.EQ.1)THEN  
            OPEN(1,FILE='FP1.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1001)IL(L),JL(L),FP1(L),FUHDYE(L),FUHDYE(LEAST(L)),  
       &          FVHDXE(L),FVHDXE(LNC(L)),QSUME(L),RIFTR(L),EVAPSW(L)  
            ENDDO  
            CLOSE(1)  
          ENDIF  
          IF(N.EQ.2)THEN  
            OPEN(1,FILE='FP2.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1001)IL(L),JL(L),FP1(L),FUHDYE(L),FUHDYE(LEAST(L)),  
       &          FVHDXE(L),FVHDXE(LNC(L)),QSUME(L),RIFTR(L),EVAPSW(L)  
            ENDDO  
            CLOSE(1)  
          ENDIF  
        ENDIF  
  C
        C1=-0.5*DELTD2*G  
        DO L=2,LA  
          CS(L)=C1*SVB(L  )*HRVO(L  )*RCY(L  )*HVTMP(L  )  
          CW(L)=C1*SUB(L  )*HRUO(L  )*RCX(L  )*HUTMP(L  )  
          CE(L)=C1*SUB(LEAST(L))*HRUO(LEAST(L))*RCX(LEAST(L))*HUTMP(LEAST(L))  
        ENDDO  
        DO L=2,LA  
          LN=LNC(L)  
          CN(L)=C1*SVB(LN )*HRVO(LN )*RCY(LN )*HVTMP(LN )  
        ENDDO  
  C
  C *** APPLY THE OPEN BOUNDARY CONDITIONS
  C
        IF(NBCSOP.GT.0) CALL SETOPENBC(DELT,DELTD2,DELTI,HUTMP,HVTMP)   
  C  
        ! *** SET THE CENTER
        DO L=2,LA  
          CC(L)=DELTI*DXYP(L)-CS(L)-CW(L)-CE(L)-CN(L)  
        ENDDO  
  C
        ! *** SCALE COEFFICIENTS IN EXTERNAL MODEL LINEAR EQUATION SYSTEM  
        CCMNM=1.E+18
        DO L=2,LA  
          CCMNM=MIN(CCMNM,CC(L))  
          FPTMP(L)=FP(L)  
        ENDDO  
        CCMNMI=1./CCMNM  
  C
  C *** APPLY THE OPEN BOUNDARY CONDITIONS FOR ADJACENT CELLS
  C
        IF(NBCSOP.GT.0) CALL SETOPENBC2   
  C
        CC(1)=1.  
        CC(LC)=1.  
  C  
  C **  SCALE BY MINIMUM DIAGONAL  
  C  
        IF(IRVEC.EQ.9)THEN  
          DO L=2,LA  
            CCS(L)=CS(L)*CCMNMI  
            CCW(L)=CW(L)*CCMNMI  
            CCE(L)=CE(L)*CCMNMI  
            CCN(L)=CN(L)*CCMNMI  
            CCC(L)=CC(L)*CCMNMI  
            FPTMP(L)=FPTMP(L)*CCMNMI  
            CCCI(L)=1./CCC(L)  
          ENDDO  
        ENDIF  
  C  
  C **  SCALE TO NORMAL FORM  
  C  
        IF(IRVEC.EQ.99)THEN  
          DO L=2,LA  
            CCS(L)=CS(L)/SQRT( CC(L)*CC(LSC(L)) )  
            CCW(L)=CW(L)/SQRT( CC(L)*CC(LWEST(L)   ) )  
            CCE(L)=CE(L)/SQRT( CC(L)*CC(LEAST(L)   ) )  
            CCN(L)=CN(L)/SQRT( CC(L)*CC(LNC(L)) )  
            CCC(L)=1.  
            FPTMP(L)=FPTMP(L)/SQRT( CC(L) )  
            P(L)=P(L)*SQRT( CC(L) )  
            CCCI(L)=1.  
          ENDDO  
        ENDIF  
  C
        ! *** CALL THE PRECONDITIONED CONJUGATE GRADIENT SOLVER 
        CALL CONGRAD (ISTL_)  
  C
        IF(IRVEC.EQ.99)THEN  
          DO L=2,LA  
            P(L)=P(L)/SQRT( CC(L) )  
          ENDDO  
        ENDIF  
  C  
  C ** DIAGNOSTICS  
  C  
        IF(ISDSOLV.GE.1.AND.DEBUG)THEN  
          OPEN(1,FILE='EQCOEF.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
          WRITE(1,1001)N,ISTL_  
          DO L=2,LA  
            SURFTMP=GI*P(L)  
            WRITE(1,1001)IL(L),JL(L),CS(L),CW(L),CC(L),CE(L),CN(L),  
       &        FP(L),SURFTMP  
          ENDDO  
          CLOSE(1)  
          IF(N.EQ.1)THEN  
            OPEN(1,FILE='EQCOEF1.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              SURFTMP=GI*P(L)  
              WRITE(1,1001)IL(L),JL(L),CS(L),CW(L),CC(L),CE(L),CN(L),  
       &          FP(L),SURFTMP  
            ENDDO  
            CLOSE(1)  
          ENDIF  
          IF(N.EQ.2)THEN  
            OPEN(1,FILE='EQCOEF2.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              SURFTMP=GI*P(L)  
              WRITE(1,1001)IL(L),JL(L),CS(L),CW(L),CC(L),CE(L),CN(L),  
       &          FP(L),SURFTMP  
            ENDDO  
            CLOSE(1)  
          ENDIF  
        ENDIF  
        IF(ISDSOLV.GE.1.AND.DEBUG)THEN  
          OPEN(1,FILE='EQTERM.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
          WRITE(1,1001)N,ISTL_  
          DO L=2,LA  
            WRITE(1,1001)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),  
       &        HRVO(L),HUTMP(L),HVTMP(L)  
          ENDDO  
          CLOSE(1)  
          IF(N.EQ.1)THEN  
            OPEN(1,FILE='EQTERM1.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1001)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),  
       &          HRVO(L),HUTMP(L),HVTMP(L)  
            ENDDO  
            CLOSE(1)  
          ENDIF  
          IF(N.EQ.2)THEN  
            OPEN(1,FILE='EQTERM2.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1001)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),  
       &          HRVO(L),HUTMP(L),HVTMP(L)  
            ENDDO  
            CLOSE(1)  
          ENDIF  
        ENDIF  
   1001 FORMAT(2I5,10(1X,E12.4))  
   1002 FORMAT(3I4,10(1X,E9.2))  
  C  
  C **  CALCULATE UHEX AND VHEX AND TOTAL DEPTHS AT TIME LEVEL (N+1)  
  C **  HRU=SUB*DYU/DXU & HRV=SVB*DXV/DYV  
  C  
        DO L=2,LA  
          LS=LSC(L)  
          UHDYE(L)=SUB(L)*( FUHDYE(L)  
       &      -DELTD2*HRUO(L)*RCX(L)*HUTMP(L)*(P(L)-P(LWEST(L))) )  
          VHDXE(L)=SVB(L)*( FVHDXE(L)  
       &      -DELTD2*HRVO(L)*RCY(L)*HVTMP(L)*(P(L)-P(LS )) ) 
        ENDDO
  C
        DO L=2,LA  
          UHE(L)=UHDYE(L)*DYIU(L)  
          VHE(L)=VHDXE(L)*DXIV(L)  
        ENDDO  
  C  
  C **  CALCULATE REVISED CELL DEPTHS BASED ON NEW HORIZONTAL  
  C **  TRANSPORTS AT (N+1)  
  C  
        IF(ISTL_.EQ.3)THEN  
          DO L=2,LA  
            LN=LNC(L)  
            HP(L)=H2P(L)+DELT*DXYIP(L)*(QSUME(L)  
       &        -0.5*(UHDYE(LEAST(L))+UHDY2E(LEAST(L))-UHDYE(L)-UHDY2E(L)  
       &             +VHDXE(LN )+VHDX2E(LN )-VHDXE(L)-VHDX2E(L)))
          ENDDO  
        ELSE  
          DO L=2,LA  
            LN=LNC(L)  
            HP(L)=H1P(L)+DELT*DXYIP(L)*(QSUME(L)  
       &        -0.5*(UHDYE(LEAST(L))+UHDY1E(LEAST(L))-UHDYE(L)-UHDY1E(L)  
       &             +VHDXE(LN )+VHDX1E(LN )-VHDXE(L)-VHDX1E(L)))  
          ENDDO
        ENDIF  
  C
        IF(ISGWIE.GE.1)THEN  
          DO L=2,LA  
            HP(L)=HP(L)-DELT*DXYIP(L)*(RIFTR(L)+EVAPSW(L))  
          ENDDO  
        ENDIF
  C
        ! *** APPLY OPEN BOUNDARYS 
        DO LL=1,NBCSOP
          L=LOBCS(LL)
          HP(L)=GI*P(L)-BELV(L)  
        ENDDO 
  C  
  C **  PERFORM FINAL UPDATES OF P,HU, AND HV  
  C  
        DO L=2,LA  
          P(L)=G*(HP(L)+BELV(L))  
        ENDDO  
        DO L=2,LA  
          LS=LSC(L)  
          HU(L)=0.5*(DXYP(L)*HP(L)+DXYP(LWEST(L))*HP(LWEST(L)))*DXYIU(L)  
          HV(L)=0.5*(DXYP(L)*HP(L)+DXYP(LS )*HP(LS ))*DXYIV(L)  
        ENDDO  
        DO L=2,LA  
          HPI(L)=1./HP(L)  
          HUI(L)=1./HU(L)  
          HVI(L)=1./HV(L)  
        ENDDO  
  C  
  C **  PERFORM UPDATE ON GROUNDWATER ELEVATION  
  C  
        IF(ISGWIE.GE.1)THEN  
          DO L=2,LA  
            QSUM(L,KC)=QSUM(L,KC)-EVAPSW(L)  
            QSUM(L,1 )=QSUM(L,1 )-RIFTR(L)  
          ENDDO  
  C  
  C       INFILTRATION STEP  
  C  
          RNPORI=1./RNPOR  
          IF(ISTL_.EQ.3)THEN  
            DO L=2,LA  
              AGWELV(L)=AGWELV2(L)+RNPORI*DELT*DXYIP(L)*RIFTR(L)  
            ENDDO  
          ELSE  
            DO L=2,LA  
              AGWELV(L)=AGWELV1(L)+RNPORI*DELT*DXYIP(L)*RIFTR(L)  
            ENDDO  
          ENDIF  
          DO L=2,LA  
            AGWELV(L)=MIN(AGWELV(L),BELV(L))  
          ENDDO  
  C  
  C       ET STEP  
  C  
          DO L=2,LA  
            IF(EVAPCVT.LT.0.)THEN  
              SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/  
       &          (1.+0.00412*TEM(L,KC))))  
             EVAPT(L)=CLEVAP(L)*0.7464E-3*WINDST(L)*(SVPW-VPA(L))/PATMT(L)  
            ENDIF  
            ETGWTMP=EVAPT(L)-EVAPSW(L)*DXYIP(L)  
            ETGWTMP=MAX(ETGWTMP,0.0)  
            ETGWAVL=RNPOR*DELTI*(AGWELV(L)-BELAGW(L))  
            ETGWAVL=MAX(ETGWAVL,0.0)  
            ETGWTMP=MIN(ETGWTMP,ETGWAVL)  
            EVAPGW(L)=ETGWTMP*DXYP(L)  
          ENDDO  
          DO L=2,LA  
            AGWELV(L)=AGWELV(L)-RNPORI*DELT*DXYIP(L)*EVAPGW(L)  
          ENDDO  
          DO L=2,LA  
            AGWELV(L)=MAX(AGWELV(L),BELAGW(L))  
          ENDDO  
        ENDIF  
  C  
  C **  CHECK FOR NEGATIVE DEPTHS  
  C  
        IF(ISNEGH.GE.1)THEN  
          INEGFLG=0  
          DO L=2,LA  
            IF(HP(L).LT.0.)THEN  
              INEGFLG=1  
              LN=LNC(L)  
              WRITE (6,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)  
              WRITE (6,6061)IL(L),JL(L),HU(L),H1U(L)  
              WRITE (6,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))  
              WRITE (6,6063)IL(L),JL(L),HV(L),H1V(L)  
              WRITE (6,6064)IL(L),JL(L),HV(LN),H1V(LN)  
              WRITE (8,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)  
              WRITE (8,6061)IL(L),JL(L),HU(L),H1U(L)  
              WRITE (8,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))  
              WRITE (8,6063)IL(L),JL(L),HV(L),H1V(L)  
              WRITE (8,6064)IL(L),JL(L),HV(LN),H1V(LN)  
            ENDIF  
          ENDDO  
          DO L=2,LA  
            IF(HU(L).LT.0.)THEN  
              INEGFLG=1  
              LN=LNC(L)  
              WRITE (6,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)  
              WRITE (6,6061)IL(L),JL(L),HU(L),H1U(L)  
              WRITE (6,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))  
              WRITE (6,6063)IL(L),JL(L),HV(L),H1V(L)  
              WRITE (6,6064)IL(L),JL(L),HV(LN),H1V(LN)  
              WRITE (8,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)  
              WRITE (8,6061)IL(L),JL(L),HU(L),H1U(L)  
              WRITE (8,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))  
              WRITE (8,6063)IL(L),JL(L),HV(L),H1V(L)  
              WRITE (8,6064)IL(L),JL(L),HV(LN),H1V(LN)  
            ENDIF  
          ENDDO  
          DO L=2,LA  
            IF(HV(L).LT.0.)THEN  
              INEGFLG=1  
              LN=LNC(L)  
              WRITE (6,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)  
              WRITE (6,6061)IL(L),JL(L),HU(L),H1U(L)  
              WRITE (6,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))  
              WRITE (6,6063)IL(L),JL(L),HV(L),H1V(L)  
              WRITE (6,6064)IL(L),JL(L),HV(LN),H1V(LN)  
              WRITE (8,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)  
              WRITE (8,6061)IL(L),JL(L),HU(L),H1U(L)  
              WRITE (8,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))  
              WRITE (8,6063)IL(L),JL(L),HV(L),H1V(L)  
              WRITE (8,6064)IL(L),JL(L),HV(LN),H1V(LN)  
            ENDIF  
          ENDDO  
        ENDIF  
        IF(ISNEGH.EQ.2)THEN  
          IF(INEGFLG.EQ.1)THEN  
            CALL RESTOUT(1)  
            OPEN(1,FILE='EQCOEF.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQCOEF.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              SURFTMP=GI*P(L)  
              WRITE(1,1001)IL(L),JL(L),CCS(L),CCW(L),CCC(L),CCE(L),CCN(L),  
       &          FPTMP(L),SURFTMP  
            ENDDO  
            CLOSE(1)  
            OPEN(1,FILE='EQTERM.OUT',STATUS='UNKNOWN')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='EQTERM.OUT',POSITION='APPEND',STATUS='UNKNOWN')  
            WRITE(1,1001)N,ISTL_  
            DO L=2,LA  
              WRITE(1,1002)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),  
       &          HRVO(L),HUTMP(L),HVTMP(L)  
            ENDDO  
            CLOSE(1)  
            OPEN(1,FILE='CFLMAX.OUT')  
            CLOSE(1,STATUS='DELETE')  
            OPEN(1,FILE='CFLMAX.OUT')  
            DO L=2,LA  
              WRITE(1,1991)IL(L),JL(L),(CFLUUU(L,K),K=1,KC)  
              WRITE(1,1992)(CFLVVV(L,K),K=1,KC)  
              WRITE(1,1992)(CFLWWW(L,K),K=1,KC)  
              WRITE(1,1992)(CFLCAC(L,K),K=1,KC)  
            ENDDO  
            CLOSE(1)  
            STOP  
          ENDIF  
        ENDIF  
   1991 FORMAT(2I5,12F7.2)  
   1992 FORMAT(10X,12F7.2)  
   6060 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HP,H1P,H2P =',3(2X,E12.4))  
   6061 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HUW,H1UW =',2(2X,E12.4))  
   6062 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HUE,H1UE =',2(2X,E12.4))  
   6063 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HVS,H1VS =',2(2X,E12.4))  
   6064 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HVN,H1VN =',2(2X,E12.4))  
   6910 FORMAT('  DRYING AT N,I,J =',I10,2I6,'  HP,H1P,H2P ='  
       &    ,3(2X,E12.4))  
   6911 FORMAT('  DRY W FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))  
   6912 FORMAT('  DRY E FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))  
   6913 FORMAT('  DRY S FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))  
   6914 FORMAT('  DRY N FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))  
   6920 FORMAT('  WETTING AT N,I,J =',I10,2I6,' HP,H1P,H2P ='  
       &    ,3(2X,E12.4))  
   6921 FORMAT('  WET S FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))  
   6922 FORMAT('  WET W FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))  
   6923 FORMAT('  WET E FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))  
   6924 FORMAT('  WET N FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))  
   6930 FORMAT('  WET BY VOL  N,I,J =',I10,2I6,' HP,H1P,H2P ='  
       &    ,3(2X,E12.4))  
   6940 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HP,H1P,H2P ='  
       &    ,3(2X,E12.4))  
   6941 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HUE,HP,H1P ='  
       &    ,3(2X,E12.4))  
   6942 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HUW,HP,H1P ='  
       &    ,3(2X,E12.4))  
   6943 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HVS,HP,H1P ='  
       &    ,3(2X,E12.4))  
   6944 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HVN,HP,H1P ='  
       &    ,3(2X,E12.4))  
   6945 FORMAT('  RESOLVE NEG,  N,I,J =',I10,2I6,' HP,H1P,H2P ='  
       &    ,3(2X,E12.4))  
   6950 FORMAT('  RESOLVE, NEG DEP N,I,J =',I10,2I6,' HP,H1P,H2P ='  
       &    ,3(2X,E12.4))  
  C     
  C **  CALCULATE THE EXTERNAL DIVERGENCE  
  C  
        IF(ISDIVEX.EQ.1)THEN  
          DIVEXMX=0.  
          DIVEXMN=1000000.  
          IF(ISTL_.EQ.3)THEN  
            DO L=2,LA  
              IF(SPB(L).NE.0)THEN  
                LN=LNC(L)  
                DIVEX=SPB(L)*(DXYP(L)*(HP(L)-H2P(L))*DELTI  
       &            +0.5*(UHDYE(LEAST(L))+UHDY2E(LEAST(L))-UHDYE(L)-UHDY2E(L)  
       &            +VHDXE(LN)+VHDX2E(LN)-VHDXE(L)-VHDX2E(L))-QSUME(L)  
       &            +RIFTR(L)+EVAPSW(L))  
                IF(DIVEX.GT.DIVEXMX)THEN  
                  DIVEXMX=DIVEX  
                  LMAX=L  
                ENDIF  
                IF(DIVEX.LT.DIVEXMN)THEN  
                  DIVEXMN=DIVEX  
                  LMIN=L  
                ENDIF  
              ENDIF  
            ENDDO  
          ELSE  
            DO L=2,LA  
              IF(SPB(L).NE.0)THEN  
                LN=LNC(L)  
                DIVEX=SPB(L)*(DXYP(L)*(HP(L)-H1P(L))*DELTI  
       &            +0.5*(UHDYE(LEAST(L))+UHDY1E(LEAST(L))-UHDYE(L)-UHDY1E(L)  
       &            +VHDXE(LN)+VHDX1E(LN)-VHDXE(L)-VHDX1E(L))-QSUME(L)  
       &            +RIFTR(L)+EVAPSW(L))  
                IF(DIVEX.GT.DIVEXMX)THEN  
                  DIVEXMX=DIVEX  
                  LMAX=L  
                ENDIF  
                IF(DIVEX.LT.DIVEXMN)THEN  
                  DIVEXMN=DIVEX  
                  LMIN=L  
                ENDIF  
              ENDIF  
            ENDDO  
          ENDIF  
          IMAX=IL(LMAX)  
          JMAX=JL(LMAX)  
          IMIN=IL(LMIN)  
          JMIN=JL(LMIN)  
          WRITE(6,6628)DIVEXMX,IMAX,JMAX  
          WRITE(6,6629)DIVEXMN,IMIN,JMIN  
        ENDIF  
    566 FORMAT('  I=',I5,3X,'J=',I5,3X,'HP=',F12.4)  
   6628 FORMAT('  DIVEXMX=',E13.5,5X,2I10)  
   6629 FORMAT('  DIVEXMN=',E13.5,5X,2I10)  
  C  
  C **  UPDATE ZERO DIMENSION VOLUME BALANCE  
  C  
        IF(ISDRY.GE.1.AND.ISTL_.EQ.3)THEN  
          VOLADD=0.  
          DO L=2,LA  
            IF(SPB(L).NE.0)THEN  
              VOLADD=VOLADD+QSUME(L)-RIFTR(L)-EVAPSW(L)  
            ENDIF  
          ENDDO  
          VOLADD=VOLADD*DT  
          VOLZERD=VOLZERD+VOLADD  
          VETZERD=VETZERD+VOLADD+DT*EVAPSW(L)  
        ENDIF  
   5303 FORMAT(2X,F10.4,2X,F10.5,3(2X,E13.5))  
        RETURN  
        END  
  
  