Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:40:14 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
613fa3996d Andr*0002       SUBROUTINE TILE (
                0003      I                 NCH, DTSTEP, ITYP, TRAINL,TRAINC, TSNOW,  UM,
                0004      I                 ETURB,  DEDQA,  DEDTC,  HSTURB, DHSDQA, DHSDTC,
                0005      I                 TM,     QM,     CD,     SUNANG, PARDIR,  PARDIF,
                0006      I                 SWNET,  HLWDWN, PSUR,    ZLAI,   GREEN,  Z2,
                0007      I                 SQSCAT, RSOIL1, RSOIL2,   RDC,    U2FAC,
                0008      I                 QSATTC, DQSDTC, ALWRAD, BLWRAD,
                0009      U                 TC, TD, QA, SWET1, SWET2, SWET3, CAPAC, SNOW,
                0010      O                 EVAP, SHFLUX, RUNOFF, BOMB,
                0011      O                 EINT,   ESOI,   EVEG,   ESNO, SMELT, HLATN,
                0012      O                 HLWUP,GDRAIN,RUNSRF,FWSOIL,
                0013      O                 STRDG1, STRDG2, STRDG3, STRDG4,
                0014      O                 STRDG5, STRDG6, STRDG7, STRDG8, STRDG9
                0015      &                )
                0016 
aab9c05604 Jean*0017 C   SCCS VERSION @(#)lsm.f      1.2 11/3/92
613fa3996d Andr*0018 
                0019 C****
                0020 C****       This subroutine computes the outgoing fluxes of sensible and
                0021 C**** latent heat from the land surface and updates the surface prognostic
                0022 C**** variables.
                0023 C****
                0024       IMPLICIT NONE
                0025 C
                0026 C****
4dad0f083a Andr*0027 #include "sibber.h"
613fa3996d Andr*0028 C****
                0029       INTEGER NCH
                0030       INTEGER ITYP(NCH)
a456aa407c Andr*0031       _RL  DTSTEP,  TRAINL(NCH), TRAINC(NCH), TSNOW(NCH), UM(NCH),
613fa3996d Andr*0032      &       ETURB(NCH),  DEDQA(NCH), HSTURB(NCH), DHSDTC(NCH),
                0033      &         TM (NCH),     CD(NCH), SUNANG(NCH), DHSDQA(NCH),
                0034      &         QM (NCH), PARDIR(NCH), PARDIF(NCH),  SWNET(NCH),
                0035      &      HLWDWN(NCH),   PSUR(NCH),   ZLAI(NCH),  GREEN(NCH),
                0036      &          Z2(NCH), SQSCAT(NCH),  DEDTC(NCH)
a456aa407c Andr*0037       _RL  RSOIL1(NCH), RSOIL2(NCH),    RDC(NCH),  U2FAC(NCH),
613fa3996d Andr*0038      &     QSATTC(NCH), DQSDTC(NCH),  ALWRAD(NCH), BLWRAD(NCH),
                0039      &          TC(NCH),     TD(NCH),     QA(NCH),   BOMB(NCH),
aab9c05604 Jean*0040      &       SWET1(NCH),  SWET2(NCH),  SWET3(NCH),  CAPAC(NCH),
613fa3996d Andr*0041      &        SNOW(NCH),   EVAP(NCH), SHFLUX(NCH), RUNOFF(NCH)
a456aa407c Andr*0042       _RL   EINT(NCH),    ESOI(NCH),   EVEG(NCH),   ESNO(NCH),
613fa3996d Andr*0043      &     STRDG1(NCH),  STRDG2(NCH), STRDG3(NCH), STRDG4(NCH),
                0044      &     STRDG5(NCH),  STRDG6(NCH), STRDG7(NCH), STRDG8(NCH),
                0045      &     STRDG9(NCH),
                0046      &      SMELT(NCH),   HLATN(NCH),  HLWUP(NCH), GDRAIN(NCH),
                0047      &     RUNSRF(NCH),  FWSOIL(NCH)
                0048 C****
                0049       INTEGER ChNo
a456aa407c Andr*0050       _RL     SWET(nch,NLAY),   VGPSAT(NTYPS),   VGCSAT (NTYPS),
613fa3996d Andr*0051      &       VGZDEP(NLAY,NTYPS),   VGSLOP(NTYPS),            DELTC,
                0052      &                    DELEA,    VGPH1(NTYPS),     VGPH2(NTYPS),
                0053      &            VGRPLN(NTYPS),    CSOIL0(NTYPS),          WSOI12,
                0054      &             VGBEE(NTYPS),    DELZ12(NTYPS),
4dad0f083a Andr*0055      &            DELZ23(NTYPS)
613fa3996d Andr*0056 C****
                0057 
a456aa407c Andr*0058       _RL PHLAY(nch,NLAY), AKLAY(nch,NLAY), SWET12(nch),
613fa3996d Andr*0059      &      CSOIL(nch),
                0060      &       RCUN(nch),     VPDSTR(nch),   ESATTX(nch),
                0061      &       VPDSTX(nch),   VGBEEX(nch)
a456aa407c Andr*0062       _RL EMAXRT(nch), VGWMAX(NLAY,NTYPS), FTEMP(nch),
613fa3996d Andr*0063      &        PHR(nch),      SOILCO(nch),    RC(nch),
                0064      &        EAX(nch),          TX(nch),   RCX(nch),
                0065      &     DRCDTC(nch),      SATCAP(nch),   PAR(nch),
                0066      &       PDIR(nch),       DUMMY(nch)
a456aa407c Andr*0067       _RL FTEMPX(nch),  DRCDEA(nch), VGPSAX(nch), VGCSAX(nch),
613fa3996d Andr*0068      &                VGZDEX(NLAY,nch), VGSLOX(nch), VGPH1X(nch),
                0069      &     VGPH2X(nch),  VGRPLX(nch)
a456aa407c Andr*0070       _RL  DEDEA(nch),  DHSDEA(nch),     EM(nch), ESATTC(nch),
613fa3996d Andr*0071      &     DESDTC(nch),      EA(nch),     RA(nch),   ALHX(nch),
                0072      &     WETEQ1(nch), WETEQ2(nch),
                0073      &        RX1(nch),     RX2(nch), SNWFRC(nch), POTFRC(nch),
                0074      &     ESNFRC(nch),  EIRFRC(nch),   FCAN(nch),  EPFRC,
                0075      &     DEFCIT, EADJST, RTBS
a456aa407c Andr*0076       _RL cmpbug
613fa3996d Andr*0077 
                0078 C****
                0079       DATA VGWMAX   /8.4,    621.6,      840.0,
                0080      2               8.4,    621.6,      840.0,
                0081      3               8.4,    621.6,      840.0,
                0082      4               8.4,    197.4,      420.0,
                0083      5               8.704,  204.5,      435.2,
                0084      6               8.4,     71.4,      420.0,
                0085      7               4.000,   4.000,    130.56,
                0086      8               4.000,   4.000,    130.56,
                0087      9               8.704,   73.984,    130.56,
                0088      1               4.000,   4.000,    130.56
                0089      &              /
aab9c05604 Jean*0090       DATA VGBEE /    4.,    4.,    4.,    4.,  1.69,    4.,
613fa3996d Andr*0091      &                4.,  1.69,  1.69, 1.69/
aab9c05604 Jean*0092 C      DATA VGBEE /    7.,    7.,    7.,    7.,  1.69,    7.,
613fa3996d Andr*0093 C     &                7.,  1.69,  1.69, 1.69/
aab9c05604 Jean*0094 c      DATA VGBEE /    4.,    4.,    4.,    4.,  0.95,    4.,
613fa3996d Andr*0095 c     &                4.,  0.95,  0.95, 0.95/
                0096 C      DATA VGBEE /7, 7, 7, 7, 2, 7, 7, 4, 2, 4/
                0097       DATA VGPSAT/    -0.281,     -0.281,       -0.281,       -0.281,
                0098      5                -0.073,     -0.281,       -0.281,       -0.073,
                0099      9                -0.073 ,    -0.073/
                0100 C      DATA VGPSAT/    -0.086,     -0.086,       -0.086,       -0.086,
                0101 C     5                -0.073,     -0.086,       -0.086,       -0.073,
                0102 C     9                -0.073 ,    -0.073/
                0103 c      DATA VGPSAT/    -0.281,     -0.281,       -0.281,       -0.281,
                0104 c     5                -0.014,     -0.281,       -0.281,       -0.014,
                0105 c     9                -0.014 ,    -0.014/
                0106 c      DATA VGCSAT /  0.00002,  0.00002,    0.00002,    0.00002,
                0107 c     5               0.0000583,  0.00002,    0.00002,    0.0000583,
                0108 c     9               0.0000583,  0.0000583
                0109 c     &              /
                0110       DATA VGCSAT /  0.0000012,  0.0000012,    0.0000012,    0.0000012,
                0111      5               0.0000583,  0.0000012,    0.0000012,    0.0000583,
                0112      9               0.0000583,  0.0000583
                0113      &              /
                0114 C**** - - - - - - - -
                0115 C**** THE FOLLOWING 3 VARIABLES MUST MAINTAIN CONSISTENT VALUES.  ZDEP12(N) =
                0116 C****  (VGZDEP(1,N)+VGZDEP(2,N))/2. ; SIMILARLY FOR ZDEP23(N).
                0117 C****
                0118       DATA VGZDEP /.02, 1.48, 2.00,
                0119      2             .02, 1.48, 2.00,
                0120      3             .02, 1.48, 2.00,
                0121      4             .02,  .47, 1.00,
                0122      5             .02,  .47, 1.00,
                0123      6             .02,  .17, 1.00,
                0124      7             .0092,  .0092,  .30,
                0125      8             .0092,  .0092,  .30,
                0126      9             .02,  .17,  .30,
                0127      1             .0092,  .0092,  .30
                0128      &            /
                0129       DATA DELZ12 / 0.75, 0.75, 0.75, 0.245, 0.245, 0.095, 0.0092,
                0130      8        0.0092, 0.095, 0.0092 /
                0131       DATA DELZ23 / 1.74, 1.74, 1.74, 0.735, 0.735, 0.585, 0.155,
                0132      8        0.155, 0.235, 0.155 /
                0133 C**** - - - - - - - -
                0134       DATA VGSLOP /    .1736,     .1736,     .1736,     .1736,
                0135      5                 1.000,     .1736,     .1736,     1.000,
                0136      &                 1.000 , 1.000 /
                0137 c      DATA VGSLOP /    .1736,     .1736,     .1736,     .1736,
                0138 c     5                 .0872,     .1736,     .1736,     .0872,
                0139 c     &                 .0872 , .0872 /
                0140 c      DATA VGSLOP / 1., 1., 1.,  1., 1., 1.,  1., 1., 1.,1./
                0141       DATA VGPH1 /     -100.,     -190.,     -200.,     -120.,
                0142      5                 -200.,     -200.,     -200.,      -10.,
                0143      9                  -10.,      -10. /
                0144       DATA VGPH2 /     -500.,     -250.,     -250.,     -230.,
                0145      5                 -400.,     -400.,     -250.,     -100.,
                0146      9                 -100.,     -100. /
                0147       DATA VGRPLN /0.245E+09, 0.245E+09, 0.245E+09, 0.250E+09,
                0148      5             0.250E+09, 0.250E+09, 0.245E+09,   0.1E+09,
                0149      9               0.1E+09, 0.100E+09 /
aab9c05604 Jean*0150 C****
613fa3996d Andr*0151       DATA DELTC /0.01/, DELEA /0.001/
                0152 
aab9c05604 Jean*0153 c Note: SNWMID & CSOIL0 parameters modified from (Koster/Suarez)
613fa3996d Andr*0154 c       to obtain improved albedo, radswg, and annual/diurnal cycle (L.Takacs 2/17/00)
                0155 c ------------------------------------------------------------------------------------
                0156       DATA CSOIL0 /175000.,175000.,175000.,175000.,175000.,
                0157      .             175000.,175000.,120000.,175000., 70000./
4dad0f083a Andr*0158 #include "snwmid.h"
613fa3996d Andr*0159 
                0160 C****
aab9c05604 Jean*0161 C**** ---------------------------------------------------------------------
613fa3996d Andr*0162 C****
                0163 CFPP$ EXPAND (TMPFAC, RCANOP, SOIL, WUPDAT, SMFAC)
                0164 C****
                0165 C**** Expand data as specified by ITYP into arrays of size NCH.
                0166 C****
                0167       DO 50 ChNo = 1, NCH
                0168       BOMB(CHNO) = 0.
                0169       VGBEEX(ChNo) = VGBEE(ITYP(ChNo))
                0170       VGPSAX(ChNo) = VGPSAT(ITYP(ChNo))
                0171       VGCSAX(ChNo) = VGCSAT(ITYP(ChNo))
                0172       VGZDEX(SFCLY ,ChNo) = VGZDEP(SFCLY ,ITYP(ChNo))
                0173       VGZDEX(ROOTLY,ChNo) = VGZDEP(ROOTLY,ITYP(ChNo))
                0174       VGZDEX(RECHLY,ChNo) = VGZDEP(RECHLY,ITYP(ChNo))
                0175       VGSLOX(ChNo) = VGSLOP(ITYP(ChNo))
                0176       VGPH1X(ChNo) = VGPH1(ITYP(ChNo))
                0177       VGPH2X(ChNo) = VGPH2(ITYP(ChNo))
                0178       VGRPLX(ChNo) = VGRPLN(ITYP(ChNo))
                0179   50  CONTINUE
                0180 
                0181 C****
                0182 C**** Pre-process input arrays as necessary:
                0183 
                0184       DO 100 ChNo = 1, NCH
                0185 
32362b8fa8 Cons*0186       DEDQA(CHNO)  = MAX(  DEDQA(CHNO), 500. _d 0/ALHE )
                0187       DEDTC(CHNO)  = MAX(  DEDTC(CHNO),   0. _d 0)
                0188       DHSDQA(CHNO) = MAX( DHSDQA(CHNO),   0. _d 0)
                0189       DHSDTC(CHNO) = MAX( DHSDTC(CHNO), -10. _d 0)
613fa3996d Andr*0190 
                0191       EM(CHNO)     = QM(CHNO) * PSUR(CHNO) / EPSILON
                0192       EA(CHNO)     = QA(CHNO) * PSUR(CHNO) / EPSILON
                0193       ESATTC(CHNO) = QSATTC(CHNO) * PSUR(CHNO) / EPSILON
                0194       DESDTC(CHNO) = DQSDTC(CHNO) * PSUR(CHNO) / EPSILON
                0195 C      DEDEA(CHNO)  = DEDQA(CHNO) * EPSILON / ( PSUR(CHNO) * ALHX(CHNO) )
                0196       DEDEA(CHNO)  = DEDQA(CHNO) * EPSILON / PSUR(CHNO)
                0197       DHSDEA(CHNO) = DHSDQA(CHNO) * EPSILON / PSUR(CHNO)
                0198 C      DEDTC(CHNO)  = DEDTC(CHNO) / ALHX(CHNO)
                0199 C      ETURB(CHNO)  = ETURB(CHNO) / ALHX(CHNO)
                0200       PAR(CHNO)    = (PARDIR(CHNO) + PARDIF(CHNO) + 1.E-20)
                0201       PDIR(CHNO)   = PARDIR(CHNO) / PAR(CHNO)
                0202       RA(CHNO)     = ONE / ( CD(CHNO) * UM(CHNO) )
                0203       SATCAP(ChNo) = 0.1 * ZLAI(ChNo)
                0204       CSOIL(CHNO)  = CSOIL0(ITYP(ChNo))
32362b8fa8 Cons*0205       SWET(ChNo,SFCLY ) = max( min(SWET1(ChNo),1. _d 0), 0. _d 0)
                0206       SWET(ChNo,ROOTLY) = max( min(SWET2(ChNo),1. _d 0), 0. _d 0)
                0207       SWET(ChNo,RECHLY) = max( min(SWET3(ChNo),1. _d 0), 0. _d 0)
                0208       CAPAC(CHNO)       = max( min(CAPAC(ChNo),SATCAP(CHNO)), 0. _d 0)
613fa3996d Andr*0209 C****
                0210 
                0211       SNWFRC(CHNO) = SNOW(CHNO) / ( SNOW(CHNO) + SNWMID(ITYP(CHNO)) )
32362b8fa8 Cons*0212       FCAN(CHNO) = MIN( 1. _d 0, MAX(0. _d 0,CAPAC(ChNo)/SATCAP(ChNo)) )
613fa3996d Andr*0213       POTFRC(CHNO)=1.-(1.-SNWFRC(CHNO))*(1.-FCAN(CHNO))
                0214 
                0215 
                0216       WSOI12=SWET(ChNo,SFCLY ) * VGWMAX(SFCLY ,ITYP(ChNo)) +
                0217      &       SWET(ChNo,ROOTLY) * VGWMAX(ROOTLY,ITYP(ChNo))
                0218       SWET12(ChNo) = WSOI12 /
                0219      &          (VGWMAX(SFCLY,ITYP(ChNo)) + VGWMAX(ROOTLY,ITYP(ChNo)))
                0220 
                0221       EMAXRT(CHNO) = ( SNOW(CHNO) + CAPAC(CHNO) + WSOI12 ) / DTSTEP
                0222 
                0223 C****
                0224       RUNOFF(CHNO) = 0.
                0225       RUNSRF(CHNO) = 0.
                0226 C**** (SMELT is zeroed in FLUXES)
                0227 C      SMELT(CHNO)  = 0.
                0228       FWSOIL(CHNO) = 0.
                0229  100  CONTINUE
                0230 
                0231 C****
                0232 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                0233 C****       STEP 1: COMPUTE EFFECTIVE RESISTANCE RC FOR ENERGY BALANCE.
                0234 C****          (RCUNST  computes the unstressed resistance,
                0235 C****           VPDFAC  computes the vapor pressure deficit stress term,
                0236 C****           TMPFAC  computes the temperature stress term,
                0237 C****           SOIL    computes soil moisture potentials and conds.,
                0238 C****           SMFAC   computes rc given a leaf water potential stress,
                0239 C****           RSURFP  computes rc given a parallel resist. from the
                0240 C****                   surface,
                0241 C****           RCANOP  computes rc corrected for snow and interception.)
                0242 C****
                0243 
                0244       CALL RCUNST (
                0245      I             NCH, ITYP, SUNANG, SQSCAT, PDIR,
                0246      I             PAR, ZLAI, GREEN,
                0247      O             RCUN
                0248      &            )
                0249 
                0250       CALL VPDFAC (
aab9c05604 Jean*0251      I             NCH,  ITYP,  ESATTC, EA,
613fa3996d Andr*0252      O             VPDSTR
                0253      &            )
                0254 
                0255       CALL TMPFAC (
                0256      I             NCH,  ITYP, TC,
                0257      O             FTEMP
                0258      &            )
                0259 
                0260       CALL SOIL (
                0261      I           NCH, ITYP, SWET12, VGBEEX, VGPSAX, VGCSAX, DELZ12,
                0262      O           PHR, SOILCO, DUMMY
                0263      &          )
                0264 
                0265       cmpbug=0.
                0266 
                0267       CALL SMFAC (
                0268      I            NCH,    ITYP,   ESATTC, EA,   PHR, SOILCO, RCUN,
                0269      I            VPDSTR, FTEMP,  TC,     PSUR, Z2,  RSOIL1, RSOIL2,
                0270      I            VGPH1X, VGPH2X, VGRPLX,
                0271      O            RC
                0272      &           )
                0273 
                0274       CALL RSURFP (
                0275      I             NCH, UM, U2FAC, Z2, RDC, SWET(1,SFCLY),
                0276      I             ESATTC, EA,
                0277      U             RC,
                0278      O             RX1, RX2
                0279      &            )
                0280 
                0281       CALL RCANOP (
                0282      I             NCH, CAPAC, SNOW, SATCAP, RA, ETURB, SNWFRC,
                0283      I             POTFRC,
                0284      U             RC
                0285      &            )
                0286 
                0287 C****
                0288 C**** -    -    -    -    -    -    -    -    -    -    -    -    -    -
                0289 C****
                0290 C**** Compute DRC/DT and DRC/DEA using temperature, v.p. perturbations:
                0291 C****
                0292 
                0293       DO 120 ChNo = 1, NCH
                0294       TX(ChNo) = TC(ChNo) + DELTC
                0295       ESATTX(ChNo) = ESATTC(ChNo) + DESDTC(CHNO) * DELTC
                0296       EAX(ChNo) = EA(ChNo) + DELEA
                0297   120 CONTINUE
                0298 C****
                0299 C**** temperature:
                0300       CALL VPDFAC (NCH, ITYP, ESATTX, EA, VPDSTX)
                0301       CALL TMPFAC (NCH, ITYP, TX, FTEMPX)
                0302       CALL SMFAC (
                0303      I            NCH,    ITYP,   ESATTX, EA,   PHR, SOILCO, RCUN,
                0304      I            VPDSTX, FTEMPX, TX,     PSUR, Z2,  RSOIL1, RSOIL2,
                0305      I            VGPH1X, VGPH2X, VGRPLX,
                0306      O            RCX
                0307      &           )
                0308       CALL RSURFP (NCH, UM, U2FAC, Z2, RDC, SWET(1,SFCLY),
                0309      I             ESATTX, EA,
                0310      &             RCX,
                0311      O             DUMMY,DUMMY
                0312      &            )
                0313       CALL RCANOP (NCH, CAPAC, SNOW, SATCAP, RA, ETURB, SNWFRC,
                0314      &               POTFRC, RCX)
                0315 C****
                0316       DO 125 ChNo = 1, NCH
                0317       DRCDTC(ChNo) = (RCX(ChNo) - RC(ChNo)) / DELTC
                0318   125 CONTINUE
                0319 C****
                0320 C**** vapor pressure:
                0321       CALL VPDFAC (NCH, ITYP, ESATTC, EAX, VPDSTX)
                0322       CALL SMFAC (
                0323      I            NCH,    ITYP,  ESATTC, EAX,  PHR, SOILCO, RCUN,
                0324      I            VPDSTX, FTEMP, TC,     PSUR, Z2,  RSOIL1, RSOIL2,
                0325      I            VGPH1X, VGPH2X, VGRPLX,
                0326      O            RCX
                0327      &           )
                0328       CALL RSURFP (NCH, UM, U2FAC, Z2, RDC, SWET(1,SFCLY),
                0329      I             ESATTC, EAX,
                0330      &             RCX,
                0331      O             DUMMY,DUMMY
                0332      &            )
                0333       CALL RCANOP (NCH, CAPAC, SNOW, SATCAP, RA, ETURB, SNWFRC,
                0334      &               POTFRC, RCX)
                0335 C****
                0336       DO 150 ChNo = 1, NCH
                0337       DRCDEA(ChNo) = (RCX(ChNo) - RC(ChNo)) / DELEA
                0338   150 CONTINUE
                0339 C****
                0340 
                0341 C****
                0342 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                0343 C****       STEP 2: Solve the energy balance at the surface.
                0344 C****
                0345 C****
                0346 C**** Determine effective latent heat of vaporization based on
                0347 C**** assumed fractional coverage of snow.  This requires the
                0348 C**** determination of the fractions of moisture taken from the various
                0349 C**** reservoirs:
                0350 C****         ESNFRC=fraction of total evap. coming from snow
                0351 C****         EIRFRC=fraction of total evap. coming from interception res.
                0352 
                0353       DO 200 CHNO=1,NCH
                0354       RTBS=RX1(CHNO)*RX2(CHNO)/(RX1(CHNO)+RX2(CHNO))
                0355       EPFRC=POTFRC(CHNO) * ( RA(CHNO) + RTBS ) /
                0356      &                    ( RA(CHNO) + POTFRC(CHNO)*RTBS )
                0357       ESNFRC(CHNO)=EPFRC*SNWFRC(CHNO)/
                0358      &                      (SNWFRC(CHNO)+FCAN(CHNO)+1.E-20)
                0359       EIRFRC(CHNO)=EPFRC*FCAN(CHNO)/(SNWFRC(CHNO)+FCAN(CHNO)+1.E-20)
                0360       ALHX(CHNO) = (1.-ESNFRC(CHNO))*ALHE + ESNFRC(CHNO)*ALHS
                0361   200 CONTINUE
                0362 
                0363       CALL FLUXES (
aab9c05604 Jean*0364      I                NCH,   ITYP, DTSTEP, ESATTC, DESDTC,   ALHX,
                0365      I              ETURB,  DEDEA,  DEDTC, HSTURB, DHSDEA, DHSDTC,
                0366      I                 RC, DRCDEA, DRCDTC,
613fa3996d Andr*0367      I              SWNET, HLWDWN, ALWRAD, BLWRAD, ESNFRC,
4dad0f083a Andr*0368      I                 TM,     EM,  CSOIL,   PSUR, EMAXRT, VGWMAX,
613fa3996d Andr*0369      U                 TC,     TD,     EA,   SWET,   SNOW,
                0370      O             RUNOFF,   EVAP, SHFLUX,  SMELT,  HLWUP, BOMB,STRDG1,
                0371      O             STRDG2, STRDG3, STRDG4, STRDG5, STRDG6, STRDG7,
                0372      O             STRDG8, STRDG9
                0373      &            )
                0374 c****
                0375 C****
                0376 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                0377 C****       STEP 3: Update the water quantities in the soil and in the
                0378 C****          interception reservoir.
                0379 C****          (WUPDAT  reduces moisture contents using calculated
                0380 C****           evap.,
                0381 C****           GWATER  computes darcian flux of water between soil
                0382 C****           layers.)
                0383 C****
                0384 
                0385       CALL WUPDAT (
                0386      I               NCH,   ITYP, DTSTEP,
                0387      I               EVAP, SATCAP, VGWMAX, TC, RA, RC,
                0388      I               RX1, RX2,ESNFRC,EIRFRC,
                0389      U               CAPAC, SNOW, SWET, RUNOFF,
                0390      O               EINT, ESOI, EVEG, ESNO, RUNSRF,FWSOIL
                0391      &              )
                0392 
aab9c05604 Jean*0393 C****
613fa3996d Andr*0394 C**** Correct any energy balance inconsistency.
                0395 C****
                0396       DO 300 CHNO=1,NCH
                0397       IF(EVAP(CHNO) .GT. 0.) THEN
                0398         EADJST=(ESNO(CHNO)/DTSTEP) - ESNFRC(CHNO)*EVAP(CHNO)
                0399         SHFLUX(CHNO)=SHFLUX(CHNO)-EADJST*(ALHS-ALHE)
                0400         ENDIF
                0401   300 CONTINUE
                0402 
                0403 
                0404 
                0405       CALL SOIL (
                0406      I             NCH,          ITYP,        SWET (1, SFCLY), VGBEEX,
                0407      I             VGPSAX,       VGCSAX,      DELZ12,
                0408      O             PHLAY(1,SFCLY), AKLAY(1,SFCLY),  DUMMY
                0409      &            )
                0410       CALL SOIL (
                0411      I             NCH,          ITYP,        SWET(1,ROOTLY), VGBEEX,
                0412      I             VGPSAX,       VGCSAX,      DELZ12,
                0413      O             PHLAY(1,ROOTLY), AKLAY(1,ROOTLY),   WETEQ1
                0414      &            )
                0415       CALL SOIL (
                0416      I             NCH,          ITYP,        SWET(1,RECHLY), VGBEEX,
                0417      I             VGPSAX,       VGCSAX,      DELZ23,
                0418      O             PHLAY(1,RECHLY), AKLAY(1,RECHLY),   WETEQ2
                0419      &            )
                0420 
                0421       CALL GWATER (
                0422      I             NCH, ITYP, VGWMAX, PHLAY, AKLAY, TC, DTSTEP,
                0423      I             VGZDEX, VGSLOX, WETEQ1, WETEQ2,
                0424      I             VGPSAX, VGCSAX,
                0425      U             SWET, RUNOFF, GDRAIN
                0426      &            )
                0427 
                0428 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                0429 C****
                0430 C****       STEP 4: Allow precipitation to fill interception reservoir
                0431 C****          and top soil layer
                0432 C****
                0433 
                0434       CALL SOIL (
                0435      I             NCH,          ITYP,        SWET(1,ROOTLY), VGBEEX,
                0436      I             VGPSAX,       VGCSAX,      DELZ12,
                0437      O             PHLAY(1,ROOTLY), AKLAY(1,ROOTLY),   WETEQ1
                0438      &            )
                0439 
                0440       CALL INTERC (
                0441      I             NCH, ITYP, DTSTEP, TRAINL, TRAINC, TSNOW, SATCAP,
                0442      I             VGWMAX, CSOIL, WETEQ1,
                0443      U             TC, CAPAC, SNOW, SWET, RUNOFF, RUNSRF,
                0444      &             SMELT, FWSOIL
                0445      &            )
                0446 
                0447 
                0448 C****
                0449 C****
                0450 C**** Process data for return to GCM:
                0451 
                0452       DO 2000 ChNo = 1, NCH
                0453       QA(CHNO) = EA(CHNO) * EPSILON / PSUR(CHNO)
                0454 C      DEDTC(CHNO) = DEDTC(CHNO) * ALHX(CHNO)
                0455 C      ETURB(CHNO) = ETURB(CHNO) * ALHX(CHNO)
                0456       HLATN (CHNO) = EVAP (CHNO) * ALHX(CHNO)
                0457       SWET1(ChNo) = SWET(ChNo,SFCLY)
                0458       SWET2(ChNo) = SWET(ChNo,ROOTLY)
                0459       SWET3(ChNo) = SWET(ChNo,RECHLY)
                0460 
                0461       IF(SWET1(ChNo).LT.1.E-10) SWET1(CHNO) = 0.0
                0462       IF(SWET2(ChNo).LT.1.E-10) SWET2(CHNO) = 0.0
                0463       IF(SWET3(ChNo).LT.1.E-10) SWET3(CHNO) = 0.0
                0464       IF(CAPAC(ChNo).LT.1.E-10) CAPAC(CHNO) = 0.0
                0465       IF(SNOW (ChNo).LT.1.E-10) SNOW (CHNO) = 0.0
                0466       IF(RUNOFF(ChNo).LT.1.E-10) RUNOFF(CHNO) = 0.0
                0467 
                0468       EINT(CHNO) = EINT(CHNO) * ALHE / DTSTEP
                0469       ESOI(CHNO) = ESOI(CHNO) * ALHE / DTSTEP
                0470       EVEG(CHNO) = EVEG(CHNO) * ALHE / DTSTEP
                0471       ESNO(CHNO) = ESNO(CHNO) * ALHS / DTSTEP
                0472 
                0473       DEFCIT = EVAP(CHNO) * ( RC(CHNO) + RA(CHNO) )
                0474 CCOOMMSTRDG1(CHNO) = DEFCIT / RA(CHNO)
                0475 CCOOMMSTRDG2(CHNO) = DEFCIT / ( RA(CHNO) + RCUN(CHNO) )
                0476 CCOOMMSTRDG3(CHNO) = DEFCIT / ( RA(CHNO) + RCUN(CHNO)/VPDSTR(CHNO) )
                0477 CCOOMMSTRDG4(CHNO) = DEFCIT / ( RA(CHNO) +
                0478 CCOOMM                   RCUN(CHNO)/(VPDSTR(CHNO)*FTEMP(CHNO)) )
                0479 
                0480  2000 CONTINUE
                0481 C****
                0482 
                0483       RETURN
                0484       END
                0485 C****
                0486 C**** [ END CHIP ]
                0487 C****
                0488 C**** -----------------------------------------------------------------
                0489 C**** /////////////////////////////////////////////////////////////////
                0490 C**** -----------------------------------------------------------------
                0491 C****
                0492 C**** [ BEGIN INTERC ]
aab9c05604 Jean*0493 C****
613fa3996d Andr*0494       SUBROUTINE INTERC (
                0495      I                   NCH, ITYP, DTSTEP, TRAINL, TRAINC,
                0496      I                   TSNOW, SATCAP, WMAX, CSOIL, WETEQ1,
                0497      U                   TC, CAPAC, SNOW,  SWET1,  RUNOFF, RUNSRF,
                0498      U                   SMELT, FWSOIL
                0499      &                  )
                0500 C****
aab9c05604 Jean*0501 C**** This routine uses the precipitation forcing to determine
613fa3996d Andr*0502 C**** changes in interception, snowcover, and soil moisture storage.
                0503 C****
                0504 C**** Note:  In this formulation, rain that falls on frozen ground
                0505 C**** runs-off, rather than freezes.
                0506 C****
                0507       IMPLICIT NONE
                0508 C****
4dad0f083a Andr*0509 #include "sibber.h"
613fa3996d Andr*0510 C****
                0511       INTEGER NCH
                0512       INTEGER ITYP(NCH), ChNo
aab9c05604 Jean*0513       _RL  TRAINL(NCH), TRAINC(NCH), TSNOW(NCH),  SATCAP(NCH),
                0514      &      WMAX(NLAY,NTYPS), TC(NCH),   CSOIL(NCH),   CAPAC(NCH),
613fa3996d Andr*0515      &      SNOW(NCH), SWET1(NCH),  RUNOFF(NCH),  SMELT(NCH),
                0516      &     RUNSRF(NCH), FWSOIL(NCH), WETEQ1(NCH), WETINT
a456aa407c Andr*0517       _RL DTSTEP, SNOWM, WATADD, CAVAIL, THRUC, WRUNC, WRUNL,
613fa3996d Andr*0518      &     TIMFRL, TIMFRC, FWETL, FWETC, THRU1, THRU2, THRUL, XTCORR,
                0519      &     WETFRC
                0520 C****
                0521 C      DATA FWET0  /0.30/
                0522       DATA FWETL /1.00/, FWETC /0.20/, TIMFRL/1.00/, TIMFRC/0.125/
                0523 C****
                0524 C**** ------------------------------------------------------------------
                0525 C**** Loop over chips:
                0526       DO 100 ChNo = 1, NCH
                0527 C****
                0528 C**** Add to snow cover.  Melt snow if necessary.
                0529       SNOW(ChNo) = SNOW(ChNo) + TSNOW(ChNo)*DTSTEP
                0530       SNOWM = 0.
                0531       IF( SNOW(CHNO).GT.0. .AND. TC(CHNO).GT.TF ) THEN
                0532         SNOWM = MIN( SNOW(ChNo),
32362b8fa8 Cons*0533      &       MAX( 0. _d 0, (TC(ChNo)-TF)*CSOIL(ChNo)/ALHM )  )
613fa3996d Andr*0534         IF( SNOWM .EQ. SNOW(CHNO) ) THEN
                0535             TC(ChNo) = TC(ChNo) - SNOWM * ALHM / CSOIL(ChNo)
                0536             SNOW(CHNO)=0.
                0537           ELSE
                0538             TC(CHNO)=TF
                0539             SNOW(ChNo) = SNOW(ChNo) - SNOWM
                0540           ENDIF
                0541         SMELT(CHNO)=SMELT(CHNO)+SNOWM/DTSTEP
                0542         ENDIF
                0543 C****
                0544 C**** =======================================================
                0545 C****
                0546 C**** Load interception reservoir.  STEP 1: Large scale condensation.
                0547 C****
                0548 C**** Determine XTCORR, the fraction of a storm that falls on a previously
                0549 C**** wet surface due to the time correlation of precipitation position.
                0550 C**** (Time scale TIMFRL for large scale storms set to one for FWETL=1
aab9c05604 Jean*0551 C**** to reflect the effective loss of "position memory" when storm
613fa3996d Andr*0552 C**** covers entire grid square.)
                0553 
e7070f537c Cons*0554         XTCORR= (1.-TIMFRL) * MIN( 1. _d 0,(CAPAC(CHNO)/SATCAP(CHNO))/
                0555      $       FWETL )
613fa3996d Andr*0556 
                0557 C****
                0558 C**** Fill interception reservoir with precipitation.
aab9c05604 Jean*0559 C**** THRU1 is first calculated as the amount falling through the
                0560 C****    canopy under the assumption that all rain falls randomly.
                0561 C****    only a fraction 1-XTCORR falls randomly, though, so the result
613fa3996d Andr*0562 C****    is multiplied by 1-XTCORR.
                0563 C****
                0564       WATADD = TRAINL(ChNo)*DTSTEP + SMELT(CHNO)*DTSTEP
                0565       CAVAIL = ( SATCAP(CHNO) - CAPAC(CHNO) ) * FWETL
                0566       WETINT = CAPAC(CHNO)/SATCAP(CHNO)
                0567       IF( WATADD*(1.-WETINT) .LT. CAVAIL ) THEN
                0568           THRU1 = WATADD*WETINT
                0569         ELSE
                0570           THRU1 = (WATADD - CAVAIL)
                0571         ENDIF
                0572       THRU1=THRU1*(1.-XTCORR)
                0573 
                0574 C**** THRU2 is the amount that falls immediately through the canopy due
                0575 C**** to 'position memory'.
                0576 
                0577       THRU2=XTCORR*WATADD
                0578 
                0579       THRUL=THRU1+THRU2
                0580       CAPAC(CHNO)=CAPAC(CHNO)+WATADD-THRU1-THRU2
                0581 C****
                0582 C**** ---------------------------------------------------
                0583 C****
                0584 C**** STEP 2: Moist convective precipitation.
                0585 C****
                0586 C**** Determine XTCORR, the fraction of a storm that falls on a previously
                0587 C**** wet surface due to the time correlation of precipitation position.
                0588 
e7070f537c Cons*0589       XTCORR= (1.-TIMFRC) * MIN( 1. _d 0,(CAPAC(CHNO)/SATCAP(CHNO))/
                0590      $     FWETC )
613fa3996d Andr*0591 
                0592 C****
                0593 C**** Fill interception reservoir with precipitation.
aab9c05604 Jean*0594 C**** THRU1 is first calculated as the amount falling through the
                0595 C****    canopy under the assumption that all rain falls randomly.
                0596 C****    only a fraction 1-XTCORR falls randomly, though, so the result
613fa3996d Andr*0597 C****    is multiplied by 1-XTCORR.
                0598 C****
                0599       WATADD = TRAINC(ChNo)*DTSTEP
                0600       CAVAIL = ( SATCAP(CHNO) - CAPAC(CHNO) ) * FWETC
                0601       WETINT = CAPAC(CHNO)/SATCAP(CHNO)
                0602       IF( WATADD*(1.-WETINT) .LT. CAVAIL ) THEN
                0603           THRU1 = WATADD*WETINT
                0604         ELSE
                0605           THRU1 = (WATADD - CAVAIL)
                0606         ENDIF
                0607       THRU1=THRU1*(1.-XTCORR)
                0608 
                0609 C**** THRU2 is the amount that falls immediately through the canopy due
                0610 C**** to 'position memory'.
                0611 
                0612       THRU2=XTCORR*WATADD
                0613 
                0614       THRUC=THRU1+THRU2
                0615       CAPAC(CHNO)=CAPAC(CHNO)+WATADD-THRU1-THRU2
                0616 C****
                0617 C***** =================================================================
                0618 C****
aab9c05604 Jean*0619 C**** Add precipitation moisture to soil, generate runoff.  if
613fa3996d Andr*0620 C**** temperature is below freezing, all precipitation runs off.
                0621 C****
                0622 C**** STEP 1: Large scale condensation:
                0623 C****
                0624       IF(SWET1(CHNO).GT.WETEQ1(CHNO) .AND. WETEQ1(CHNO).NE.1.) THEN
                0625           WETFRC=(SWET1(CHNO)-WETEQ1(CHNO))/(1.-WETEQ1(CHNO))
                0626         ELSE
                0627           WETFRC=0.
                0628         ENDIF
                0629 
                0630       CAVAIL = ( 1.-WETFRC)*FWETL*WMAX(1,ITYP(ChNo))*(1-WETEQ1(CHNO))
                0631       IF ( THRUL * (1-WETFRC) .LT. CAVAIL ) THEN
                0632           WRUNL = THRUL * WETFRC
aab9c05604 Jean*0633           SWET1(ChNo) = SWET1(ChNo)
613fa3996d Andr*0634      *    + THRUL * ( 1. - WETFRC) / WMAX(1,ITYP(ChNo))
                0635         ELSE
                0636           WRUNL = THRUL - CAVAIL
                0637           SWET1(CHNO) = SWET1(CHNO) + CAVAIL / WMAX(1,ITYP(ChNo))
                0638         ENDIF
                0639 C****
                0640 
                0641 C**** STEP 2: Moist convective precipitation:
                0642 C****
                0643       IF(SWET1(CHNO).GT.WETEQ1(CHNO) .AND. WETEQ1(CHNO).NE.1.) THEN
                0644           WETFRC=(SWET1(CHNO)-WETEQ1(CHNO))/(1.-WETEQ1(CHNO))
                0645         ELSE
                0646           WETFRC=0.
                0647         ENDIF
                0648 
                0649       CAVAIL = (1.-WETFRC)*FWETC*WMAX(1,ITYP(ChNo))*(1-WETEQ1(CHNO))
                0650       IF ( THRUC * (1-WETFRC) .LT. CAVAIL ) THEN
                0651           WRUNC = THRUC * WETFRC
aab9c05604 Jean*0652           SWET1(ChNo) = SWET1(ChNo)
613fa3996d Andr*0653      *    + THRUC * ( 1. - WETFRC) / WMAX(1,ITYP(ChNo))
                0654         ELSE
                0655           WRUNC = THRUC - CAVAIL
                0656           SWET1(CHNO) = SWET1(CHNO) + CAVAIL / WMAX(1,ITYP(ChNo))
                0657         ENDIF
                0658 C****
                0659       RUNOFF(ChNo) = RUNOFF(ChNo) + (WRUNC+WRUNL)/DTSTEP
                0660       RUNSRF(ChNo) = RUNSRF(ChNo) + (WRUNC+WRUNL)/DTSTEP
                0661       FWSOIL(CHNO) = FWSOIL(CHNO) + (THRUC+THRUL-WRUNC-WRUNL)/DTSTEP
                0662 C****
                0663  100  CONTINUE
                0664 C****
                0665       RETURN
                0666       END
                0667 C****
                0668 C**** [ END INTERC ]
                0669 C****
                0670 C**** -----------------------------------------------------------------
                0671 C**** /////////////////////////////////////////////////////////////////
                0672 C**** -----------------------------------------------------------------
                0673 C****
                0674 C**** [ BEGIN RCUNST ]
                0675 C****
                0676       SUBROUTINE RCUNST (
                0677      I                   NCH, ITYP, SUNANG, SQSCAT, PDIR,
                0678      I                   PAR, ZLAI, GREEN,
                0679      O                   RCUN
                0680      &                  )
                0681 C****
                0682 C****     This subroutine calculates the unstressed canopy resistance.
                0683 C**** (p. 1353, Sellers 1985.)  Extinction coefficients are computed first.
                0684 C****
                0685       IMPLICIT NONE
                0686 C****
4dad0f083a Andr*0687 #include "sibber.h"
613fa3996d Andr*0688 C****
                0689       INTEGER NCH
                0690       INTEGER ITYP(NCH), ChNo
                0691 
a456aa407c Andr*0692       _RL  SUNANG(NCH),   PDIR(NCH),   PAR(NCH),   ZLAI(NCH),
613fa3996d Andr*0693      &      SQSCAT(NCH),  GREEN(NCH),  RCUN(NCH)
                0694 
aab9c05604 Jean*0695       _RL   VGCHIL(NTYPS), VGZMEW(NTYPS),
613fa3996d Andr*0696      &       VGRST1(NTYPS), VGRST2(NTYPS), VGRST3(NTYPS)
                0697 
a456aa407c Andr*0698       _RL            RHO4,         EXTK1,         EXTK2,
613fa3996d Andr*0699      &               RCINV,         GAMMA,          EKAT,    DUM1,
                0700      &                DUM2,          DUM3,            AA,      BB,
                0701      &                  ZK,            CC
                0702 
                0703 
                0704       DATA VGCHIL /        0.1,        0.25,        0.01,        -0.3,
                0705      5                    0.01,        0.20,         0.0,         0.0,
                0706      9                     0.0,        0.0 /
                0707 
                0708       DATA VGZMEW/      0.9809,      0.9638,      0.9980,      1.0773,
                0709      5                  0.9980,      0.9676,       1.000,       1.000,
                0710      9                   1.000,      1.0000 /
                0711 
                0712       DATA VGRST1 /     2335.9,      9802.2,      2869.7,      2582.0,
                0713      5                 93989.4,      9802.2,         0.0,         0.0,
                0714      9                     0.0 ,        0.0/
                0715 
                0716       DATA VGRST2 /        0.0,        10.6,         3.7,         1.1,
                0717      5                    0.01,        10.6,         0.0,         0.0,
                0718      9                     0.0 ,        0.0/
                0719 
                0720       DATA VGRST3 /      153.5,       180.0,       233.0,       110.0,
                0721      5                   855.0,       180.0,         1.0,         1.0,
                0722      9                     1.0,         1.0 /
                0723 
                0724 
                0725 
                0726 
                0727       DO 100 ChNo = 1, NCH
                0728 
                0729 C**** First compute optical parameters.
                0730 C**** (Note: CHIL is constrained to be >= 0.01, as in SiB calcs.)
                0731 
                0732       AA = 0.5 - (0.633 + 0.330*VGCHIL(ITYP(ChNo)))*VGCHIL(ITYP(ChNo))
                0733       BB = 0.877 * ( ONE - 2.*AA )
                0734       CC =  ( AA + BB*SUNANG(ChNo) ) / SUNANG(ChNo)
                0735 
                0736       EXTK1 =  CC * SQSCAT(ChNo)
                0737       EXTK2 = (ONE / VGZMEW(ITYP(ChNo))) * SQSCAT(ChNo)
                0738 
                0739       DUM1 =      PDIR(ChNo)  *   CC
                0740       DUM2 = (ONE-PDIR(ChNo)) * ( BB*(ONE/3.+PIE/4.) + AA*1.5 )
                0741 
                0742 C**** Bound extinction coefficient by 50./ZLAI:
                0743 
32362b8fa8 Cons*0744       ZK =     PDIR(ChNo) *MIN( EXTK1, 50. _d 0/ZLAI(ChNo) ) +
                0745      &    (ONE-PDIR(ChNo))*MIN( EXTK2, 50. _d 0/ZLAI(ChNo) )
613fa3996d Andr*0746 
                0747 C**** Now compute unstressed canopy resistance:
                0748 
                0749       GAMMA = VGRST1(ITYP(ChNo)) / VGRST3(ITYP(ChNo)) +
                0750      &        VGRST2(ITYP(ChNo))
                0751 
                0752       EKAT = EXP( ZK*ZLAI(ChNo) )
                0753       RHO4 = GAMMA / (PAR(ChNo) * (DUM1 + DUM2))
                0754 
                0755       DUM1 = (VGRST2(ITYP(ChNo)) - GAMMA) / (GAMMA + 1.E-20)
                0756       DUM2 = (RHO4 * EKAT + ONE) / (RHO4 + ONE)
                0757       DUM3 = ZK * VGRST3(ITYP(ChNo))
                0758       RCINV = ( DUM1*LOG(DUM2) + ZK*ZLAI(ChNo) ) / DUM3
                0759 
                0760       RCUN(ChNo) = ONE / (RCINV * GREEN(ChNo) + 1.E-10)
                0761 
                0762  100  CONTINUE
                0763 
                0764 
                0765       RETURN
                0766       END
                0767 C****
                0768 C**** [ END RCUNST ]
                0769 C****
                0770 C**** -----------------------------------------------------------------
                0771 C**** /////////////////////////////////////////////////////////////////
                0772 C**** -----------------------------------------------------------------
                0773 C****
                0774 C**** [ BEGIN SOIL ]
                0775 C****
                0776       SUBROUTINE SOIL (
                0777      I         NCH, ITYP,  WET, VGBEEX, VGPSAX, VGCSAX, DELZ,
                0778      O         PHR, SOILCO, WETEQ
                0779      &                )
                0780 C****
                0781 C**** This subroutine returns soil moisture potential and conductivity.
                0782 
                0783       IMPLICIT NONE
                0784 C****
4dad0f083a Andr*0785 #include "sibber.h"
613fa3996d Andr*0786 C****
                0787       INTEGER NCH
                0788       INTEGER ITYP(NCH), ChNo
                0789 
a456aa407c Andr*0790       _RL    WET(NCH),    PHR(NCH), SOILCO(NCH), VGPSAX(NCH),
613fa3996d Andr*0791      &     VGCSAX(NCH), VGBEEX(NCH),  DELZ(NTYPS),  WETEQ(NCH),
                0792      &     WEXPB,  WET0,  PHEQ
                0793 
                0794       DO 100 ChNo = 1, NCH
                0795 
32362b8fa8 Cons*0796       WET0  = MAX(WET(CHNO),0.01 _d 0)
613fa3996d Andr*0797       WEXPB = WET0**VGBEEX(ChNo)
                0798 
                0799       PHR(ChNo)    = VGPSAX(ChNo) / WEXPB
                0800       SOILCO(ChNo) = VGCSAX(ChNo) * WEXPB * WEXPB * WET0 * WET0 * WET0
                0801 
                0802       PHEQ = PHR(CHNO) - DELZ(ITYP(CHNO))
                0803       WETEQ(CHNO) = ( PHEQ/VGPSAX(CHNO) ) ** ( -1/VGBEEX(CHNO) )
                0804 
                0805  100  CONTINUE
                0806 
                0807       RETURN
                0808       END
                0809 C****
                0810 C**** [ END SOIL ]
                0811 C****
                0812 C**** -----------------------------------------------------------------
                0813 C**** /////////////////////////////////////////////////////////////////
                0814 C**** -----------------------------------------------------------------
                0815 C****
                0816 C**** [ BEGIN FLUXES ]
                0817 C****
                0818       SUBROUTINE FLUXES (
                0819      I            NCH,   ITYP, DTSTEP, ESATTC, DESDTC,   ALHX,
                0820      I          ETURB,  DEDEA,  DEDTC, HSTURB, DHSDEA, DHSDTC,
                0821      I             RC, DRCDEA, DRCDTC,
                0822      I          SWNET, HLWDWN, ALWRAD, BLWRAD, ESNFRC,
                0823      I             TM,     EM,  CSOIL,   PSUR, EMAXRT,   WMAX,
                0824      U             TC,     TD,     EA,  SWET1,   SNOW,
                0825      O         RUNOFF,   EVAP, SHFLUX,  SMELT,  HLWUP, BOMB, STRDG1,
                0826      O         STRDG2, STRDG3, STRDG4, STRDG5, STRDG6, STRDG7,
                0827      O         STRDG8, STRDG9
                0828      &                  )
                0829 C****
                0830 C**** This subroutine computes the fluxes of latent and sensible heat
                0831 C**** from the surface through an energy balance calculation.
                0832 C****
                0833       IMPLICIT NONE
                0834 C****
4dad0f083a Andr*0835 #include "sibber.h"
613fa3996d Andr*0836 C****
                0837       INTEGER NCH
                0838       INTEGER ITYP(NCH), ChNo
                0839 
a456aa407c Andr*0840       _RL       DTSTEP, ESATTC(NCH), DESDTC(NCH),   ALHX(NCH),
613fa3996d Andr*0841      &       ETURB(NCH),  DEDEA(NCH),  DEDTC(NCH),
                0842      &      HSTURB(NCH), DHSDEA(NCH), DHSDTC(NCH),
                0843      &          RC(NCH), DRCDEA(NCH), DRCDTC(NCH),
                0844      &       SWNET(NCH), HLWDWN(NCH), ALWRAD(NCH), BLWRAD(NCH),
                0845      &          TM(NCH),     EM(NCH),  CSOIL(NCH),   PSUR(NCH),
                0846      &      EMAXRT(NCH),         WMAX(NLAY,NTYPS),
                0847      &          TC(NCH),     TD(NCH),     EA(NCH), ESNFRC(NCH),
                0848      &               SWET1(NCH,NLAY),   SNOW(NCH),
                0849      &      RUNOFF(NCH),   EVAP(NCH), SHFLUX(NCH),  SMELT(NCH),
                0850      &       HLWUP(NCH),   BOMB(NCH)
a456aa407c Andr*0851       _RL STRDG1(NCH),STRDG2(NCH),STRDG3(NCH),STRDG4(NCH)
                0852       _RL STRDG5(NCH),STRDG6(NCH),STRDG7(NCH),STRDG8(NCH)
                0853       _RL STRDG9(NCH)
613fa3996d Andr*0854 
a456aa407c Andr*0855       _RL     HLWTC,  CDEEPS,      Q0,  RHOAIR,   CONST,  DHLWTC,
613fa3996d Andr*0856      &        EPLANT,     A11,     A12,     A21,     A22,      F0,
                0857      &           DEA,     DTC,  SNLEFT,     Q0X,  Q0SNOW,
                0858      &         EANEW,  ESATNW,  EHARMN, DETERM, DENOM
                0859 
87329dfe43 Andr*0860       LOGICAL CHOKE
a456aa407c Andr*0861       _RL deepfac(ntyps)
613fa3996d Andr*0862       DATA deepfac /1.,1.,1.,1.,1.,1.,1.,1.,1.,1./
                0863 C****
                0864 C**** -------------------------------------------------------------------
                0865 
                0866       DO 200 ChNo = 1, NCH
                0867 C****
                0868       HLWTC = ALWRAD(CHNO) + BLWRAD(CHNO) * TC(CHNO)
                0869       CDEEPS = PIE * CSOIL(ChNo) * 2. / 86400.
                0870       RHOAIR = PSUR(ChNo) * 100. / (RGAS * TC(ChNo))
                0871       CONST = RHOAIR * EPSILON / PSUR(ChNo)
                0872       DHLWTC = BLWRAD(CHNO)
                0873 C****
                0874 C**** Compute matrix elements A11, A22, AND Q0 (energy balance equation).
                0875 C****
44d19a8651 Andr*0876 
613fa3996d Andr*0877       A11 = CSOIL(ChNo)/DTSTEP +
                0878      &        DHLWTC +
                0879      &        DHSDTC(ChNo) +
                0880      &        ALHX(CHNO)*DEDTC(ChNo) +
                0881      &        CDEEPS
                0882       A12 = DHSDEA(ChNo) + ALHX(CHNO) * DEDEA(ChNo)
                0883       Q0 =  SWNET(ChNo) +
                0884      &        HLWDWN(ChNo) -
                0885      &        HLWTC -
                0886      &        HSTURB(ChNo) -
                0887      &        ALHX(CHNO) * ETURB(ChNo) -
                0888      &        CDEEPS * (TC(ChNo) - TD(ChNo))
                0889 
                0890           STRDG3(ChNo)=Q0/A11
                0891           STRDG4(ChNo)=(SWNET(ChNo)+HLWDWN(ChNo)-HLWTC)/A11
                0892           STRDG5(ChNo)=HSTURB(ChNo)/A11
                0893           STRDG6(ChNo)=ALHX(CHNO)*ETURB(ChNo)/A11
                0894           STRDG7(ChNo)=CDEEPS*(TC(ChNo) - TD(ChNo))/A11
                0895           STRDG9(ChNo)=A11
                0896 
                0897 C****
aab9c05604 Jean*0898 C**** Compute matrix elements A21, A22, and F0 (canopy water budget
613fa3996d Andr*0899 C**** equation) and solve for fluxes.  Three cases are considered:
                0900 C****
                0901 C**** 1. Standard case: RC>0.
                0902 C**** 2. RC = 0.  Can only occur if CIR is full or ETURB is negative.
                0903 C****
                0904       CHOKE = .TRUE.
                0905 
aab9c05604 Jean*0906       IF( RC(CHNO) .GT. 0.) THEN
613fa3996d Andr*0907           EPLANT = CONST * (ESATTC(ChNo) - EA(ChNo)) / RC(ChNo)
                0908           IF(EPLANT*ETURB(ChNo).GE.0.) THEN
                0909               EHARMN = 2.*EPLANT*ETURB(CHNO) / (EPLANT + ETURB(ChNo))
                0910             ELSE
                0911               EHARMN=0.
                0912             ENDIF
                0913 C****
                0914 C****            Some limitations to A21 and A22 are applied:
                0915 C****            we assume that the increase in plant evaporation
aab9c05604 Jean*0916 C****            due to an increase in either TC or EA balances
613fa3996d Andr*0917 C****            or outweighs any decrease due to RC changes.
                0918 C****
                0919 
                0920           A21 =  -DEDTC(ChNo)*RC(ChNo) +
32362b8fa8 Cons*0921      &      max(0. _d 0, CONST*DESDTC(ChNo) - EHARMN*DRCDTC(ChNo) )
613fa3996d Andr*0922           A22 = -( RC(ChNo)*DEDEA(ChNo) +
32362b8fa8 Cons*0923      &               max( 0. _d 0, CONST + EHARMN*DRCDEA(ChNo) )   )
613fa3996d Andr*0924 
                0925           F0 = RC(ChNo) * (ETURB(ChNo) - EPLANT)
32362b8fa8 Cons*0926           DETERM = MIN( A12*A21/(A11*A22) - 1., -0.1 _d 0)
613fa3996d Andr*0927           DEA = ( Q0*A21 - A11*F0 ) / ( DETERM * A11*A22 )
                0928           DTC = ( Q0 - A12*DEA ) / A11
aab9c05604 Jean*0929 
613fa3996d Andr*0930           STRDG1(ChNo)=DTC
                0931           STRDG2(ChNo)=DEA
                0932           STRDG8(ChNo)=-A12*DEA/A11
                0933 
                0934           EVAP(ChNo) = ETURB(ChNo) + DEDEA(ChNo)*DEA + DEDTC(ChNo)*DTC
                0935           SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                0936      &                                            DHSDTC(ChNo)*DTC
                0937           DENOM = DETERM * A11*A22
                0938         ELSE
                0939           CHOKE = .FALSE.
                0940           A21 = -DESDTC(ChNo)
                0941           A22 = 1.
                0942           F0 = ESATTC(ChNo) - EA(ChNo)
                0943           DEA = ( Q0*A21 - A11*F0 ) / ( A12*A21 - A11*A22 )
                0944           DTC = ( Q0 - A12*DEA ) / A11
                0945 
                0946           STRDG1(ChNo)=DTC
                0947           STRDG2(ChNo)=DEA
                0948           STRDG8(ChNo)=-A12*DEA/A11
                0949 
                0950           EVAP(ChNo) = ETURB(ChNo) + DEDEA(ChNo)*DEA + DEDTC(ChNo)*DTC
                0951           SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                0952      &                                            DHSDTC(ChNo)*DTC
                0953           DENOM = A12 * A21 - A11*A22
                0954         ENDIF
                0955 
                0956 
                0957 C**** - - - - - - - - - - - - - - - - - - - -
                0958 C**** Account for snowmelt, if necessary.
                0959 C**** - - - - - - - - - - - - - - - - - - - -
                0960       SMELT(CHNO)=0.
                0961       SNLEFT=SNOW(CHNO)-EVAP(CHNO)*DTSTEP*ESNFRC(CHNO)
                0962       IF(TC(CHNO)+DTC .GT. TF .AND. SNLEFT.GT.0.) THEN
                0963 
                0964 C****      First re-calculate energy balance under assumption that all
                0965 C****      snow is melted.
                0966 
                0967         Q0X=Q0-ALHM*SNLEFT/DTSTEP
                0968         SMELT(CHNO)=SNLEFT/DTSTEP
aab9c05604 Jean*0969 
613fa3996d Andr*0970             DEA = ( Q0X*A21 - A11*F0 ) / DENOM
                0971             DTC = ( Q0X - A12*DEA ) / A11
                0972 
                0973 C****      If TC+DTC is now less than TF, too much snow has melted.  Now
                0974 C****      solve for balance assuming only some of the snow has melted.
                0975 C****      Set TC to TF.
                0976 
                0977         IF(TC(CHNO)+DTC .LT. TF) THEN
                0978           DTC=TF-TC(CHNO)
                0979           DEA=(F0-A21*DTC)/A22
                0980           Q0SNOW=A11*DTC+A12*DEA
                0981           SMELT(CHNO)=(Q0-Q0SNOW)/ALHM
                0982           ENDIF
                0983 
                0984           STRDG1(ChNo)=DTC
                0985           STRDG2(ChNo)=DEA
                0986           STRDG8(ChNo)=-A12*DEA/A11
aab9c05604 Jean*0987 
613fa3996d Andr*0988         EVAP(ChNo) = ETURB(ChNo) + DEDEA(ChNo)*DEA +
                0989      &                   DEDTC(ChNo)*DTC
                0990         SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                0991      &                                            DHSDTC(ChNo)*DTC
                0992 
                0993         ENDIF
                0994 
                0995 C**** - - - - - - - - - - - - - - - - - - - - - - -
                0996 C**** Adjustments
                0997 
                0998 C**** 1. Adjust deltas and fluxes if all available water evaporates
                0999 C****    during time step:
                1000 C****
                1001       IF( EVAP(CHNO) .GT. EMAXRT(CHNO) ) THEN
                1002         CHOKE = .FALSE.
                1003         DEA = EM(CHNO) - EA(CHNO)
aab9c05604 Jean*1004         DTC =
4dad0f083a Andr*1005      &   (Q0 + ALHX(CHNO)*(ETURB(ChNo)-EMAXRT(CHNO)) - DHSDEA(CHNO)*DEA)
613fa3996d Andr*1006      &          /  ( A11 - ALHX(CHNO)*DEDTC(ChNo) )
aab9c05604 Jean*1007 
613fa3996d Andr*1008           STRDG1(ChNo)=DTC
                1009           STRDG2(ChNo)=DEA
                1010           STRDG8(ChNo)=0.
aab9c05604 Jean*1011 
613fa3996d Andr*1012         EVAP(CHNO) = EMAXRT(CHNO)
                1013         SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                1014      &                                            DHSDTC(ChNo)*DTC
                1015         SMELT(CHNO)=0.
                1016         ENDIF
                1017 C****
                1018 C**** Adjust DEA and DTC if solutions were pathological:
                1019 C****
                1020       ESATNW = ESATTC(CHNO)+DESDTC(CHNO)*DTC
                1021       EANEW = EA(CHNO) + DEA
                1022 
                1023 
                1024 
aab9c05604 Jean*1025 C**** 2. Pathological cases.
613fa3996d Andr*1026 
                1027 C**** Case 1: EVAP is positive in presence of negative gradient.
aab9c05604 Jean*1028 C**** Case 2: EVAP and ETURB have opposite sign, implying that
613fa3996d Andr*1029 C**** "virtual effects" derivatives are meaningless and thus that we
3daafce20b Jean*1030 C**** do not know the proper tendency terms.
613fa3996d Andr*1031 C**** In both cases, assume zero evaporation for the time step.
                1032 
                1033 C      IF( ( EVAP(CHNO) .LT. 0. .AND. EM(CHNO).LT.ESATNW )
aab9c05604 Jean*1034 C     &    .OR. ( EVAP(CHNO)*ETURB(CHNO) .LT. 0. ) )  THEN
                1035       IF( EVAP(CHNO) .LT. 0. .AND. EM(CHNO).LT.ESATNW )  THEN
613fa3996d Andr*1036         CHOKE = .FALSE.
                1037         DEA = EM(CHNO) - EA(CHNO)
aab9c05604 Jean*1038         DTC = ( Q0 + ALHX(CHNO)*ETURB(ChNo) - DHSDEA(CHNO)*DEA ) /
613fa3996d Andr*1039      &            ( A11 - ALHX(CHNO)*DEDTC(ChNo) )
aab9c05604 Jean*1040 
613fa3996d Andr*1041           STRDG1(ChNo)=DTC
                1042           STRDG2(ChNo)=DEA
                1043           STRDG8(ChNo)=0.
                1044 
                1045         EVAP(CHNO) = 0.
                1046         SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                1047      &                                            DHSDTC(ChNo)*DTC
                1048 
                1049         IF(TC(CHNO)+DTC .GT. TF .AND. SNOW(CHNO).GT.0.) THEN
                1050           Q0X=Q0-ALHM*SNOW(CHNO)/DTSTEP
                1051           SMELT(CHNO)=SNOW(CHNO)/DTSTEP
aab9c05604 Jean*1052           DTC = ( Q0X + ALHX(CHNO)*ETURB(ChNo) - DHSDEA(CHNO)*DEA ) /
613fa3996d Andr*1053      &              ( A11 - ALHX(CHNO)*DEDTC(ChNo) )
                1054           IF(TC(CHNO)+DTC .LT. TF) THEN
                1055             DTC=TF-TC(CHNO)
                1056             Q0SNOW=A11*DTC+A12*DEA
                1057             SMELT(CHNO)=(Q0-Q0SNOW)/ALHM
                1058             ENDIF
                1059           SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                1060      &                                            DHSDTC(ChNo)*DTC
                1061           ENDIF
                1062         ENDIF
                1063 
                1064 
                1065 
                1066 C**** 3. Exceesive dea change: apply "choke".
                1067 
                1068       IF( CHOKE .AND. ABS(DEA) .GT. 0.5*EA(CHNO) ) THEN
                1069         DEA = SIGN(.5*EA(CHNO),DEA)
                1070         DTC = ( Q0 - A12*DEA ) / A11
aab9c05604 Jean*1071 
613fa3996d Andr*1072           STRDG1(ChNo)=DTC
                1073           STRDG2(ChNo)=DEA
                1074           STRDG8(ChNo)=-A12*DEA/A11
                1075 
                1076         EVAP(ChNo) = ETURB(ChNo) + DEDEA(ChNo)*DEA + DEDTC(ChNo)*DTC
                1077         SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                1078      &                                            DHSDTC(ChNo)*DTC
                1079         IF(TC(CHNO)+DTC .GT. TF .AND. SNOW(CHNO).GT.0.) THEN
                1080           SNLEFT=SNOW(CHNO)-EVAP(CHNO)*DTSTEP*ESNFRC(CHNO)
                1081           Q0X=Q0-ALHM*SNLEFT/DTSTEP
                1082           SMELT(CHNO)=SNLEFT/DTSTEP
                1083           DTC = ( Q0X - A12*DEA ) / A11
                1084           IF(TC(CHNO)+DTC .LT. TF) THEN
                1085             DTC=TF-TC(CHNO)
                1086             Q0SNOW=A11*DTC+A12*DEA
                1087             SMELT(CHNO)=(Q0-Q0SNOW)/ALHM
                1088             ENDIF
                1089           EVAP(ChNo) = ETURB(ChNo) + DEDEA(ChNo)*DEA + DEDTC(ChNo)*DTC
                1090           SHFLUX(ChNo) = HSTURB(ChNo) + DHSDEA(ChNo)*DEA +
                1091      &                                            DHSDTC(ChNo)*DTC
                1092           ENDIF
                1093         ENDIF
                1094 
                1095 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                1096 
                1097       TC(ChNo) = TC(ChNo) + DTC
                1098       EA(ChNo) = EA(ChNo) + DEA
                1099       TD(ChNo) = TD(ChNo) +
                1100 c  CHANGED THIS: deep layer 2 times deeper
aab9c05604 Jean*1101      &     DTSTEP * CDEEPS * (TC(ChNo) - TD(ChNo)) /
613fa3996d Andr*1102      .       (CSOIL(ChNo)*67.73*deepfac(ityp(ChNo)))
                1103       HLWUP(CHNO) = HLWTC + DHLWTC*DTC
                1104 
                1105       SNOW(CHNO)=SNOW(CHNO)-SMELT(CHNO)*DTSTEP
                1106 
                1107 C**** Make sure EA remains positive
                1108 
32362b8fa8 Cons*1109       EA(CHNO) = MAX(EA(CHNO), 0.0 _d 0)
613fa3996d Andr*1110 
                1111   200 CONTINUE
                1112 
                1113 
                1114 
                1115       RETURN
                1116       END
                1117 C****
                1118 C**** [ END FLUXES ]
                1119 C****
                1120 C**** -----------------------------------------------------------------
                1121 C**** /////////////////////////////////////////////////////////////////
                1122 C**** -----------------------------------------------------------------
                1123 C****
                1124 C**** [ BEGIN VPDFAC ]
                1125 C****
                1126       SUBROUTINE VPDFAC (
                1127      I                   NCH, ITYP, ESATTC, EA,
                1128      O                   VPDSTR
                1129      &                  )
                1130 C****
                1131 C**** This subroutine computes the vapor pressure deficit stress.
                1132 C****
                1133       IMPLICIT NONE
                1134 C****
4dad0f083a Andr*1135 #include "sibber.h"
613fa3996d Andr*1136 C****
                1137       INTEGER NCH
                1138       INTEGER ITYP(NCH), ChNo
a456aa407c Andr*1139       _RL    ESATTC(NCH), EA(NCH),  VPDSTR(NCH)
                1140 c     _RL  VGDFAC(NTYPS)
613fa3996d Andr*1141 C****
87329dfe43 Andr*1142 c     DATA VGDFAC /   .0273,    .0357,    .0310,    .0238,
                1143 c    5                .0275,    .0275,       0.,       0.,
                1144 c    9                   0.,      0. /
613fa3996d Andr*1145 C****
                1146 C**** -----------------------------------------------------------------
                1147 
                1148       DO 100 ChNo = 1, NCH
                1149 C****
                1150 c      VPDSTR(ChNo) = 1. - (ESATTC(ChNo)-EA(ChNo)) * VGDFAC(ITYP(ChNo))
32362b8fa8 Cons*1151 c      VPDSTR (ChNo) = MIN( 1. _d 0, MAX( VPDSTR(ChNo), 1. _d -10 ) )
613fa3996d Andr*1152       VPDSTR(CHNO) = 1.
                1153 C****
                1154  100  CONTINUE
                1155 C****
                1156       RETURN
                1157       END
                1158 C****
                1159 C**** [ END VPDFAC ]
                1160 C****
                1161 C**** -----------------------------------------------------------------
                1162 C**** /////////////////////////////////////////////////////////////////
                1163 C**** -----------------------------------------------------------------
                1164 C****
                1165 C**** [ BEGIN TMPFAC ]
                1166 C****
                1167       SUBROUTINE TMPFAC (
                1168      I                   NCH,  ITYP, TC,
                1169      O                   FTEMP
                1170      &                  )
                1171 C****
                1172 C**** Compute temperature stress factor.
                1173 C****
                1174       IMPLICIT NONE
                1175 C****
4dad0f083a Andr*1176 #include "sibber.h"
613fa3996d Andr*1177 C****
                1178       INTEGER NCH
                1179       INTEGER ITYP(NCH), ChNo, TypPtr
a456aa407c Andr*1180       _RL TC(NCH), FTEMP(NCH)
                1181       _RL    VGTLL(MemFac*NTYPS),    VGTU(MemFac*NTYPS),
613fa3996d Andr*1182      &       VGTCF1(MemFac*NTYPS),  VGTCF2(MemFac*NTYPS),
                1183      &       VGTCF3(MemFac*NTYPS)
                1184 C****
                1185       DATA VGTLL /MemFac*273., MemFac*273., MemFac*268., MemFac*283.,
aab9c05604 Jean*1186      5            MemFac*283., MemFac*273., MemFac*  0., MemFac*  0.,
613fa3996d Andr*1187      9            MemFac*  0., MemFac*  0. /
                1188       DATA VGTU /MemFac*318., MemFac*318., MemFac*313., MemFac*328.,
                1189      5           MemFac*323., MemFac*323., MemFac*  0., MemFac*  0.,
                1190      9           MemFac*  0. , MemFac*  0./
                1191       DATA VGTCF1 / MemFac*-1.43549E-06,  MemFac*-6.83584E-07,
                1192      3              MemFac* 1.67699E-07,  MemFac*-1.43465E-06,
                1193      5              MemFac*-2.76097E-06,  MemFac*-1.58094E-07,
                1194      7                       MemFac* 0.,           MemFac* 0.,
                1195      9                       MemFac* 0. , MemFac*  0./
                1196       DATA VGTCF2 / MemFac* 7.95859E-04,  MemFac* 3.72064E-04,
                1197      3              MemFac*-7.65944E-05,  MemFac* 8.24060E-04,
aab9c05604 Jean*1198      5              MemFac* 1.57617E-03,  MemFac* 8.44847E-05,
613fa3996d Andr*1199      7                       MemFac* 0.,           MemFac* 0.,
                1200      9                       MemFac* 0. , MemFac*  0./
                1201       DATA VGTCF3 / MemFac*-1.11575E-01,  MemFac*-5.21533E-02,
                1202      3              MemFac* 6.14960E-03,  MemFac*-1.19602E-01,
                1203      5              MemFac*-2.26109E-01,  MemFac*-1.27272E-02,
                1204      7                       MemFac* 0.,           MemFac* 0.,
                1205      9                       MemFac* 0. , MemFac*  0./
                1206 C****
                1207 C**** ----------------------------------------------------------------
                1208 
                1209       DO 100 ChNo = 1, NCH
                1210 C****
                1211       TypPtr = MOD(ChNo,MemFac) + (ITYP(ChNo)-1)*MemFac + 1
                1212       FTEMP(ChNo) = (TC(ChNo) - VGTLL(TypPtr)) *
                1213      &              (TC(ChNo) - VGTU(TypPtr)) *
                1214      &                    ( VGTCF1(TypPtr)*TC(ChNo)*TC(ChNo) +
                1215      &                      VGTCF2(TypPtr)*TC(ChNo) +
                1216      &                      VGTCF3(TypPtr) )
                1217       IF ( TC(ChNo) .LE. VGTLL(TypPtr) .OR. TC(ChNo) .GE. VGTU(TypPtr) )
                1218      &      FTEMP (ChNo) = 1.E-10
32362b8fa8 Cons*1219       FTEMP(CHNO) = MIN( 1. _d 0, MAX( FTEMP(ChNo), 1. _d -10 ) )
613fa3996d Andr*1220 C****
                1221  100  CONTINUE
                1222 C****
                1223       RETURN
                1224       END
                1225 C****
                1226 C**** [ END TMPFAC ]
                1227 C****
                1228 C**** -----------------------------------------------------------------
                1229 C**** /////////////////////////////////////////////////////////////////
                1230 C**** -----------------------------------------------------------------
                1231 C****
                1232 C**** [ BEGIN SMFAC ]
                1233 C****
                1234       SUBROUTINE SMFAC (
                1235      I                  NCH,    ITYP,   ESATTC, EA, PHR,  SOILCO,
                1236      I                  RCUN,   VPDSTR, FTEMP,  TC, PSUR, Z2,
                1237      I                  RSOIL1, RSOIL2, VGPH1X, VGPH2X, VGRPLX,
                1238      O                  RC
                1239      &                 )
                1240 C****
aab9c05604 Jean*1241 C**** This subroutine estimates RC after computing the
613fa3996d Andr*1242 C**** leaf water potential stress.
                1243 C****
                1244       IMPLICIT NONE
                1245 C****
4dad0f083a Andr*1246 #include "sibber.h"
613fa3996d Andr*1247 C****
                1248       INTEGER NCH
                1249       INTEGER ITYP(NCH), ChNo
a456aa407c Andr*1250       _RL  ESATTC(NCH),      EA(NCH),     PHR(NCH),  SOILCO(NCH),
613fa3996d Andr*1251      &        RCUN(NCH),  VPDSTR(NCH),   FTEMP(NCH),      TC(NCH),
                1252      &        PSUR(NCH),      Z2(NCH),  RSOIL1(NCH),  RSOIL2(NCH),
                1253      &      VGPH1X(NCH),  VGPH2X(NCH),  VGRPLX(NCH),      RC(NCH)
a456aa407c Andr*1254       _RL   RCUNTD,  RHOAIR,   CONST,     DEF,     D12,     DR2,
613fa3996d Andr*1255      &        RSOIL,      R0,    EEST,  RSTFAC
                1256 C****
                1257 C**** -----------------------------------------------------------------
                1258 
                1259       DO 100 ChNo = 1, NCH
                1260 C****
                1261       RCUNTD = RCUN(ChNo) / ( VPDSTR(ChNo)*FTEMP(ChNo) )
                1262       RHOAIR = PSUR(ChNo) * 100. / ( RGAS*TC(ChNo) )
                1263 
                1264 C****
                1265       CONST = RHOAIR * EPSILON / PSUR(ChNo)
                1266       DEF = ( ESATTC(ChNo) - EA(ChNo) ) * CONST
                1267       D12 = VGPH1X(ChNo) - VGPH2X(ChNo)
                1268       DR2 = PHR(ChNo) - Z2(ChNo) - VGPH2X(ChNo)
                1269       RSOIL = RSOIL1(ChNo) + RSOIL2(ChNo)/SOILCO(ChNo)
                1270       R0 = ( VGRPLX(ChNo) + RSOIL ) / RHOW
                1271       EEST = DEF*DR2 / ( RCUNTD*D12 + DEF*R0 )
                1272       RSTFAC = ( DR2 - R0*EEST ) / D12
32362b8fa8 Cons*1273       RSTFAC = MIN( 1. _d 0, MAX( 0.001 _d 0, RSTFAC ) )
613fa3996d Andr*1274       RC(ChNo) = RCUNTD / RSTFAC
                1275 C****
                1276  100  CONTINUE
                1277 C****
                1278       RETURN
                1279       END
                1280 C****
                1281 C**** [ END SMFAC ]
                1282 C****
                1283 C**** -----------------------------------------------------------------
                1284 C**** /////////////////////////////////////////////////////////////////
                1285 C**** -----------------------------------------------------------------
                1286 C****
                1287 C**** [ BEGIN RSURFP ]
                1288 C****
                1289       SUBROUTINE RSURFP (
                1290      I                   NCH, UM, U2FAC, Z2, RDC, WET, ESATTC, EA,
                1291      U                   RC,
                1292      O                   RX1, RX2
                1293      &                  )
                1294 C****
                1295       IMPLICIT NONE
                1296       INTEGER NCH
                1297       INTEGER ChNo
aab9c05604 Jean*1298       _RL     UM(NCH),   U2FAC(NCH),     Z2(NCH),    RDC(NCH),
613fa3996d Andr*1299      &        WET(NCH),  ESATTC(NCH),     EA(NCH),
                1300      &         RC(NCH),     RX1(NCH),    RX2(NCH)
a456aa407c Andr*1301       _RL  U2, RSURF, HESAT
613fa3996d Andr*1302 C****
                1303 C**** -----------------------------------------------------------------
                1304 
                1305       DO 100 ChNo = 1, NCH
                1306 C****
                1307       U2 = UM(ChNo) * U2FAC(ChNo)
                1308 c      RSURF = RDC(ChNo) / U2 + 30. / (1.E-20 + WET(ChNo))
                1309       RSURF = RDC(ChNo) / U2 + 26. + 6. / (1.E-20 + WET(ChNo))**2
                1310 
                1311 C**** Account for subsaturated humidity at soil surface:
                1312 C****
32362b8fa8 Cons*1313       HESAT = ESATTC(CHNO) * MIN( 1. _d 0, WET(CHNO)*2. _d 0)
613fa3996d Andr*1314       IF( EA(CHNO) .LT. HESAT ) THEN
                1315           RSURF=RSURF*( 1. + (ESATTC(CHNO)-HESAT)/(HESAT-EA(CHNO)) )
                1316         ELSE
                1317           RSURF=1.E10
                1318         ENDIF
                1319 
                1320 
                1321       RX1(CHNO)=RC(CHNO)
                1322       RX2(CHNO)=RSURF
                1323 
                1324       RC(ChNo) = RC(CHNO) * RSURF / ( RC(ChNo) + RSURF )
                1325 C****
                1326  100  CONTINUE
                1327 C****
                1328       RETURN
                1329       END
                1330 C****
                1331 C**** [ END RSURFP ]
                1332 C****
                1333 C**** -----------------------------------------------------------------
                1334 C**** /////////////////////////////////////////////////////////////////
                1335 C**** -----------------------------------------------------------------
                1336 C****
                1337 C**** [ BEGIN RCANOP ]
                1338 C****
                1339       SUBROUTINE RCANOP (
                1340      I                   NCH, CAPAC, SNOW, SATCAP, RA, ETURB, SNWFRC,
                1341      I                   POTFRC,
                1342      U                   RC
                1343      &                  )
                1344 C****
aab9c05604 Jean*1345 C**** The effective latent heat resistance RC depends on the quantity
613fa3996d Andr*1346 C**** of interception reservoir water and the snow cover.  POTFRC
                1347 C**** is the fraction of the tile from which potential evaporation
                1348 C**** occurs.
                1349 C****
                1350       IMPLICIT NONE
                1351       INTEGER NCH
                1352       INTEGER ChNo
a456aa407c Andr*1353       _RL CAPAC(NCH), SNOW(NCH), SATCAP(NCH), RA(NCH), ETURB(NCH),
613fa3996d Andr*1354      &          RC(NCH), SNWFRC(NCH), POTFRC(NCH)
a456aa407c Andr*1355       _RL ETCRIT,RAMPFC
613fa3996d Andr*1356 
                1357 C**** (Note: ETCRIT arbitrarily set to ~-5 W/m2, or -2.e-6 mm/sec.)
                1358       DATA ETCRIT/ -2.E-6 /
                1359 C****
                1360 C**** -----------------------------------------------------------------
                1361 
                1362       DO 100 ChNo = 1, NCH
                1363 
aab9c05604 Jean*1364 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
613fa3996d Andr*1365 C**** Case 1: Vegetation present (SATCAP GT .001).  Potential evap. from
                1366 c**** both interception reservoir and snow.
                1367 
                1368       IF(SATCAP(CHNO).GT..001) THEN
                1369         RC(ChNo)=RC(ChNo)*(1.-POTFRC(CHNO))/
                1370      &                ( 1.+POTFRC(CHNO)*RC(ChNo)/RA(ChNo) )
                1371         ENDIF
                1372 
aab9c05604 Jean*1373 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
613fa3996d Andr*1374 C**** Case 2: Vegetation absent (SATCAP LE .001).  Potential evap. from
                1375 c**** snow only.
                1376 
                1377       IF(SATCAP(CHNO) .LE. .001) THEN
                1378         RC(ChNo)=RC(ChNo)*(1.-SNWFRC(CHNO))/
                1379      &                  ( 1.+SNWFRC(CHNO)*RC(ChNo)/RA(ChNo) )
                1380         ENDIF
                1381 
                1382 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                1383 C**** Assume RC=0 for condensation (dew).
                1384 C**** RAMPFC is used to ensure continuity in RC.
                1385 
                1386       RAMPFC=ETURB(CHNO)/ETCRIT
                1387       IF ( RAMPFC .GE. 0. ) RC(CHNO) = RC(CHNO)*(1.-RAMPFC)
                1388       IF ( RAMPFC .GT. 1. ) RC(ChNo) = 0.
                1389 C****
                1390  100  CONTINUE
                1391 C****
                1392       RETURN
                1393       END
                1394 C****
                1395 C**** [ END RCANOP ]
                1396 C****
                1397 C**** -----------------------------------------------------------------
                1398 C**** /////////////////////////////////////////////////////////////////
                1399 C**** -----------------------------------------------------------------
                1400 C****
                1401 C***** [ BEGIN WUPDAT ]
                1402 C****
                1403       SUBROUTINE WUPDAT (
                1404      I                     NCH,   ITYP, DTSTEP,
                1405      I                     EVAP, SATCAP, VGWMAX, TC, RA, RC,
                1406      I                     RX1, RX2,ESNFRC,EIRFRC,
                1407      U                     CAPAC, SNOW, SWET, RUNOFF,
                1408      O                     EINT, ESOI, EVEG, ESNO, RUNSRF, FWSOIL
                1409      &                    )
                1410 C****
                1411 C**** This subroutine allows evapotranspiration to adjust the water
                1412 C**** contents of the interception reservoir and the soil layers.
                1413 C****
                1414       IMPLICIT NONE
                1415 C****
4dad0f083a Andr*1416 #include "sibber.h"
613fa3996d Andr*1417 C****
                1418       INTEGER NCH
                1419       INTEGER ITYP(NCH), ChNo
a456aa407c Andr*1420       _RL    EVAP(NCH),  SATCAP(NCH),  VGWMAX(NLAY,NTYPS),
613fa3996d Andr*1421      &          TC(NCH),      RA(NCH),             RC(NCH),
                1422      &       CAPAC(NCH),    SNOW(NCH),    SWET(nch,NLAY),
                1423      &      RUNOFF(NCH),     RX1(NCH),
                1424      &         RX2(NCH),  RUNSRF(NCH),  FWSOIL(NCH),
                1425      &      ESNFRC(NCH),  EIRFRC(NCH)
a456aa407c Andr*1426       _RL EINT(NCH), ESOI(NCH), EVEG(NCH), ESNO(NCH)
                1427       _RL  DTSTEP, EGRO, FWS, THRU, DEWRUN,
613fa3996d Andr*1428      &     WTOTAL,WLAY1,WLAY2,ELAY1,ELAY2,EGROI
                1429 C****
                1430 C**** -----------------------------------------------------------------
                1431 
                1432       DO 100 ChNo = 1, NCH
                1433 C****
                1434 C**** Partition evap between interception, snow, and ground reservoirs.
                1435 C****
                1436       WLAY1 = SWET(ChNo,SFCLY) * VGWMAX(SFCLY,ITYP(ChNo))
                1437       WLAY2 = SWET(ChNo,ROOTLY) * VGWMAX(ROOTLY,ITYP(ChNo))
                1438       WTOTAL = WLAY1 + WLAY2
                1439 C****
                1440       ESNO(CHNO)=ESNFRC(CHNO)*EVAP(CHNO)*DTSTEP
                1441       EINT(CHNO)=EIRFRC(CHNO)*EVAP(CHNO)*DTSTEP
                1442       EGRO = EVAP(CHNO)*DTSTEP - ESNO(CHNO) - EINT(CHNO)
                1443 
                1444 C**** Ensure that individual capacities are not exceeded.
                1445 
                1446       IF(ESNO(CHNO) .GT. SNOW(CHNO)) THEN
                1447         EINT(CHNO)=EINT(CHNO)+(ESNO(CHNO)-SNOW(CHNO))
                1448         ESNO(CHNO)=SNOW(CHNO)
                1449         ENDIF
                1450       EGROI=EGRO+EINT(CHNO)
                1451       IF(EGROI .GT. CAPAC(CHNO)+WTOTAL) THEN
                1452         ESNO(CHNO)=ESNO(CHNO)+EGROI-(CAPAC(CHNO)+WTOTAL)
                1453         EGROI=CAPAC(CHNO)+WTOTAL
                1454         ENDIF
                1455 
                1456       EINT(CHNO)=EGROI-EGRO
                1457       IF(EINT(CHNO) .GT. CAPAC(CHNO)) THEN
                1458         EGRO=EGRO+EINT(CHNO)-CAPAC(CHNO)
                1459         EINT(CHNO)=CAPAC(CHNO)
                1460         ENDIF
                1461       IF(EGRO .GT. WTOTAL) THEN
                1462         EINT(CHNO)=EINT(CHNO)+EGRO-WTOTAL
                1463         EGRO=WTOTAL
                1464         ENDIF
                1465 C****
                1466 C**** Separate egro into surface-evaporation/transpiration components:
                1467 C****
                1468       IF( RX1(CHNO)+RX2(CHNO) .NE. 0. ) THEN
                1469           ESOI(CHNO)=EGRO*RX1(CHNO)/(RX1(CHNO)+RX2(CHNO))
                1470           EVEG(CHNO)=EGRO - ESOI(CHNO)
                1471         ELSE
                1472           ESOI(CHNO)=EGRO/2.
                1473           EVEG(CHNO)=EGRO/2.
                1474         ENDIF
                1475 C****
                1476 C**** Translate ESOI and EVEG into evaporation fluxes from layers 1 and 2:
                1477 C****
                1478       IF( WTOTAL .GT. 0. ) THEN
aab9c05604 Jean*1479           FWS = EVEG(CHNO) / WTOTAL
613fa3996d Andr*1480         ELSE
                1481           FWS = 1.
                1482         ENDIF
                1483       ELAY1 = WLAY1*FWS + ESOI(CHNO)
                1484       ELAY2 = WLAY2*FWS
                1485 C****
                1486 C**** Ensure that enough soil water is available in each layer:
                1487 C****
                1488       IF( ELAY1 .GT. WLAY1 ) THEN
                1489         ELAY2 = ELAY2 + (ELAY1 - WLAY1)
                1490         ELAY1 = WLAY1
                1491         ENDIF
                1492       IF( ELAY2 .GT. WLAY2 ) THEN
                1493         ELAY1 = ELAY1 + (ELAY2 - WLAY2)
                1494         ELAY2 = WLAY2
                1495         ENDIF
                1496 
                1497 C**** Special case for condensation:
                1498       IF(EVAP(CHNO) .LT. 0.) THEN
                1499           EINT(CHNO)=(1.-ESNFRC(CHNO))*EVAP(CHNO)*DTSTEP
                1500           ELAY1=0.
                1501           ELAY2=0.
                1502           ESOI(CHNO)=0.
                1503           EVEG(CHNO)=0.
                1504           ESNO(CHNO)=ESNFRC(CHNO)*EVAP(CHNO)*DTSTEP
                1505         ENDIF
                1506 
                1507 C****
                1508 C**** Remove moisture from reservoirs:
                1509 C****
                1510       SNOW(ChNo) = SNOW(ChNo) - ESNO(CHNO)
                1511       CAPAC(ChNo) = CAPAC(ChNo) - EINT(CHNO)
aab9c05604 Jean*1512       SWET(ChNo,SFCLY) = (WLAY1 - ELAY1) / VGWMAX(SFCLY,ITYP(ChNo))
613fa3996d Andr*1513       SWET(ChNo,ROOTLY) = (WLAY2 - ELAY2) / VGWMAX(ROOTLY,ITYP(CHNO))
                1514 C****
                1515 C**** Ensure against numerical precision problems:
e7070f537c Cons*1516       SWET(ChNo,SFCLY)  = MIN( 1. _d 0, MAX( 0. _d 0, SWET(ChNo,SFCLY) )
                1517      $     )
                1518       SWET(ChNo,ROOTLY) = MIN( 1. _d 0, MAX( 0. _d 0, SWET(ChNo,ROOTLY)
                1519      $     ) )
613fa3996d Andr*1520 C****
                1521 C****
                1522 C**** -------------------------------------------------
                1523 C**** DEWFALL:
aab9c05604 Jean*1524 C****
                1525 C**** If dewfall adds to cir, insure that it does not fill
                1526 C**** beyond capacity.  If resulting throughfall adds to top soil layer,
3daafce20b Jean*1527 C**** insure that it also does not fill beyond capacity.
aab9c05604 Jean*1528 C****
613fa3996d Andr*1529       IF( CAPAC(ChNo) .GT. SATCAP(ChNo) ) THEN
                1530         THRU = CAPAC(ChNo) - SATCAP(ChNo)
                1531         CAPAC(ChNo) = SATCAP(ChNo)
aab9c05604 Jean*1532         SWET(ChNo,SFCLY) = SWET(ChNo,SFCLY) +
613fa3996d Andr*1533      &                         THRU / VGWMAX(SFCLY,ITYP(ChNo))
                1534         DEWRUN=0.
                1535         IF ( SWET(ChNo,SFCLY) .GT. 1. ) THEN
                1536           DEWRUN = ( SWET(ChNo,SFCLY) - 1. ) * VGWMAX(SFCLY,ITYP(ChNo))
                1537           SWET(ChNo,SFCLY) = 1.
                1538           RUNOFF(ChNo) = RUNOFF(ChNo) + DEWRUN/DTSTEP
                1539           RUNSRF(ChNo) = RUNSRF(ChNo) + DEWRUN/DTSTEP
                1540           ENDIF
aab9c05604 Jean*1541         FWSOIL(CHNO) = FWSOIL(CHNO) + (THRU-DEWRUN)/DTSTEP
613fa3996d Andr*1542         ENDIF
                1543 C****
                1544  100  CONTINUE
                1545 C****
                1546       RETURN
                1547       END
                1548 C****
                1549 C**** [ END WUPDAT ]
                1550 C****
                1551 C**** -----------------------------------------------------------------
                1552 C**** /////////////////////////////////////////////////////////////////
                1553 C**** -----------------------------------------------------------------
                1554 C****
                1555 C**** [ BEGIN GWATER ]
                1556 C****
                1557       SUBROUTINE GWATER (
                1558      I                     NCH,    ITYP,   WSMAX,  PHLAY, AKLAY, TC,
                1559      I                     DTSTEP, VGZDEX, VGSLOX, WETEQ1, WETEQ2,
                1560      I                     PHSAT,   AKSAT,
                1561      U                     SWET,   RUNOFF,  GDRAIN
                1562      &                    )
                1563 C****
                1564 C**** This subroutine computes diffusion between soil layers.
                1565 C****
                1566       IMPLICIT NONE
                1567 C****
4dad0f083a Andr*1568 #include "sibber.h"
613fa3996d Andr*1569 C****
                1570       INTEGER NCH
                1571       INTEGER ITYP(NCH), ChNo
a456aa407c Andr*1572       _RL  VGSLOX(NCH), RUNOFF(NCH), GDRAIN(NCH)
                1573       _RL  ZDEP12,   AKAVE,  GWFLUX,  ZDEP23,  HALFMX,  DHDZ,
613fa3996d Andr*1574      &    FAREA, TFM2, FRAMP
                1575 
a456aa407c Andr*1576       _RL   WSMAX(NLAY,NTYPS),   PHLAY(nch,NLAY),
613fa3996d Andr*1577      &       AKLAY(nch,NLAY),             TC(NCH),
                1578      &                  DTSTEP,    SWET(nch,NLAY),
                1579      &      VGZDEX(NLAY,nch),  WETEQ1(NCH),  WETEQ2(NCH),
                1580      &      PHSAT(NCH), AKSAT(NCH)
                1581 
                1582 C****
                1583 C**** ----------------------------------------------------------------
                1584 
                1585       DO 100 ChNo = 1, NCH
                1586 C****
                1587 C**** Diffusion between layer 1 and 2:
                1588       ZDEP12 = VGZDEX(SFCLY,ChNo) + VGZDEX(ROOTLY,ChNo)
                1589       IF( SWET(CHNO,SFCLY) .GT. WETEQ1(CHNO) ) THEN
                1590           FAREA=(SWET(CHNO,SFCLY) - WETEQ1(CHNO)) / (1. - WETEQ1(CHNO))
                1591           DHDZ = 1. + 2.*(PHSAT(CHNO)-PHLAY(CHNO,ROOTLY))/ZDEP12
                1592           AKAVE = AKSAT(CHNO)
                1593         ELSE
                1594           FAREA = 1.
                1595           DHDZ = 1. + 2.*(PHLAY(ChNo,SFCLY)-PHLAY(ChNo,ROOTLY))/ZDEP12
                1596           AKAVE = AKLAY(ChNo,ROOTLY)
                1597         ENDIF
                1598 C****
                1599       GWFLUX = 1000. * DTSTEP * AKAVE * DHDZ * FAREA
                1600 C****
                1601 C**** Test for limits on water holding capacity.
                1602 C****
                1603       HALFMX=0.5*ABS( SWET(CHNO,SFCLY)-WETEQ1(CHNO) )
                1604      &                               * WSMAX(SFCLY,ITYP(CHNO))
                1605       IF (GWFLUX .GE. 0.) THEN
aab9c05604 Jean*1606           GWFLUX = MIN( GWFLUX,
613fa3996d Andr*1607      &               SWET(ChNo,SFCLY) * WSMAX(SFCLY,ITYP(ChNo)) )
                1608           GWFLUX = MIN( GWFLUX,    WSMAX(ROOTLY,ITYP(ChNo)) -
                1609      &              SWET(ChNo,ROOTLY) * WSMAX(ROOTLY,ITYP(ChNo)) )
                1610           GWFLUX = MIN( GWFLUX, HALFMX )
                1611         ELSE
                1612           GWFLUX = -MIN( ABS(GWFLUX),
                1613      &                SWET(ChNo,ROOTLY) * WSMAX(ROOTLY,ITYP(ChNo)) )
                1614           GWFLUX = -MIN( ABS(GWFLUX),   WSMAX(SFCLY,ITYP(ChNo)) -
                1615      &                SWET(ChNo,SFCLY) * WSMAX(SFCLY,ITYP(ChNo)) )
                1616           GWFLUX = -MIN( ABS(GWFLUX),  HALFMX )
                1617         ENDIF
                1618 C****
                1619 C**** Prevent diffusion when ground is frozen (INCLUDE RAMPING):
                1620       TFM2=TF-2.
                1621       FRAMP=(TC(CHNO)-TFM2)/2.
32362b8fa8 Cons*1622       FRAMP=MIN(1. _d 0, MAX(0. _d 0,FRAMP) )
613fa3996d Andr*1623       GWFLUX=GWFLUX*FRAMP
                1624 C****
                1625 C**** Update water contents
                1626       SWET(ChNo,SFCLY) = SWET(ChNo,SFCLY) -
                1627      &                        GWFLUX / WSMAX(SFCLY,ITYP(ChNo))
                1628       SWET(ChNo,ROOTLY) = SWET(ChNo,ROOTLY) +
                1629      &                        GWFLUX / WSMAX(ROOTLY,ITYP(ChNo))
                1630 C****
aab9c05604 Jean*1631 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
613fa3996d Andr*1632 C**** Diffusion between root and recharge layers:
                1633       ZDEP23 = VGZDEX(ROOTLY,ChNo) + VGZDEX(RECHLY,ChNo)
                1634       DHDZ=1. + 2. * (PHLAY(ChNo,ROOTLY) - PHLAY(ChNo,RECHLY)) / ZDEP23
                1635       IF(DHDZ.GE.0.) THEN
                1636           AKAVE = AKLAY(ChNo,ROOTLY)
                1637         ELSE
                1638           AKAVE = AKLAY(ChNo,RECHLY)
                1639         ENDIF
                1640       GWFLUX = 1000. * DTSTEP * AKAVE * DHDZ
                1641 C****
                1642 C**** Test for limits on water holding capacity:
                1643       HALFMX=0.5*ABS( SWET(CHNO,ROOTLY)-WETEQ2(CHNO) )
                1644      &                               * WSMAX(ROOTLY,ITYP(CHNO))
                1645       IF (GWFLUX .GE. 0.) THEN
                1646           GWFLUX = MIN( GWFLUX,
                1647      &                    SWET(ChNo,ROOTLY) * WSMAX(ROOTLY,ITYP(ChNo)) )
                1648           GWFLUX = MIN( GWFLUX,    WSMAX(RECHLY,ITYP(ChNo)) -
                1649      &                    SWET(ChNo,RECHLY) * WSMAX(RECHLY,ITYP(ChNo)) )
                1650           GWFLUX = MIN( GWFLUX, HALFMX )
                1651         ELSE
                1652           GWFLUX = -MIN( ABS(GWFLUX),
                1653      &                    SWET(ChNo,RECHLY) * WSMAX(RECHLY,ITYP(ChNo)) )
                1654           GWFLUX = -MIN( ABS(GWFLUX), WSMAX(ROOTLY,ITYP(ChNo)) -
                1655      &                    SWET(ChNo,ROOTLY) * WSMAX(ROOTLY,ITYP(ChNo)) )
                1656           GWFLUX = -MIN( ABS(GWFLUX), HALFMX )
                1657         ENDIF
                1658 C****
                1659 C**** Prevent diffusion when ground is frozen (INCLUDE RAMPING):
                1660       FRAMP=(TC(CHNO)-TFM2)/2.
32362b8fa8 Cons*1661       FRAMP=MIN(1. _d 0, MAX(0. _d 0,FRAMP) )
613fa3996d Andr*1662       GWFLUX=GWFLUX*FRAMP
                1663 C****
                1664 C**** Update water contents
                1665       SWET(ChNo,ROOTLY) = SWET(ChNo,ROOTLY) -
                1666      &                            GWFLUX / WSMAX(ROOTLY,ITYP(ChNo))
                1667       SWET(ChNo,RECHLY) = SWET(ChNo,RECHLY) +
                1668      &                            GWFLUX / WSMAX(RECHLY,ITYP(ChNo))
                1669       GDRAIN(CHNO)=GWFLUX/DTSTEP
                1670 C****
                1671 C**** - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                1672 C**** Flux of moisture out of lower soil layer to water table
                1673 C**** (approximation to SiB)
                1674 C****
                1675       GWFLUX = VGSLOX(ChNo) * AKLAY(ChNo,RECHLY) * 1000. * DTSTEP
                1676 
                1677 C**** Prevent diffusion when ground is frozen (INCLUDE RAMPING):
                1678       FRAMP=(TC(CHNO)-TFM2)/2.
32362b8fa8 Cons*1679       FRAMP=MIN(1. _d 0, MAX(0. _d 0,FRAMP) )
613fa3996d Andr*1680       GWFLUX=GWFLUX*FRAMP
                1681 
                1682       GWFLUX = MIN( GWFLUX,
                1683      &                 SWET(ChNo,RECHLY) * WSMAX(RECHLY,ITYP(ChNo)) )
                1684       SWET(ChNo,RECHLY) = SWET(ChNo,RECHLY) -
                1685      &                                GWFLUX / WSMAX(RECHLY,ITYP(ChNo))
                1686 C****
                1687       RUNOFF (ChNo) = RUNOFF (ChNo) + GWFLUX/DTSTEP
                1688 C****
                1689  100  CONTINUE
                1690 C****
                1691       RETURN
                1692       END
                1693 C****
                1694 C**** [ END GWATER ]
                1695
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!