1. C**** 2. C**** C073H.S Fortran source code for ground Hydrology 2003/11/19 3. C**** 4. C**** NASA/GISS Climate Model III Frank Abramopoulos 5. C**** 6. C**** NASA/Goddard Space Flight Center 7. C**** Institute for Space Studies 8. C**** 2880 Broadway, New York, NY 10025 9. C**** U.S.A. 10. C**** 10.3 C**** C073H: WEARTH is calculated correctly in GHIJ 10.4 C**** H041AM9: correct line 2028 in subroutine FLHG 10.5 C**** H041M9: when snow masks vegetation, it receives radiation directly 10.6 C**** H040M9: snow evaporation corrected in subroutine QSBAL 10.7 C**** H038M9: snow melts before ice in layer 1 in SNMLT. 10.9 C**** H036M9: VDATA now contians 10 types: crops and second bare ground 11. C**** H035M9: soils57 but adds GHIJ interface routine for Gary's GCM 12. C**** SOILS45 10/4/93 13. C**** USES BEDROCK AS A SOIL TEXTURE. SOIL DEPTH OF 3.5M 14. C**** EVERYWHERE, WHERE LAYERS CAN HAVE BEDROCK. 15. C**** REQUIRES SM693.DATA INSTEAD OF SM691.DATA. 16. C**** SDATA NEEDS TO BE CHANGED IN CALLING PROGRAM. 17. C**** SOILS44B 8/25/93 18. C**** USES SNOW CONDUCTIVITY OF .088 W M-1 C-1 INSTEAD OF .3 19. C**** SOILS44 8/16/93 20. C**** ADDS BEDROCK FOR HEAT CALCULATIONS, TO FILL OUT THE 21. C**** NUMBER OF LAYERS TO NGM. 22. C**** SOILS43 6/25/93 23. C**** COMMENTS OUT CALL TO FHLMT HEAT FLUX LIMITS. 24. C**** USES GHINIJ TO RETURN WFC1, ELIMINATES REWFC. 25. C**** SOILS42 6/15/93 26. C**** ADDS SNOW INSULATION 27. C**** SOILS41 5/24/93 28. C**** USES SNOW MASKING DEPTH FROM VEGETATION HEIGHT TO DETERMINE 29. C**** FRACTION OF SNOW THAT IS EXPOSED. 30. C**** THE @PROCESS DPC DIRECTIVES ARE USED INSTEAD OF -qdpc=e 31. C**** FOR DOUBLE PRECISION RUNS. THEY SHOULD BE CHANGED 32. C**** TO THE @PROCESS NODPC DIRECTIVES (OR COMMENTED OUT) 33. C**** FOR SINGLE PRECISION RUNS. 34. C**** RETH MUST BE CALLED PRIOR TO RETP. 35. C**** SOILS40 5/10/93 36. C**** REMOVES SNOW FROM CANOPY AND PLACES IT ON VEGETATED SOIL. 37. C**** SOILS39 4/19/93 38. C**** MODIFICATIONS FOR REAL*8 OR REAL*4 RUNS. COMMON BLOCK 39. C**** ORDERING CHANGED FOR EFFICIENT ALIGNMENT. SDATA,FDATA, 40. C**** AND VDATA ARE EXPLICITLY REAL*4. ON IBM RS/6000, SHOULD 41. C**** BE COMPILED WITH -qdpc=e OPTION FOR REAL*8 OPERATION. 42. C**** TO RUN REAL*4, CHANGE IMPLICIT STATEMENT IN INCLUDE FILE. 43. C**** SOILS38 2/9/93 44. C**** ADDS HEAT FLUX CORRECTION TO HANDLE VARYING COEFFICIENTS 45. C**** OF DRAG. 46. C**** SOILS37 1/25/93 47. C**** CHANGES SOIL CRUSTING PARAMETER KU/D FROM .05 PER HOUR TO .1, 48. C**** TO AGREE WITH MORIN ET AL. 49. C**** SOILS36 11/12/92 50. C**** CALCULATES HEAT CONDUCTIVITY OF SOILS USING DEVRIES METHOD. 51. C**** CHANGES LOAM MATERIAL HEAT CAPACITY AND CONDUCTIVITY 52. C**** TO MINERAL VALUES. 53. C**** SOILS35 10/27/92 54. C**** INCLUDES EFFECT OF SOIL CRUSTING FOR INFILTRATION BY 55. C**** MODIFYING HYDRAULIC CONDUCTIVITY CALCULATION OF LAYER 56. C**** 1 IN HYDRA. 57. C**** SOILS34 8/28/92 58. C**** USES EFFECTIVE LEAF AREA INDEX ALAIE FOR PURPOSES OF 59. C**** CANOPY CONDUCTANCE CALCULATION. 60. C**** SOILS33 8/9/92 61. C**** CHANGES CANOPY WATER STORAGE CAPACITY TO .1MM PER LAI FROM 1. 62. C**** SOILS32 7/15/92 63. C**** 1) FOR PURPOSES OF INFILTRATION ONLY, REDUCES SOIL CONDUCTIVITY 64. C**** BY (1-THETR*FICE) INSTEAD OF (1-FICE). 65. C**** 2) BETAD IS REDUCED BY FRACTION OF ICE IN EACH LAYER. 66. C**** 3) TRANSPIRED WATER IS REMOVED BY BETAD FRACTION IN EACH LAYER, 67. C**** INSTEAD OF BY FRACTION OF ROOTS. PREVENTS NEGATIVE RUNOFF. 68. C**** 4) SPEEDS UP HYDRA BY USING DO LOOP INSTEAD OF IF CHECK, 69. C**** BY USING INTERPOLATION POINT FROM BISECTION INSTEAD OF LOGS, 70. C**** AND BY AVOIDING UNNECESSARY CALLS TO HYDRA. ALSO ELIMATES CALL 71. C**** TO HYDRA IN MA89EZM9.F. 72. C**** SOILS31 7/1/92 73. C**** 1) FIXES FRACTION OF ROOTS WHEN SOIL DEPTH IS LESS THAN ROOT 74. C**** DEPTH, THUS FIXING NON-CONSERVATION OF WATER. 75. C**** SOILS30 6/4/92 76. C**** 1) USES ACTUAL FINAL SNOW DEPTH IN FLUX LIMIT CALCULATIONS, 77. C**** INSTEAD OF UPPER AND LOWER LIMITS. FIXES SPURIOUS DRYING 78. C**** OF FIRST LAYER. 79. C**** 101. SUBROUTINE RETH 102. C**** REVISES VALUES OF THETA BASED UPON W. 103. C**** INPUT: 104. C**** W - WATER DEPTH, M 105. C**** WS - SATURATED WATER DEPTH, M 106. C**** DZ - LAYER THICKNESS, M 107. C**** THETS - SATURATED THETA 108. C**** SNOWD - SNOW DEPTH, WATER EQUIVALENT M 109. C**** SNOWM - SNOW MASKING DEPTH, WATER EQUIVALENT M 110. C**** OUTPUT: 111. C**** THETA - WATER SATURATION 112. C**** FW - FRACTION OF WET CANOPY 113. C**** FD - FRACTION OF DRY CANOPY 114. C**** FM - FRACTION OF SNOW THAT IS EXPOSED, OR MASKING. 115. INCLUDE 'SOILS035.COM' 116. C**** SOILS28 Common block 9/25/90 117. ONE=1. 118. DO 10 IBV=1,2 119. DO 10 L=1,N 120. THETA(L,IBV)=W(L,IBV)/DZ(L) 121. 10 CONTINUE 122. C**** DO CANOPY LAYER 123. C**** HERE THETA IS THE FRACTION OF CANOPY COVERED BY WATER 124. IF(WS(0,2).GT.0.)THEN 125. THETA(0,2)=(W(0,2)/WS(0,2))**(2./3.) 126. ELSE 127. THETA(0,2)=0. 128. ENDIF 129. THETA(0,2)=MIN(THETA(0,2),ONE) 130. C**** CORRECT FIRST SOIL LAYER FOR SNOW 131. DO 20 IBV=1,2 132. THETA(1,IBV)=THETA(1,IBV)-SNOWD(IBV)/DZ(1) 133. 20 CONTINUE 134. C**** FRACTION OF WET CANOPY FW 135. FW=THETA(0,2) 136. C**** DETERMINE FM FROM SNOWD DEPTH AND MASKING DEPTH 137. FM=1.-EXP(-SNOWD(2)/(SNOWM+1.E-12)) 138. C**** CORRECT FRACTION OF WET CANOPY BY SNOW FRACTION 139. FW=FW+FM*(1.-FW) 140. FD=1.-FW 141. RETURN 142. END 201. SUBROUTINE HYDRA 202. C ROUTINE TO RETURN THE EQULIBRIUM VALUE OF H IN A MIXED SOIL 203. C LAYER. THE H IS SUCH THAT EACH SOIL TEXTURE HAS THE SAME 204. C VALUE OF H, BUT DIFFERING VALUES OF THETA. 205. C HYDRA ALSO CALCULATES THE CONDUCTIVITY XK AND DIFFUSSIVITY D. 206. C**** INPUT: 207. C**** THETA(L,IBV) - VOLUMETRIC WATER CONCENTRATION 208. C**** THETM(L,IBV) - MINIMUM THETA 209. C**** THETS(L,IBV) - MAXIMUM THETA 210. C**** NTH - NUMBER OF H0 INTERVALS IN TABLE, A POWER OF TWO. 211. C**** HLM(J) - TABLE OF H VALUES, FROM 0 (AT J=0) TO HMIN (AT J=NTH) 212. C**** THM(J,I) - VALUE OF RELATIVE THETA AT HLM(J) IN TEXTURE I, 213. C**** RANGING BETWEEN THETS(L,IBV) AT J=0 TO THETM(L,IBV) AT J=NTH. 214. C**** OUTPUT: 215. C**** H - POTENTIAL, M, INCLUDING BOTH MATRIC AND GRAVITATIONAL 216. C**** D - DIFFUSIVITY, DL. 217. C**** XK - CONDUCTIVITY M/S. 218. INCLUDE 'SOILS035.COM' 219. C**** SOILS28 Common block 9/25/90 220. C SOLVE FOR H USING BISECTION 221. C WE ASSUME THAT IF J1.LT.J2 THEN HLM(J1).GT.HLM(J2) 222. C AND THM(J1,I).GT.THM(J2,I). 223. C 224. C ALGDEL=ALOG(1.+ALPH0) 225. C 226. ZERO=0. 227. XKUD=2.78E-5 228. JCM=NINT(LOG(FLOAT(NTH))/LOG(2.)) 229. DO 600 IBV=1,2 230. XK(N+1,IBV)=0.0 231. XKU(0,IBV)=0. 232. DO 550 L=1,N 233. J1=0 234. J2=NTH 235. THR1=THETS(L,IBV) 236. THR2=THETM(L,IBV) 237. THR0=THETA(L,IBV) 238. THR0=MIN(THR1,THR0) 239. THR0=MAX(THR2,THR0) 240. DO 400 JC=1,JCM 241. J=(J1+J2)/2 242. THR=0. 243. DO 20 I=1,IMT-1 244. THR=THR+THM(J,I)*Q(I,L) 245. 20 CONTINUE 246. IF(THR-THR0) 100,300,200 247. 100 CONTINUE 248. C HERE THR IS TOO SMALL, BISECT ON LOW J END 249. J2=J 250. THR2=THR 251. GO TO 400 252. 200 CONTINUE 253. C HERE THR IS TOO LARGE, BISECT ON HIGH J END 254. J1=J 255. THR1=THR 256. GO TO 400 257. 300 CONTINUE 258. C HERE THR IS EQUAL TO THR0 259. HL=HLM(J) 260. J1=J 261. THR1=THR0 262. C THE STRANGE VALUE FOR THR2 BELOW IS ONLY FOR CALCULATING TEMP 263. THR2=-10. 264. GO TO 500 265. 400 CONTINUE 266. C HERE THETA IS BETWEEN TWO ADJACENT THR'S. INTERPOLATE. 267. HL=(HLM(J1)*(THR0-THR2)+HLM(J2)*(THR1-THR0))/(THR1-THR2) 268. 500 CONTINUE 269. C**** ONLY FILLING HL ARRAY WITH MATRIC POTENTIAL (GRAVITATIONAL TO BE 270. C**** ADDED LATER) 271. H(L,IBV)=HL 272. HZ=HL 273. C**** CALCULATE DIFFUSIVITY 274. ITH=J1 275. TEMP=(THR1-THR0)/(THR1-THR2) 276. D1=0. 277. D2=0. 278. XKU1=0. 279. XKU2=0. 280. DO 541 I=1,IMT-1 281. D1=D1+Q(I,L)*DLM(ITH,I) 282. D2=D2+Q(I,L)*DLM(ITH+1,I) 283. XKU1=XKU1+Q(I,L)*XKLM(ITH,I) 284. XKU2=XKU2+Q(I,L)*XKLM(ITH+1,I) 285. 541 CONTINUE 286. DL=(1.-TEMP)*D1+TEMP*D2 287. DL=(1.-FICE(L,IBV))*DL 288. D(L,IBV)=DL 289. C**** CALCULATE CONDUCTIVITY 290. XKLU=(1.-TEMP)*XKU1+TEMP*XKU2 291. XKLU=(1.-FICE(L,IBV))*XKLU 292. XKU(L,IBV)=XKLU 293. IF(L.EQ.1) THEN 294. XK1=0. 295. DO 641 I=1,IMT-1 296. XK1=XK1+QK(I,1)*XKLM(0,I) 297. 641 CONTINUE 298. XKL=XK1 299. XKL=XKL/(1.+XKL/(-ZC(1)*XKUD)) 300. XKL=(1.-FICE(1,IBV)*THETA(1,IBV)/THETS(1,IBV))*XKL 301. XKL=MAX(ZERO,XKL) 302. XK(1,IBV)=XKL 303. ELSE 304. XK(L,IBV)=SQRT(XKU(L-1,IBV)*XKU(L,IBV)) 305. END IF 306. 550 CONTINUE 307. 600 CONTINUE 308. C ADD GRAVITATIONAL POTENTIAL TO HL 309. DO 700 L=1,N 310. DO 680 IBV=1,2 311. H(L,IBV)=H(L,IBV)+ZC(L) 312. 680 CONTINUE 313. 700 CONTINUE 314. RETURN 315. END 401. SUBROUTINE HL0 402. C**** HL0 SETS UP A TABLE OF THETA VALUES AS A FUNCTION OF MATRIC 403. C**** POTENTIAL, H. H IS TABULATED IN A GEOMETRIC SERIES FROM 404. C**** 0 TO HMIN, WITH A FIRST STEP OF DELH1. THE THETA VALUES 405. C**** DEPEND NOT ONLY ON THE MATRIC POTENTIAL, BUT ALSO ON THE 406. C**** SOIL TEXTURE. WE SOLVE A CUBIC EQUATION TO DETERMINE 407. C**** THETA AS A FUNCTION OF H. HL0 ALSO OUTPUTS THE CONDUCTIVITY 408. C**** AND DIFFUSIVITY TABLES. 409. C**** INPUT: 410. C**** A - MATRIC POTENTIAL FUNCTION PARAMETERS 411. C**** B - HYDRAULIC CONDUCTIVITY FUNCTION PARAMETERS 412. C**** P - HYDRAULIC DIFFUSIVITY FUNCTION PARAMETERS 413. C**** SAT - SATURATED THETAS 414. C**** OUTPUT: 415. C**** THM(J,I) - THETA AT J'TH H POINT FOR TEXTURE I 416. C**** XKLM(J,I) - CONDUCTIVITY AT J'TH H POINT FOR TEXTURE I 417. C**** DLM(J,I) - DIFFUSIVITY AT J'TH H POINT FOR TEXTURE I 418. INCLUDE 'SOILS035.COM' 419. C**** SOILS28 Common block 9/25/90 420. PARAMETER (NEXP=6) 421. PARAMETER (C=2.3025851) 422. DIMENSION ARG(IMT-1) 423. DIMENSION SAT(IMT-1) 424. DATA SAT/.394,.537,.577,.885/ 425. SXTN=16. 426. NTH=2**NEXP 427. HLM(0)=0.0 428. DELH1=-0.00625 429. HMIN=-1000. 430. DELHN=DELH1 431. C SOLVE FOR ALPH0 IN S=((1+ALPH0)**N-1)/ALPH0 432. S=HMIN/DELH1 433. ALPH0=1./8. 434. 10 ALPH0O=ALPH0 435. ALPH0=(S*ALPH0+1.)**(1./NTH)-1. 436. IF(ABS(ALPH0O-ALPH0).GE.1.E-8) GO TO 10 437. ALPLS1=1.0+ALPH0 438. ALGDEL=LOG(1.+ALPH0) 439. DO 100 J=1,NTH 440. HLM(J)=HLM(J-1)+DELHN 441. DELHN=ALPLS1*DELHN 442. 100 CONTINUE 443. MMAX=100 444. XTOL=1.E-6 445. DO 200 I=1,IMT-1 446. THM(0,I)=1.00 447. DO 150 J=1,NTH 448. HS=-EXP(C*(A(1,I)+A(2,I)+A(3,I)+A(4,I))) 449. A1=A(3,I)/A(4,I) 450. A2=(A(2,I)-(LOG(-HLM(J)-HS))/C)/A(4,I) 451. A3=A(1,I)/A(4,I) 452. TESTH=THM(J-1,I) 453. DO 130 M=1,MMAX 454. FUNC=(TESTH**3)+(A1*(TESTH**2))+(A2*(TESTH))+A3 455. DFUNC=(3*TESTH**2)+(2*A1*TESTH)+A2 456. DIFF=FUNC/DFUNC 457. TESTH=TESTH-DIFF 458. IF(ABS(DIFF).LT.XTOL) GO TO 140 459. 130 CONTINUE 460. PRINT *,'MAX # ITERATIONS:',MMAX 461. 140 THM(J,I)=TESTH 462. 150 CONTINUE 463. 200 CONTINUE 464. DO 280 J=0,NTH 465. DO 245 I=1,IMT-1 466. XKLM(J,I)=0. 467. ARG(I)=0. 468. DO 240 L=-1,2 469. ARG(I)=ARG(I)+B(L+2,I)*THM(J,I)**L 470. 240 CONTINUE 471. ARG(I)=MIN(ARG(I),SXTN) 472. ARG(I)=MAX(ARG(I),-SXTN) 473. XKLM(J,I)=EXP(C*ARG(I)) 474. 245 CONTINUE 475. DO 265 I=1,IMT-1 476. DLM(J,I)=0. 477. ARG(I)=0. 478. DO 260 L=-1,2 479. ARG(I)=ARG(I)+P(L+2,I)*THM(J,I)**L 480. 260 CONTINUE 481. ARG(I)=MIN(ARG(I),SXTN) 482. ARG(I)=MAX(ARG(I),-SXTN) 483. DLM(J,I)=EXP(C*ARG(I)) 484. 265 CONTINUE 485. 280 CONTINUE 486. DO 350 J=0,NTH 487. DO 310 K=1,IMT-1 488. THM(J,K)=THM(J,K)*SAT(K) 489. 310 CONTINUE 490. 350 CONTINUE 491. RETURN 492. END 601. SUBROUTINE FL 602. C**** EVALUATES THE FLUX BETWEEN LAYERS. 603. C**** INPUT: 604. C**** H - SOIL POTENTIAL OF LAYERS, M 605. C**** XK - CONDUCTIVITY OF LAYERS, M S-1 606. C**** ZC - LAYER CENTERS, M 607. C**** OUTPUT: 608. C**** F - FLUXES BETWEEN LAYERS, M S-1 609. C**** XINFC - INFILTRATION CAPACITY, M S-1 610. INCLUDE 'SOILS035.COM' 611. C**** SOILS28 Common block 9/25/90 612. C**** 613. DO 15 IBV=1,2 614. 15 F(N+1,IBV)=0. 615. C**** 616. DO 20 IBV=1,2 617. DO 20 L=2,N 618. F(L,IBV)=-XK(L,IBV)*(H(L-1,IBV)-H(L,IBV))/(ZC(L-1)-ZC(L)) 619. 20 CONTINUE 620. C**** PUT INFILTRATION MAXIMUM INTO XINFC 621. DO 25 IBV=1,2 622. XINFC(IBV)=XK(1,IBV)*H(1,IBV)/ZC(1) 623. 25 CONTINUE 624. RETURN 625. END 701. SUBROUTINE QSBAL 702. C**** FINDS QS THAT BALANCES FLUXES. 703. C**** OBTAINS QS BY SUCCESSIVE APPROXIMATION. 704. C**** CALCULATES EVAPORATION. 705. C**** INPUT: 706. C**** CH - HEAT CONDUCTIVITY COEFFICIENT FROM GROUND TO SURFACE 707. C**** VSM - SURFACE LAYER WIND SPEED, M S-1 708. C**** RHO - AIR DENSITY, KG M-3 709. C**** EDDY - TRANSFER COEFFICIENT FROM SURFACE TO FIRST ATMOSPHERE 710. C**** THETA - WATER SATURATION OF LAYERS AND CANOPY 711. C**** TP - TEMPERATURES OF LAYERS AND CANOPY, C 712. C**** PRES - ATMOSPHERIC PRESSURE AT GROUND 713. C**** Z1 - HEIGHT OF FIRST LAYER, M 714. C**** ZS - HEIGHT OF SURFACE LAYER, M 715. C**** DZ - LAYER THICKNESSES, M 716. C**** SNOWD - SNOW DEPTHS, EQUIVALENT WATER M 717. C**** PR - PRECIPITATION, M S-1 718. C**** Q1 - MIXING RATIO OF FIRST LAYER 719. C**** FR - FRACTION OF ROOTS IN LAYER 720. C**** FB - FRACTION OF BARE SOIL 721. C**** FV - FRACTION OF VEGETATED SOIL 722. C**** HW - WILTING POINT, M 723. C**** OUTPUT: 724. C**** QS - MIXING RATIO AT SURFACE LAYER 725. C**** EVAP - EVAPORATION FROM BARE AND VEGETATED REGIONS, M S-1 726. C**** EVAPW - EVAPORATION FROM WET CANOPY, M S-1, INCLUDING FROM SNOW 727. C**** EVAPD - EVAPORATION FROM DRY CANOPY, M S-1 728. C**** EVAPS - EVAPORATION FROM SNOW FROM CANOPY, M S-1 729. C**** BETAD - DRY CANOPY BETA, BASED ON ROOTS 730. C**** 731. C**** USES: COND 732. INCLUDE 'SOILS035.COM' 733. C**** SOILS28 Common block 9/25/90 734. C**** MIXING RATIO FOR GIVEN TEMPERATURE TX(C) AND PRESSURE PX(MB) 735. QSAT(TX,PX)=.622 * 33.8639*( (0.00738*TX+0.8072)**8 - 736. * 0.000019*ABS(1.8*TX+48)+0.001316 )/PX 737. ZERO=0. 738. C**** QL HAS MASS OF WATER VAPOR IN FIRST ATMOSPHERE LAYER, KG M-2 738.1 QM1 = 0 738.2 QL=QM1 739. IF(IGCM.EQ.-1) QL=1.E+7 740. QLDT=.001*QL/DT 741. C CNA IS THE CONDUCTANCE OF THE ATMOSPHERE 742. CNA=CH*VSM 743. RHO3=.001*RHO 744. IF(IGCM.GE.0 .AND. IGCM.LE.3) XL=EDDY/(Z1-ZS) 745. C CALCULATE BARE SOIL AND CANOPY MIXING RATIOS 746. QB = QSAT(TP(1,1),PRES) 747. QC = QSAT(TCS,PRES) 748. C ON FIRST ITERATION, ASSUME BETA'S = 1 749. BETAB=1. 750. BETAV=1. 751. IF(IGCM.GE.0 .AND. IGCM.LE.3) 752. & QS=(FB*BETAB*CNA*QB+FV*BETAV*CNA*QC+XL*Q1) 753. & /(FB*BETAB*CNA+FV*BETAV*CNA+XL+1.E-12) 754. EPS=5.E-5 755. C BARE SOILS DIFFUSIVITY DD 756. DD=D(1,1) 757. C BETAD IS THE THE ROOT BETA FOR TRANSPIRATION. 758. C HW IS THE WILTING POINT. 759. C FR(L) IS THE FRACTION OF ROOTS IN LAYER L 760. BETAD=0. 761. DO 30 L=1,N 762. BETADL(L)=(1.-FICE(L,2))*FR(L)*MAX((HW-H(L,2))/HW,ZERO) 763. BETAD=BETAD+BETADL(L) 764. 30 CONTINUE 765. ABETAD=BETAD 766. C CANOPY CONDUCTIVITY CNC 767. CALL COND 768. C SURFACE LAYER MIXING RATIO TO BALANCE FLUXES 769. ITR=1 770. 10 QSO=QS 771. C POTENTIAL EVAPORATION FOR BARE SOIL AND CANOPY 772. EPB=RHO3*CNA*(QB-QS) 773. EPC=RHO3*CNA*(QC-QS) 774. C BARE SOIL CORRECTION 775. C DIFFUSION LIMITED FLUX ED 776. ED=2.467*DD*(THETA(1,1)-THETM(1,1))/DZ(1) 777. C CORRECT IF NOT POTENTIAL EVAPORATION 778. C EVAP(1) IS EVAPORATION FROM BARE SOIL 779. IF(SNOWD(1).GT.0.)THEN 780. EVAP(1)=EPB 781. ELSE 782. EVAP(1)=MIN(EPB,ED+PR) 783. ENDIF 784. C VEGETATED SOIL CORRECTION 785. C EVAP(2) IS EVAPORATION FROM VEGETATED LAND 786. C EVAPD IS DRY EVAPORATION (TRANSPIRATION) FROM CANOPY 787. C EVAPW IS WET EVAPORATION FROM CANOPY (FROM INTERCEPTION) 788. IF(EPC.GT.0) THEN 789. EVAPW=EPC*FW 790. C**** LIMIT THE WET CANOPY EVAPORATION TO CANOPY WATER 791. EVAPW=MIN(EVAPW,W(0,2)/DT+PR) 792. EVAPD=EPC*FD 793. C**** IF POSITIVE EVAPORATION, LIMIT DRY CANOPY EVAPORATION TO TRANS 794. BETAT=CNC/(CNC+CNA+1.E-12) 795. EVAPD=MIN(EVAPD,EVAPD*BETAT) 796. ELSE 797. EVAPW=MAX(EPC,-QLDT) 798. EVAPD=0. 799. END IF 800. C**** EVAPORATION FROM VEGETATED SNOW REGION IS FROM THAT PART 801. C**** OF THE WET CANOPY THAT REPRESENTS SNOW 802. EVAPS = EPC*FM 803. C**** RESTRICT CONDENSATION TO WATER AVAILABLE IN FIRST ATMOSPHERE 804. EVAP(1)=MAX(EVAP(1),-QLDT) 805. EVAP(2)=EVAPW+EVAPD 806. C**** CALCULATE BETAS AND Q OF SURFACE LAYER 807. IF(EPB.LE.0.) THEN 808. BETAB=1.0 809. ELSE 810. BETAB=EVAP(1)/EPB 811. END IF 812. IF(EPC.LE.0.) THEN 813. BETAV=1.0 814. ELSE 815. BETAV=EVAP(2)/EPC 816. END IF 817. C**** FOR OVERALL BETA, USE WEIGHTED AVERAGE OF BETAB AND BETAV. 818. C**** DON'T USE TOTAL EVAP OVER TOTAL POTENTIAL EVAP. THIS AVOIDS 819. C**** THE POSSIBILITY OF NEGATIVE BETA. 820. BETA=FB*BETAB+FV*BETAV 821. ABETAV=BETAV 822. ABETAT=BETAT 823. ABETAB=BETAB 824. ABETA=BETA 825. ACNA=CNA 826. ACNC=CNC 827. IF(IGCM.GE.0 .AND. IGCM.LE.3) 828. & QS=(FB*BETAB*CNA*QB+FV*BETAV*CNA*QC+XL*Q1) 829. & /(FB*BETAB*CNA+FV*BETAV*CNA+XL+1.E-12) 830. 70 CONTINUE 831. C LOOP BACK UNTIL QS CONVERGED 832. IF(ITR.GE.60)THEN 833. WRITE(*,*)'QSBAL:1',ID,ITR,QS,QSO 834. WRITE(*,*)'QSBAL:2',FB,BETAB,CNA,QB 835. WRITE(*,*)'QSBAL:3',FV,BETAV,QC,XL 836. WRITE(*,*)'QSBAL:4',EVAP(1),EVAP(2),EPB,EPC 837. WRITE(*,*)'QSBAL:5',EVAPW,EVAPD,ED,PR 838. WRITE(*,*)'QSBAL:6',W(0,2),QLDT,DT,CNC 839. WRITE(*,*)'QSBAL:7',BETAD,Q1,ALAI,RS 840. WRITE(*,*)'QSBAL:8',SRHT,TP(1,1),TCS,TS 841. ENDIF 842. IF(ITR.GE.64)THEN 843. CALL OUTW(0) 844. STOP 'QSBAL' 845. ENDIF 846. ITR=ITR+1 847. IF(ABS(QSO-QS).GT.EPS)GO TO 10 848. DO 100 IBV=1,2 849. L=2-IBV 850. SNSH(IBV)=SHA*RHO*CNA*(TP(L,IBV)-TS+TFRZ) 851. XLTH(IBV)=EVAP(IBV)*ELH 852. 100 CONTINUE 853. RETURN 854. END 1001. SUBROUTINE FLG 1002. C**** CALCULATES THE GROUND WATER FLUXES (TO THE SURFACE) 1003. C**** INPUT: 1004. C**** EVAP - EVAPORATION FROM BARE AND VEGETATED REGIONS, M S-1 1005. C**** EVAPW - EVAPORATION FROM WET CANOPY, M S-1 1006. C**** PR - PRECIPITATION, M S-1 1007. C**** HTPR - HEAT OF PRECIPITATION 1008. C**** FSN - HEAT OF FUSION 1009. C**** PRE - EXTRA PRECIPITATION, I.E. SMOWMELT, M S-1 1010. C**** PRFR - FRACTION BY AREA OF PRECIPITATION 1011. C**** OUTPUT: 1012. C**** F - WATER FLUXES FROM GROUND AND CANOPY 1013. C**** SNOWF - SNOW FALL, EQUIVALENT WATER M S-1 1014. C**** DR - CANOPY DRIP, M S-1 1015. INCLUDE 'SOILS035.COM' 1016. C**** SOILS28 Common block 9/25/90 1017. ZERO=0. 1018. C CALCULATE SNOW FALL. SNOWF IS SNOW FALL, M S-1 OF WATER DEPTH. 1019. SNOWF=0. 1020. IF(HTPR.LT.0.)SNOWF=MIN(-HTPR/FSN,PR) 1021. C SNOWFS IS THE LARGE SCALE SNOW FALL. 1022. SNOWFS=0. 1023. IF(HTPRS.LT.0.)SNOWFS=MIN(-HTPRS/FSN,PRS) 1024. C BARE SOIL 1025. C UPWARD FLUX FROM FIRST LAYER IS EVAPORATION LESS PRECIPITATION 1026. F(1,1)=-PR+EVAP(1) 1027. C UPWARD FLUX FROM WET CANOPY, INCLUDING EVAPORATION FROM SNOW. 1028. F(0,2)=-PR+EVAPW 1029. PTMPS=PRS-SNOWFS 1030. PTMPS=PTMPS-EVAPW 1031. PTMP=PR-PRS-(SNOWF-SNOWFS) 1032. C USE EFFECTS OF SUBGRID SCALE PRECIPITATION TO CALCULATE DRIP 1033. PM=1.E-6 1034. PMAX=FD*PM 1035. DRS=MAX(PTMPS-PMAX,ZERO) 1036. DR=DRS 1037. IF(PTMP.GT.0.)THEN 1038. PFAC=(PMAX-PTMPS)*PRFR/PTMP 1039. IF(PFAC.GE.0.)THEN 1040. IF(PFAC.LT.30.000) DR=PTMP*EXP(-PFAC) 1041. ELSE 1042. DR=PTMP+PTMPS-PMAX 1043. ENDIF 1044. ENDIF 1045. C VEGETATED SOIL 1046. C UPWARD FLUX FROM SOIL SURFACE IS MINUS DRIP LESS SNOWFALL 1047. C PLUS THE EVAPORATION FROM SNOW 1048. F(1,2)=-DR-SNOWF+EVAPS 1049. RETURN 1050. END 1101. SUBROUTINE COND 1102. C**** CALCULATES THE CANOPY CONDUCTANCE 1103. C**** INPUT: 1104. C**** BETAD - BETA DUE TO ROOTS 1105. C**** ALAIE - EFFECTIVE LEAF AREA INDEX 1106. C**** RS - MINIMUM STOMATAL RESISTANCE, M S-1 1107. C**** XINC - INCOMING SOLAR RADIATION, W M-2 1108. C**** TP - TEMPERATURE OF CANOPY, C 1109. C**** TFRZ - FREEZING POINT OF WATER, 0 C IN K 1110. C**** OUTPUT: 1111. C**** CNC - CANOPY CONDUCTANCE, M S-1 1112. INCLUDE 'SOILS035.COM' 1113. C**** SOILS28 Common block 9/25/90 1114. ZERO=0. 1115. C**** ADJUST CANOPY CONDUCTANCE FOR SOIL WATER POTENTIAL 1116. CNC=BETAD*ALAIE/RS 1117. C**** ADJUST CANOPY CONDUCTANCE FOR INCOMING SOLAR RADIATION 1118. SRHT0=MAX(SRHT,ZERO) 1119. CNC=CNC*(SRHT0/C1)/(1.+SRHT0/C1) 1120. CNC=CNC/(1.+((TP(0,2)+TFRZ-296.)/15.)**4) 1121. RETURN 1122. END 1201. SUBROUTINE RUNOFF 1202. C**** CALCULATES SURFACE AND UNDERGROUND RUNOFFS. 1203. C**** INPUT: 1204. C**** PRE - EFFECTIVE PRECIPITATION, M S-1 1205. C**** SNOWF - SNOW FALL, EQUIVALENT WATER M S-1 1206. C**** EVAP - EVAPORATION, M S-1 1207. C**** DR - CANOPY DRIP, M S-1 1208. C**** XINFC - INFILTRATION CAPACITY, M S-1 1209. C**** PRFR - FRACTION OF PRECIPITATION 1210. C**** XK - CONDUCTIVITY, M S-1 1211. C**** DZ - LAYER THICKNESSES, M 1212. C**** SL - SLOPE 1213. C**** SDSTNC - INTERSTREAM DISTANCE, M 1214. C**** OUTPUT: 1215. C**** RNF - SURFACE RUNOFF 1216. C**** RNFF - UNDERGROUND RUNOFF, M S-1 1217. INCLUDE 'SOILS035.COM' 1218. C**** SOILS28 Common block 9/25/90 1219. C USE EFFECTS OF SUBGRID SCALE RAIN 1220. C USE PRECIPITATION THAT INCLUDES SMOW MELT 1221. DIMENSION PTMP(2),PTMPS(2) 1222. ZERO=0. 1223. PTMPS(1)=PRS+PRE(1)-SNOWFS 1224. PTMP(1)=PR-PRS-(SNOWF-SNOWFS) 1225. PTMPS(2)=DRS+PRE(2)-SNOWFS 1226. PTMP(2)=DR-DRS-(SNOWF-SNOWFS) 1227. DO 10 IBV=1,2 1228. RNFS=MAX(PTMPS(IBV)-XINFC(IBV),ZERO) 1229. RNF(IBV)=RNFS 1230. IF(PTMP(IBV).GT.0.)THEN 1231. PRFAC=(XINFC(IBV)-PTMPS(IBV))*PRFR/PTMP(IBV) 1232. IF(PRFAC.GE.0.)THEN 1233. IF(PRFAC.LT.30.000) RNF(IBV)=PTMP(IBV)*EXP(-PRFAC) 1234. ELSE 1235. RNF(IBV)=PTMP(IBV)+PTMPS(IBV)-XINFC(IBV) 1236. ENDIF 1237. ENDIF 1238. 10 CONTINUE 1239. C UNDERGROUND RUNOFF 1240. C SL IS THE SLOPE, SDSTNC IS THE INTERSTREAM DISTANCE 1241. DO 20 IBV=1,2 1242. DO 20 L=1,N 1243. RNFF(L,IBV)=XKU(L,IBV)*SL*DZ(L)/SDSTNC 1244. 20 CONTINUE 1245. RETURN 1246. END 1301. SUBROUTINE SINK 1302. C**** CALCULATES WATER SINKS FROM EACH SOIL LAYER 1303. C**** INPUT: 1304. C**** RNF - SURFACE RUNOFF, M S-1 1305. C**** RNFF - UNDERGROUND RUNOFF, M S-1 1306. C**** EVAPD - EVAPORATION FROM DRY CANOPY, M S-1 1307. C**** FR - FRACTION OF ROOTS IN LAYERS 1308. C**** OUTPUT: 1309. C**** SNK - WATER SINK FROM LAYERS, M S-1 1310. INCLUDE 'SOILS035.COM' 1311. C**** SOILS28 Common block 9/25/90 1312. C**** UNDERGROUND RUNOFF IS A SINK 1313. DO 10 IBV=1,2 1314. DO 10 L=1,N 1315. 10 SNK(L,IBV)=RNFF(L,IBV) 1316. C**** REMOVE TRANSPIRED WATER FROM SOIL LAYERS 1317. DO 15 L=1,N 1318. 15 SNK(L,2)=SNK(L,2)+EVAPD*BETADL(L)/(BETAD+1.E-12) 1319. C**** INCLUDE EFFECTS OF SURFACE RUNOFF IN SINK FROM FIRST SOIL LAYERS 1320. DO 20 IBV=1,2 1321. 20 SNK(1,IBV)=SNK(1,IBV)+RNF(IBV) 1322. SNK(0,2)=0. 1323. RETURN 1324. END 1401. SUBROUTINE FLLMT 1402. C**** PLACES LIMITS ON THE SOIL WATER FLUXES 1403. C**** INPUT: 1404. C**** W - WATER IN LAYERS, M 1405. C**** WS - SATURATED WATER IN LAYERS, M 1406. C**** DTS - CURRENT TIME STEP SIZE, S 1407. C**** F - WATER FLUXES, M S-1 1408. C**** SNK - WATER SINK FROM LAYERS, M S-1 1409. C**** RNF - SURFACE RUNOFF, M S-1 1410. C**** SNOWD - SNOW DEPTH, EQUIVALENT WATER M 1411. C**** SNOWF - SNOW FALL, EQUIVALENT WATER M S-1 1412. C**** OUTPUT: 1413. C**** F - LIMITED WATER FLUXES, M S-1 1414. C**** SNK - LIMITED WATER SINKS, M S-1 1415. C**** RNF - LIMITED SURFACE RUNOFF, M S-1 1416. C**** TEMP VARIABLES: 1417. C**** SNOWDU - THE UPPER BOUND ON THE SNOW DEPTH AT END OF TIME STEP 1418. C**** SNOWDL - THE LOWER BOUND ON THE SNOW DEPTH AT END OF TIME STEP 1419. C**** TRUNC - FIX FOR TRUNCATION ON IBM MAINFRAMES 1420. INCLUDE 'SOILS035.COM' 1421. C**** SOILS28 Common block 9/25/90 1422. DIMENSION SNOWDU(2) 1423. ZERO=0. 1424. TRUNC=1.E-6 1425. TRUNC=1.E-12 1426. C PREVENT OVER/UNDERSATURATION OF LAYERS 2-N 1427. SNOWDU(1)=SNOWD(1) 1428. SNOWDU(2)=SNOWD(2) 1429. IF(HT(1,1).LT.0)SNOWDU(1)=SNOWDU(1)+(SNOWF-EVAP(1))*DTS 1430. IF(HT(1,2).LT.0)SNOWDU(2)=SNOWDU(2)+(SNOWF-EVAPS )*DTS 1431. DO 10 IBV=1,2 1432. LL=2-IBV 1433. SNOWDU(IBV)=MAX(ZERO,SNOWDU(IBV)) 1434. DO 10 L=N,2,-1 1435. FLMT=(W(L,IBV)-WS(L,IBV)+TRUNC)/DTS+F(L+1,IBV)-SNK(L,IBV) 1436. F(L,IBV)=MAX(F(L,IBV),FLMT) 1437. FLMT=(W(L,IBV)-DZ(L)*THETM(L,IBV)-TRUNC)/DTS+F(L+1,IBV)-SNK(L,IBV) 1438. F(L,IBV)=MIN(F(L,IBV),FLMT) 1439. 10 CONTINUE 1440. C PREVENT OVER/UNDERSATURATION OF FIRST LAYER 1441. C W(1) CAN INCLUDE SNOW LAYER. 1442. C BARE SOIL 1443. FLMT=(W(1,1)-WS(1,1)-SNOWDU(1)+TRUNC)/DTS+F(2,1)-SNK(1,1) 1444. DRNF=MAX(ZERO,FLMT-F(1,1)) 1445. RNF(1)=RNF(1)+DRNF 1446. SNK(1,1)=SNK(1,1)+DRNF 1447. FLMT=(W(1,1)-SNOWDU(1)-DZ(1)*THETM(1,1)-TRUNC)/DTS+F(2,1)-SNK(1,1) 1448. DRNF=MIN(ZERO,FLMT-F(1,1)) 1449. RNF(1)=RNF(1)+DRNF 1450. SNK(1,1)=SNK(1,1)+DRNF 1451. C PREVENT OVER/UNDERSATURATION OF CANOPY LAYER 1452. FLMT=(WS(0,2)-W(0,2)-TRUNC)/DTS+F(0,2)+SNK(0,2) 1453. F(1,2)=MIN(FLMT,F(1,2)) 1454. FLMT=(-W(0,2)+TRUNC)/DTS+F(0,2)+SNK(0,2) 1455. F(1,2)=MAX(FLMT,F(1,2)) 1456. DR=-F(1,2) 1457. DR=MAX(ZERO,DR) 1458. C PREVENT OVER/UNDERSATURATION OF FIRST LAYER 1459. C VEGETATED SOIL 1460. FLMT=(W(1,2)-WS(1,2)-SNOWDU(2)+TRUNC)/DTS+F(2,2)-SNK(1,2) 1461. DRNF=MAX(ZERO,FLMT-F(1,2)) 1462. RNF(2)=RNF(2)+DRNF 1463. SNK(1,2)=SNK(1,2)+DRNF 1464. FLMT=(W(1,2)-SNOWDU(2)-DZ(1)*THETM(1,2)-TRUNC)/DTS+F(2,2)-SNK(1,2) 1465. DRNF=MIN(ZERO,FLMT-F(1,2)) 1466. RNF(2)=RNF(2)+DRNF 1467. SNK(1,2)=SNK(1,2)+DRNF 1468. RETURN 1469. END 1601. SUBROUTINE FHLMT 1602. C**** MODIFIES SOIL HEAT FLUXES TO ELIMINATE POSSIBLE 1603. C**** OSCILLATION IN PRESENCE OF VARYING COEFFICIENT OF DRAG. 1604. C**** INPUT: 1605. C**** FH - HEAT FLUXES 1606. C**** SHC - HEAT CAPACITIES 1607. C**** W - WATER AMOUNTS 1608. C**** FICE - ICE FRACTION 1609. C**** DT - EXTERNAL TIME STEP 1610. C**** OUTPUT: 1611. C**** FH - CORRECTED HEAT FLUXES 1612. C**** PARAMETER: 1613. C**** DTPL - THE MAX TEMPERATURE CHANGE IN A TIME STEP 1614. C**** 1615. C**** ADD EXCESS FLUX TO FLUX OF LAYER BELOW 1616. INCLUDE 'SOILS035.COM' 1617. C**** SOILS28 Common block 9/25/90 1618. DTPL=1.0 1619. DO 100 IBV=1,2 1620. LL=2-IBV 1621. DO 100 L=LL,N-1 1622. CC=(SHC(L,IBV)+W(L,IBV)*(FICE(L,IBV)*SHI+(1.-FICE(L,IBV))*SHW)) 1623. GM=CC*DTPL/DT 1624. DFH=FH(L+1,IBV)-FH(L,IBV) 1625. IF(ABS(DFH).GT.GM)FH(L+1,IBV)=FH(L,IBV)+SIGN(GM,DFH) 1626. 100 CONTINUE 1627. RETURN 1628. END 1701. SUBROUTINE XKLH 1702. C**** EVALUATES THE HEAT CONDUCTIVITY BETWEEN LAYERS 1703. C**** USES THE METHOD OF DEVRIES. 1704. C**** INPUT: 1705. C**** ZB - SOIL LAYER BOUNDARIES, M 1706. C**** ZC - SOIL LAYER CENTERS, M 1707. C**** THETA - SOIL WATER SATURATION 1708. C**** FICE - FRACTION OF ICE IN LAYERS 1709. C**** ALAMI - ICE HEAT CONDUCTIVITY 1710. C**** ALAMW - WATER HEAT CONDUCTIVITY 1711. C**** TP - TEMPERATURE OF LAYERS, C 1712. C**** SHW - SPECIFIC HEAAT OF WATER 1713. C**** SHI - SPECIFIC HEAT OF ICE 1714. C**** SHC - HEAT CAPACITY OF SOIL LAYERS 1715. C**** DZ - LAYER THICKNESSES 1716. C**** L,IBV - SOIL LAYER 1717. C**** DTS - THE CURRENT TIME STEP 1718. C**** OUTPUT: 1719. C**** XKH(L,IBV) - HEAT CONDUCTIVITIES IN EACH OF THE SOIL LAYERS 1720. C**** XKHM(L,IBV) - AVERAGE HEAT CONDUCTIVITY BETWEEN LAYER L AND L-1 1721. C**** 1722. INCLUDE 'SOILS035.COM' 1723. C**** SOILS28 Common block 9/25/90 1724. DIMENSION XSHA(NG,2),XSH(NG,2),GABC(3),HCWT(IMT-1) 1725. C 1726. C CALCULATE WITH CHANGING GA FOR AIR. GA IS THE DEPOLARIZATION 1727. C FACTOR FOR AIR, CALCULATED BY LINEAR INTERPOLATION FROM .333 1728. C AT SATURATION TO .035 AT 0 WATER, FOLLOWING DEVRIES. 1729. DO 2 IBV=1,2 1730. DO 2 L=1,N 1731. GAA=.298*THETA(L,IBV)/(THETS(L,IBV)+1.E-6)+.035 1732. GCA=1.-2.*GAA 1733. HCWTA=(2./(1.+BA*GAA)+1./(1.+BA*GCA))/3. 1734. C XW,XI,XA ARE THE VOLUME FRACTIONS. DON'T COUNT SNOW IN SOIL LAYER 1 1735. XW=W(L,IBV)*(1.-FICE(L,IBV))/DZ(L) 1736. XI=W(L,IBV)*FICE(L,IBV)/DZ(L) 1737. IF(L.EQ.1)XI=XI-SNOWD(IBV)/DZ(L) 1738. XA=(THETS(L,IBV)-THETA(L,IBV)) 1739. XB=Q(IMT,L) 1740. XNUM=XW*HCWTW*ALAMW+XI*HCWTI*ALAMI+XA*HCWTA*ALAMA+XSHA(L,IBV) 1741. & + XB*HCWTB*ALAMBR 1742. XDEN=XW*HCWTW+XI*HCWTI+XA*HCWTA+XSH(L,IBV)+XB*HCWTB 1743. XKH(L,IBV)=XNUM/XDEN 1744. 2 CONTINUE 1745. C PUT THE SNOW CONDUCTIVITY IN SERIES WITH THE FIRST GROUND LAYER. 1746. DO 3 IBV=1,2 1747. SNDP=SNOWD(IBV)/SPGSN 1748. XKH(1,IBV)=(DZ(1)+SNDP)/(DZ(1)/XKH(1,IBV)+SNDP/ALAMSN) 1749. 3 CONTINUE 1750. C GET THE AVERAGE CONDUCTIVITY BETWEEN LAYERS 1751. DO 4 IBV=1,2 1752. DO 4 L=2,N 1753. XKHM(L,IBV)=((ZB(L)-ZC(L-1))*XKH(L,IBV) 1754. & + (ZC(L)-ZB(L))*XKH(L-1,IBV) 1755. & )/(ZC(L)-ZC(L-1)) 1756. 4 CONTINUE 1757. C**** 1758. RETURN 1759. ENTRY XKLH0 1760. C GABC'S ARE THE DEPOLARIZATION FACTORS, OR RELATIVE SPHEROIDAL AXES. 1761. GABC(1)=.125 1762. GABC(2)=GABC(1) 1763. GABC(3)=1.-GABC(1)-GABC(2) 1764. C HCWT'S ARE THE HEAT CONDUCTIVITY WEIGHTING FACTORS 1765. HCWTW=1. 1766. HCWTI=0. 1767. HCWTB=1. 1768. DO 10 I=1,IMT-1 1769. HCWT(I)=0. 1770. 10 CONTINUE 1771. DO 20 J=1,3 1772. HCWTI=HCWTI+1./(1.+(ALAMI/ALAMW-1.)*GABC(J)) 1773. DO 20 I=1,IMT-1 1774. HCWT(I)=HCWT(I)+1./(1.+(ALAMS(I)/ALAMW-1.)*GABC(J)) 1775. 20 CONTINUE 1776. HCWTI=HCWTI/3. 1777. DO 30 I=1,IMT-1 1778. HCWT(I)=HCWT(I)/3. 1779. 30 CONTINUE 1780. DO 25 IBV=1,2 1781. DO 25 L=1,N 1782. XSHA(L,IBV)=0. 1783. XSH(L,IBV)=0. 1784. DO 25 I=1,IMT-1 1785. XS=(1.-THM(0,I))*Q(I,L) 1786. XSHA(L,IBV)=XSHA(L,IBV)+XS*HCWT(I)*ALAMS(I) 1787. XSH(L,IBV)=XSH(L,IBV)+XS*HCWT(I) 1788. 25 CONTINUE 1789. BA=ALAMA/ALAMW-1. 1790. C 1791. RETURN 1792. END 1901. SUBROUTINE FLH 1902. C**** EVALUATES THE HEAT FLUX BETWEEN LAYERS 1903. C**** SUBROUTINE FL MUST BE CALLED FIRST 1904. C**** INPUT: 1905. C**** ZB - SOIL LAYER BOUNDARIES, M 1906. C**** ZC - SOIL LAYER CENTERS, M 1907. C**** THETA - SOIL WATER SATURATION 1908. C**** FICE - FRACTION OF ICE IN LAYERS 1909. C**** ALAMI - ICE HEAT CONDUCTIVITY 1910. C**** ALAMW - WATER HEAT CONDUCTIVITY 1911. C**** TP - TEMPERATURE OF LAYERS, C 1912. C**** SHW - SPECIFIC HEAAT OF WATER 1913. C**** OUTPUT: 1914. C**** FH - HEAT FLUX BETWEEN LAYERS 1915. INCLUDE 'SOILS035.COM' 1916. C**** SOILS28 Common block 9/25/90 1917. C**** 1918. DO 20 IBV=1,2 1919. FH(N+1,IBV)=0. 1920. C TOTAL HEAT FLUX IS HEAT CARRIED BY WATER FLOW PLUS HEAT CONDUCTION 1921. DO 20 L=2,N 1922. C TAKE INTO ACCOUNT SNOW DEPTH FOR FIRST LAYER 1923. IF(L.EQ.2)THEN 1924. SNC=.5*SNOWD(IBV)/SPGSN 1925. ELSE 1926. SNC=0. 1927. ENDIF 1928. FH(L,IBV)=-XKHM(L,IBV)*(TP(L-1,IBV)-TP(L,IBV))/(SNC+ZC(L-1)-ZC(L)) 1929. IF(F(L,IBV).GT.0)THEN 1930. FH(L,IBV)=FH(L,IBV)+F(L,IBV)*TP(L,IBV)*SHW 1931. ELSE 1932. FH(L,IBV)=FH(L,IBV)+F(L,IBV)*TP(L-1,IBV)*SHW 1933. ENDIF 1934. 20 CONTINUE 1935. RETURN 1936. END 2001. SUBROUTINE FLHG 2002. C**** CALCULATES THE GROUND HEAT FLUXES (TO THE SURFACE) 2003. C**** INPUT: 2004. C**** OUTPUT: 2005. C**** FH - HEAT FLUXES FROM BARE SOIL SURFACE, AND FROM CANOPY, 2006. C**** AND BETWEEN CANOPY AND VEGETATED SOIL. 2007. C**** AFHG - HEAT FLUX FROM GROUND TO CANOPY 2008. INCLUDE 'SOILS035.COM' 2009. C**** SOILS28 Common block 9/25/90 2010. PARAMETER (STBO=5.67032E-8) 2011. C 2012. C**** BARE SOIL FLUXES 2013. IBV=1 2014. L=2-IBV 2015. THRM(IBV)=STBO*(TP(L,IBV)+TFRZ)**4 2016. FH(L,IBV)=XLTH(IBV)+SNSH(IBV) 2017. FH(L,IBV)=FH(L,IBV)-HTPR 2018. FH(L,IBV)=FH(L,IBV)+THRM(IBV)-SRHT-TRHT 2019. C 2020. C**** CANOPY FLUXES, AND FLUXES FROM MASKING SNOW 2021. IBV=2 2022. L=2-IBV 2023. THRM(IBV) = STBO*(TCS+TFRZ)**4 2024. FH(L,IBV) = XLTH(IBV) + SNSH(IBV) 2025. FH(L,IBV) = FH(L,IBV) - HTPR 2026. FH(L,IBV) = FH(L,IBV) + (THRM(IBV)-SRHT-TRHT) 2027. C FH(L+1,IBV) = FM*FH(L,IBV) 2028. FH(L+1,IBV) = FM*(THRM(IBV)-SRHT-TRHT) 2029. C**** 2030. FH(1,2)=FH(1,2)-SHW*DR*TP(0,2) 2031. FH(1,2)=FH(1,2)-STBO*((TP(0,2)+TFRZ)**4-(TP(1,2)+TFRZ)**4) 2032. RETURN 2033. END 2101. SUBROUTINE SINKH 2102. C**** CALCULATES THE HEAT REMOVAL FROM EACH LAYER 2103. C**** INPUT: 2104. C**** SHW - SPECIFIC HEAT OF WATER 2105. C**** TP - TEMPERATURE OF LAYERS, C 2106. C**** SNK - SOIL WATER SINK IN LAYERS, M S-1 2107. C**** OUTPUT: 2108. C**** SNKH - HEAT SINK FROM SOIL LAYERS 2109. INCLUDE 'SOILS035.COM' 2110. C**** SOILS28 Common block 9/25/90 2111. DO 10 IBV=1,2 2112. DO 10 L=1,N 2113. 10 SNKH(L,IBV)=SHW*TP(L,IBV)*SNK(L,IBV) 2114. RETURN 2115. END 2201. SUBROUTINE RETP 2202. C**** EVALUATES THE TEMPERATURES IN THE SOIL LAYERS BASED ON THE 2203. C**** HEAT VALUES. ALSO EXECUTES SNOW MELT. 2204. C**** INPUT: 2205. C**** W - WATER IN SOIL LAYERS, M 2206. C**** HT - HEAT IN SOIL LAYERS 2207. C**** FSN - HEAT OF FUSION OF WATER 2208. C**** SHC - SPECIFIC HEAT CAPACITY OF SOIL 2209. C**** SHI - SPECIFIC HEAT CAPACITY OF ICE 2210. C**** SHW - SPECIFIC HEAT CAPCITY OF WATER 2211. C**** SNOWD - SNOW DEPTH, EQUIVALENT WATER M 2212. C**** FB - FRACTION OF BARE SOIL 2213. C**** FV - FRACTION OF VEGETATION 2214. C**** FM - SNOW VEGETATION MASKING FRACTION (REQUIRES RETH CALLED FIRST) 2215. C**** OUTPUT: 2216. C**** TP - TEMPERATURE OF LAYERS, C 2217. C**** FICE - FRACTION OF ICE OF LAYERS 2218. C**** TBCS - TEMPERATURE OF BARE SOIL, CANOPY, AND SNOW AS SEEN 2219. C**** BY ATMOSPHERE, C. ALSO CALLED GROUND TEMPERATURE. 2220. C**** TCS - TEMPERATURE OF CANOPY AND SNOW, C. 2221. INCLUDE 'SOILS035.COM' 2223. DO 30 IBV=1,2 2224. DO 30 L=2-IBV,N 2225. C**** Ground temperature is > 0 2226. IF(HT(L,IBV).lt.0.) GO TO 10 2227. TP(L,IBV) = HT(L,IBV)/(SHC(L,IBV)+W(L,IBV)*SHW) 2228. FICE(L,IBV) = 0. 2229. GO TO 30 2230. C**** Ground temperature is = 0 2231. 10 IF(FSN*W(L,IBV)+HT(L,IBV).lt.0.) GO TO 20 2232. TP(L,IBV) = 0. 2233. FICE(L,IBV) = -HT(L,IBV)/(FSN*W(L,IBV)) 2234. GO TO 30 2235. C**** Ground temperature is < 0 2236. 20 TP(L,IBV) = (HT(L,IBV)+W(L,IBV)*FSN)/(SHC(L,IBV)+W(L,IBV)*SHI) 2237. FICE(L,IBV) = 1. 2238. 30 CONTINUE 2240. TCS = (1.-FM)*TP(0,2) + FM*TP(1,2) 2241. TBCS = FB*TP(1,1) + FV*TCS 2242. C**** 2243. IF(TP(1,1).le.100. .and. TP(0,2).le.100.) RETURN 2244. WRITE (6,*) 'RETP TP BOUNDS ERROR' 2245. WRITE (6,*) 'ID',ID 2246. CALL RETH 2247. CALL HYDRA 2248. CALL OUTW (1) 2249. STOP 'TP' 2250. END 2301. SUBROUTINE SNMLT 2302. INCLUDE 'SOILS035.COM' 2303. C**** SOILS28 Common block 9/25/90 2304. C**** EXECUTE SMOW MELT - snow melts before ice in layer 1 2305. C**** ADD MELT TO EFFECTIVE PRECIPITATION, FOR PURPOSES OF RUNOFF 2306. DO 100 IBV=1,2 2307. L=1 2308. PRE(IBV)=0. 2309. IF(FICE(L,IBV).lt.1. .and. SNOWD(IBV).gt.0.) THEN 2309.1 XLW = (1.-FICE(L,IBV))*W(L,IBV) 2309.2 XLW = MIN (XLW, SNOWD(IBV)) 2310. PRE(IBV) = XLW/DTS 2311. SNOWD(IBV) = SNOWD(IBV) - XLW 2312. ENDIF 2313. 100 CONTINUE 2314. RETURN 2315. END 2401. SUBROUTINE ADVNC 2402. C**** ADVANCES QUANTITIES BY ONE TIME STEP. 2403. C**** INPUT: 2404. C**** DT - TIME STEP, S 2405. C**** DZ - LAYER THICKNESS, M 2406. C**** TP - LAYER TEMPERATURES, C 2407. C**** TFRZ - FREEZING POINT OF WATER, K 2408. C**** W - SOIL WATER IN LAYERS, M 2409. C**** SNOWD - SNOW DEPTH, M 2410. C**** F - WATER FLUX, M S-1 2411. C**** SNK - WATER SINKS, M S-1 2412. C**** HT - HEAT IN SOIL LAYERS 2413. C**** FH - HEAT FLUX IN SOIL LAYERS 2414. C**** SNKH - HEAT SINK IN LAYERS 2415. C**** SNOWF - SNOW FALL, M S-1 OF EQUIVALENT WATER 2416. C**** EVAP - EVAPORATION, M S-1 2417. C**** OUTPUT: 2418. C**** W - UPDATER WATER IN SOIL LAYERS, M S-1 2419. C**** HT - UPDATED HEAT IN SOIL LAYERS 2420. C**** SNOWD - UPDATED SNOW DEPTH, M S-1 OF EQUIVALENT WATER 2421. C**** RUS - OVERALL SURFACE RUNOFF, M S-1 REPLACED BY ARUNS 2422. C**** ARUNS - OVERALL SURFACE RUNOFF, KG M-2 2423. C**** AERUNS - OVERALL SURFACE HEAT RUNOFF, J M-2 2424. C**** AERUNU - UNDERGROUND HEAT RUNOFF, J M-2 2425. C**** USES: 2426. C**** RETP,RETH,FL,FLG,RUNOFF,SINK,SINKH,FLLMT,FLH,FLHG. 2427. C**** ALSO USES SURF WITH ITS REQUIRED VARIABLES. 2428. INCLUDE 'SOILS035.COM' 2429. PARAMETER (AMASS1=500., RHOW=1000.) 2430. ZERO=0. 2431. LIMIT=1800 2432. NIT=0 2433. DTR=DT 2434. DR=0. 2435. NINTEG=0 2436. CALL RETH 2437. CALL RETP 2438. C CALL HYDRA 2439. C**** 2440. TB0=TP(1,1) 2441. TC0=TP(0,2) 2442. C**** 2443. 20 CONTINUE 2444. C**** 2445. NIT=NIT+1 2446. IF(NIT.GT.LIMIT)GO TO 900 2447. CALL HYDRA 2448. CALL WTAB 2449. CALL QSBAL 2450. CALL XKLH 2451. CALL GDTM(DTM) 2452. DTS=MIN(DTR,DTM) 2453. DTR=DTR-DTS 2454. NINTEG=NINTEG+1 2455. CALL SNMLT 2456. CALL FL 2457. CALL FLG 2458. CALL RUNOFF 2459. CALL SINK 2460. CALL FLLMT 2461. CALL SINKH 2462. CALL FLH 2463. CALL FLHG 2464. C CALL FHLMT 2464.1 TS = TS + (FB*SNSH(1)+FV*SNSH(2) )*DTS / (SHA*AMASS1) 2464.2 QS = QS + (FB*EVAP(1)+FV*(EVAPW+EVAPD))*DTS*RHOW / AMASS1 2465. C**** 2466. IF(HT(1,1).LT.0.)SNOWD(1)=SNOWD(1)+(SNOWF-EVAP(1))*DTS 2467. IF(HT(1,2).LT.0.)SNOWD(2)=SNOWD(2)+(SNOWF-EVAPS )*DTS 2468. DO 70 IBV=1,2 2469. LL=2-IBV 2470. SNOWD(IBV)=MAX(ZERO,SNOWD(IBV)) 2471. DO 60 L=LL,N 2472. W(L,IBV)=W(L,IBV)+(F(L+1,IBV)-F(L,IBV)-SNK(L,IBV))*DTS 2473. 60 HT(L,IBV)=HT(L,IBV)+(FH(L+1,IBV)-FH(L,IBV)-SNKH(L,IBV))*DTS 2474. 70 CONTINUE 2475. C**** 2476. W(0,2)=MAX(W(0,2),ZERO) 2477. CALL ACCM 2478. CALL RETH 2479. CALL RETP 2480. C CALL HYDRA 2481. IF(DTR.gt.0.) GO TO 20 2484. CALL ACCMF 2485. C CALL WTAB 2486. CALL OUTGH 2487. C**** 2488. RETURN 2489. 900 CONTINUE 2490. WRITE(6,*)'LIMIT EXCEEDED' 2491. WRITE(6,*)'DTR,DTM,DTS',DTR,DTM,DTS 2492. WRITE(6,*)'TB0,TC0',TB0,TC0 2493. CALL OUTW(2) 2494. STOP 'ADVNC' 2495. END 2601. SUBROUTINE ACCM 2602. C**** ACCUMULATES GCM DIAGNOSTICS 2603. INCLUDE 'SOILS035.COM' 2604. C**** SOILS28 Common block 9/25/90 2605. C**** MIXING RATIO FOR GIVEN TEMPERATURE TX(C) AND PRESSURE PX(MB) 2606. QSAT(TX,PX)=.622 * 33.8639*( (0.00738*TX+0.8072)**8 - 2607. * 0.000019*ABS(1.8*TX+48)+0.001316 )/PX 2608. C**** 2609. C**** THE FOLLOWING LINES WERE ORIGINALLY CALLED BEFORE RETP, 2610. C**** RETH, AND HYDRA. 2611. AERUNS=AERUNS+SHW*(FB*TP(1,1)*RNF(1)+FV*TP(1,2)*RNF(2))*DTS 2612. ADIFS=ADIFS-DTS*(F(2,1)*FB+F(2,2)*FV) 2613. DEDIFS=F(2,1)*TP(2,1) 2614. IF(F(2,1).LT.0.) DEDIFS=F(2,1)*TP(1,1) 2615. AEDIFS=AEDIFS-DTS*SHW*DEDIFS*FB 2616. DEDIFS=F(2,2)*TP(2,2) 2617. IF(F(2,2).LT.0.) DEDIFS=F(2,2)*TP(1,2) 2618. AEDIFS=AEDIFS-DTS*SHW*DEDIFS*FV 2619. C ALHG=ALHG+DTS*(EVAP(1)*(ELH+TP(1,1)*SHV)*FB+ 2620. C * EVAP(2)*(ELH+TP(0,2)*SHV)*FV) 2621. ALHG=ALHG+(XLTH(1)*FB+XLTH(2)*FV)*DTS 2622. C**** 2623. C**** THE FOLOWING LINES WERE ORIGINALLY CALLED AFTER RETP, 2624. C**** RETH, AND HYDRA. 2625. AF0DT=AF0DT-DTS*(FB*FH(1,1)+FV*FH(0,2)+HTPR) 2626. AF1DT=AF1DT-DTS*(FB*FH(2,1)+FV*FH(2,2)) 2627. ARUNS=ARUNS+(FB*RNF(1)+FV*RNF(2))*DTS 2628. DO 74 L=1,N 2629. ARUNU=ARUNU+(RNFF(L,1)*FB+RNFF(L,2)*FV)*DTS 2630. AERUNU=AERUNU+(SNKH(L,1)*FB+SNKH(L,2)*FV)*DTS 2631. 74 CONTINUE 2632. ASHG=ASHG+(SNSH(1)*FB+SNSH(2)*FV)*DTS 2633. ATRG=ATRG+(THRM(1)*FB+THRM(2)*FV)*DTS 2634. C**** 2635. AEVAPW=AEVAPW+(EVAPW*FV*DTS) 2636. AEVAPD=AEVAPD+(EVAPD*FV*DTS) 2637. AEVAPB=AEVAPB+(EVAP(1)*FB*DTS) 2638. AEPC=AEPC+(EPC*FV*DTS) 2639. AEPB=AEPB+(EPB*FB*DTS) 2640. AFHG=AFHG+(FH(1,2)*FV*DTS) 2641. C 2642. RETURN 2643. ENTRY ACCMF 2644. C PROVIDES ACCUMULATION UNITS FIXUPS, AND CALCULATES 2645. C PENMAN EVAPORATION. SHOULD BE CALLED ONCE AFTER 2646. C ACCUMULATIONS ARE COLLECTED. 2647. ONE=1. 2648. ZERO=0. 2649. ARUNS=1000.0*ARUNS 2650. ARUNU=1000.0*ARUNU 2651. AEVAPW=1000.0*AEVAPW 2652. AEVAPD=1000.0*AEVAPD 2653. AEVAPB=1000.0*AEVAPB 2654. AEPC=1000.0*AEPC 2655. AEPB=1000.0*AEPB 2656. ADIFS=1000.*ADIFS 2657. AF1DT=AF1DT-AEDIFS 2658. C**** CALCULATION OF PENMAN VALUE OF POTENTIAL EVAPORATION, AEPP 2659. H0=FB*(SNSH(1)+XLTH(1))+FV*(SNSH(2)+XLTH(2)) 2660. C H0=-ATRG/DT+SRHT+TRHT 2661. CCC H0=-THRM(2)+SRHT+TRHT 2662. EL0=ELH*1.E-3 2663. CNA=CH*VSM 2664. CPFAC=SHA*RHO*CNA 2665. T0=TS-TFRZ 2666. EDELT=100.*PRES*(QSAT(T0,PRES)-QS)/0.622 2667. GAMMA=SHA*100.*PRES/(0.622*EL0) 2668. IF(1.8*T0+48.0 .LT. 0.) THEN 2669. DELT=33.8639*(8.*0.00738*(0.00738*T0+0.8072)**7.+0.000019*1.8) 2670. * *100.0 2671. ELSE 2672. DELT=33.8639*(8.*0.00738*(0.00738*T0+0.8072)**7.-0.000019*1.8) 2673. * *100.0 2674. END IF 2675. EPEN=(DELT*H0+CPFAC*EDELT)/(EL0*(DELT+GAMMA)) 2676. AEPP=EPEN*DT 2677. ABETAP=1. 2678. IF(AEPP.GT.0.)ABETAP=(AEVAPW+AEVAPD+AEVAPB)/AEPP 2679. ABETAP=MIN(ABETAP,ONE) 2680. ABETAP=MAX(ABETAP,ZERO) 2681. RETURN 2682. END 2801. SUBROUTINE GDTM(DTM) 2802. C**** CALCULATES THE MAXIMUM TIME STEP ALLOWED BY STABILITY 2803. C**** CONSIDERATIONS. 2804. INCLUDE 'SOILS035.COM' 2805. C**** SOILS28 Common block 9/25/90 2806. PARAMETER (STBO=5.67032E-8) 2807. DIMENSION QG(2),XK2(2),AK2(2) 2808. T450=450. 2809. T0=TS-TFRZ 2810. IF(1.8*T0+48.0 .LT. 0.) THEN 2811. DELT=33.8639*(8.*0.00738*(0.00738*T0+0.8072)**7.+0.000019*1.8) 2812. * *100.0 2813. ELSE 2814. DELT=33.8639*(8.*0.00738*(0.00738*T0+0.8072)**7.-0.000019*1.8) 2815. * *100.0 2816. END IF 2817. DQDT=.622*DELT/(100.*PRES) 2818. QG(1)=QB 2819. QG(2)=QC 2820. C**** 2821. C**** FIRST CALCULATE TIMESTEP FOR WATER MOVEMENT IN SOIL. 2822. SGMM=1.0 2823. DLDZ2=0. 2824. DO 30 IBV=1,2 2825. DO 30 L=1,N 2826. DLDZ2=MAX(DLDZ2,D(L,IBV)/DZ(L)**2) 2827. 30 CONTINUE 2828. DTM=SGMM/(DLDZ2+1.E-12) 2829. IF(Q(4,1).GT.0.)DTM=MIN(DTM,T450) 2830. DTM1=DTM 2831. C**** 2832. C**** NEXT CALCULATE TIMESTEP FOR HEAT MOVEMENT IN SOIL. 2833. DO 40 IBV=1,2 2834. DO 40 L=1,N 2835. XK1=XKH(L,IBV) 2836. AK1=(SHC(L,IBV)+((1.-FICE(L,IBV))*SHW+FICE(L,IBV)*SHI) 2837. & *W(L,IBV))/DZ(L) 2838. DTM=MIN(DTM,.5*AK1*DZ(L)**2/(XK1+1.E-12)) 2839. 40 CONTINUE 2840. DTM2=DTM 2841. C**** 2842. C**** FINALLY, CALCULATE MAX TIME STEP FOR TOP LAYER BARE SOIL 2843. C**** AND CANOPY INTERACTION WITH SURFACE LAYER. 2844. DO 50 IBV=1,2 2845. L=2-IBV 2846. XK2(IBV)=ABS(SNSH(IBV))/(ABS(TP(L,IBV)-TS+TFRZ)+.01) 2847. & + ABS(XLTH(IBV))/(ABS(QG(IBV)-QS)+1.E-5)*DQDT 2848. & + 8.*STBO*(TP(L,IBV)+TFRZ)**3 2849. AK2(IBV)=SHC(L,IBV)+((1.-FICE(L,IBV))*SHW+FICE(L,IBV)*SHI) 2850. & *W(L,IBV) 2851. DTM=MIN(DTM,AK2(IBV)/(XK2(IBV)+1.E-12)) 2852. IF(IBV.EQ.1)DTM3=DTM 2853. IF(IBV.EQ.2)DTM4=DTM 2854. 50 CONTINUE 2855. IF(DTM.LT.1.)THEN 2856. WRITE(*,*)'DTM',DTM1,DTM2,DTM3,DTM4 2857. WRITE(*,*)'XK2',XK2 2858. WRITE(*,*)'AK2',AK2 2859. WRITE(*,*)'SNSH',SNSH 2860. WRITE(*,*)'XLTH',XLTH 2861. WRITE(*,*)'DQDT',DQDT 2862. WRITE(*,*)'TS,TFRZ',TS,TFRZ 2863. WRITE(*,*)'DLT',TP(1,1)-TS+TFRZ,TP(0,2)-TS+TFRZ 2863.1 WRITE(*,*)'QS',QS 2863.2 WRITE(*,*)'QG',QG 2863.3 WRITE(*,*)'FICE',FICE(1,1),FICE(0,2) 2864. ENDIF 2865. C**** 2866. RETURN 2867. END 2901. SUBROUTINE GHINIT (FDATA,VDATA,DTSURF,SHCLC0) 2902. PARAMETER (IM=72,JM=46) 2903. INCLUDE 'SOILS035.COM' 2904. C**** SOILS28 Common block 9/25/90 2905. PARAMETER (STBO=5.67032E-8) 2909. REAL*8 FDATA(IM,JM,3),VDATA(IM,JM,10) 2912. COMMON/SOILS2/SDATA(IM,JM,11*NGM+1),ALA(3,IM,JM),ACS(3,IM,JM) 2913. * ,AFB(IM,JM),AFR(IM,JM,NGM),AVH(IM,JM) 2914. REAL*8 RSAR(8),ALAMIN(8),ALAMAX(8),LADAY(8),AROOT(8),BROOT(8) 2915. REAL*8 VHGHT(8) 2916. DATA AROOT/ 12.5, 0.9, 0.8, 0.25, 0.25, 0.25, 1.1, .9/ 2917. DATA BROOT/ 1.0, 0.9, 0.4, 2.00, 2.00, 2.00, 0.4, .9/ 2918. DATA ALAMAX/ 1.5, 2.0, 2.5, 4.0, 6.0, 10.0, 8.0, 4.5/ 2919. DATA ALAMIN/ 1.0, 1.0, 1.0, 1.0, 1.0, 8.0, 6.0, 1.0/ 2920. DATA LADAY/ 196, 196, 196, 196, 196, 196, 196, 196/ 2921. DATA RSAR/100., 100., 200., 200., 200., 300., 250., 125./ 2922. DATA VHGHT/ 0.1, 1.5, 5., 15., 20., 30., 25., 1.75/ 2923. C**** TUNDR GRASS SHRUB TREES DECID EVRGR RAINF CULTI 2924. C**** 2925. C**** LADAY(veg type, lat belt) = day of peak LAI 2926. C OLD PEAK LAI: 2ND LINE IS FOR LATITUDES < 23.5 DEG 2927. C**** 1 temperate latitudes 2928. C**** 2 non-temperate latitudes 2929. C DATA LADAY/ 196, 196, 196, 196, 196, 196, 105, 2930. C DATA LADAY/ 196, 288, 288, 288, 288, 196, 105, 2931. C**** 2932. C**** CONTENTS OF ALA(K,I,J), LAI coefficients 2933. C**** 1 AVERAGE LEAF AREA INDEX 2934. C**** 2 REAL AMPLITUDE OF LEAF AREA INDEX 2935. C**** 3 IMAGINARY AMPLITUDE OF LEAF AREA INDEX 2936. C**** 2937. C**** CONTENTS OF ACS(K,I,J), CS coefficients 2938. C**** 1 AVERAGE STOMATAL CONDUCTANCE 2939. C**** 2 REAL AMPLITUDE OF STOMATAL CONDUCTANCE 2940. C**** 3 IMAGINARY AMPLITUDE OF STOMATAL CONDUCTANCE 2941. C**** 2942. C**** CONTENTS OF SDATA(I,J,K): 2943. C**** 1 - NGM DZ(NGM) 2944. C**** NGM+1 - 6*NGM Q(IS,NGM) 2945. C**** 6*NGM+1 - 11*NGM QK(IS,NGM) 2946. C**** 11*NGM+1 SL 2947. C 2948. ONE=1. 2949. 10 CONTINUE 2950. C**** 2951. C**** INITIALIZE CONSTANTS 2952. C**** 2953. C**** Time step for ground hydrology 2954. DT=DTSURF 2955. C**** UNITS ARE MKS 2956. C**** WATER QUANTITIES ARE DENSITY TIMES USUAL VALUES IN MKS 2957. C**** TO GET VOLUMETRIC UNITS 2958. C**** 1M WATER = 1000 KG M-2; 1M3 WATER = 1000 KG 2959. C FSN IS THE HEAT OF FUSION 2960. FSN=3.34 E+8 2961. C ELH IS THE HEAT OF VAPORIZATION 2962. ELH=2.50 E+9 2963. C THE SH'S ARE THE SPECIFIC HEAT CAPACATIES 2964. SHW=4.185 E+6 2965. SHI=2.060 E+6 2966. SHA=1003.4965 2967. SHV=1911. 2968. C THE ALAM'S ARE THE HEAT CONDUCTIVITIES 2969. ALAMW=.573345 2970. ALAMI=2.1762 2971. ALAMA=.025 2972. ALAMSN=0.088 2973. ALAMBR=2.9 2974. ALAMS(1)=8.8 2975. ALAMS(2)=2.9 2976. ALAMS(3)=2.9 2977. ALAMS(4)=.25 2978. C HW IS THE WILTING POINT IN METERS 2979. HW=-100 2980. C TFRZ IS 0 C IN K 2981. TFRZ=273.16 2982. C ZHTB IS DEPTH FOR COMBINING HEAT LAYERS FOR STABILITY 2983. IF(Q(4,1).LT..01)THEN 2984. ZHTB=6. 2985. ELSE 2986. ZHTB=6. 2987. ENDIF 2988. C SPGSN IS THE SPECIFIG GRAVITY OF SNOW 2989. SPGSN=.1 2990. C 2991. CR CALL SOLDAT (FDATA,VDATA) 2992. C**** 2993. C**** Initialize global arrays ALA, ACS, AFB, AFR 2994. C**** 2995. TWOPI=6.283185 2996. FJEQ=(JM+1)/2. 2997. DO 210 I=1,IM*JM*(8+NGM) 2998. 210 ALA(I,1,1)=0. 2999. DO 220 J=1,JM 3000. DO 220 I=1,IM 3001. 220 ACS(1,I,J)=.01 3002. DO 400 J=1,JM 3003. DO 400 I=1,IM 3004. PEARTH=FDATA(I,J,2)*(1.-FDATA(I,J,3)) 3005. AFB(I,J) = VDATA(I,J,1) + VDATA(I,J,10) 3005.1 IF(AFB(I,J).gt..999) AFB(I,J) = 1. 3006. IF(PEARTH.le.0. .or. AFB(I,J).ge.1.) GO TO 400 3007. C**** CALCULATE LAI, CS COEFFICICENTS 3008. SFV=0. 3009. SLA0=0. 3010. SLRE=0. 3011. SLIM=0. 3012. SCS0=0. 3013. SCSRE=0. 3014. SCSIM=0. 3015. SVH=0. 3016. DO 250 IV=1,8 3017. PHASE=TWOPI*LADAY(IV)/365. 3018. IF(J.LT.FJEQ) PHASE=PHASE+TWOPI/2. 3019. FV=VDATA(I,J,IV+1) 3020. SFV=SFV+FV 3021. SVH=SVH+FV*VHGHT(IV) 3022. DIF=(ALAMAX(IV) - ALAMIN(IV)) 3023. SLA0=SLA0+FV*(ALAMAX(IV) + ALAMIN(IV)) 3024. SLRE=SLRE+FV*DIF*COS(PHASE) 3025. SLIM=SLIM+FV*DIF*SIN(PHASE) 3026. SCS0=SCS0+FV*(ALAMAX(IV) + ALAMIN(IV))/RSAR(IV) 3027. SCSRE=SCSRE+FV*DIF*COS(PHASE)/RSAR(IV) 3028. 250 SCSIM=SCSIM+FV*DIF*SIN(PHASE)/RSAR(IV) 3029. ALA(1,I,J)=.5/SFV*SLA0 3030. ALA(2,I,J)=.5/SFV*SLRE 3031. ALA(3,I,J)=.5/SFV*SLIM 3032. ACS(1,I,J)=.5/SFV*SCS0 3033. ACS(2,I,J)=.5/SFV*SCSRE 3034. ACS(3,I,J)=.5/SFV*SCSIM 3035. AVH(I,J)=SVH/SFV 3036. C**** CALCULATE ROOT FRACTION AFR AVERAGED OVER VEGETATION TYPES 3037. DO 310 N=1,NGM 3038. DZ(N)=SDATA(I,J,N) 3039. IF(DZ(N).LE.0.) GO TO 320 3040. 310 CONTINUE 3041. 320 N=N-1 3042. DO 350 IV=1,8 3043. FV=VDATA(I,J,IV+1) 3044. Z=0. 3045. FRUP=0. 3046. DO 350 L=1,N 3047. Z=Z+DZ(L) 3048. FRDN=AROOT(IV)*Z**BROOT(IV) 3049. FRDN=MIN(FRDN,ONE) 3050. IF(L.EQ.N)FRDN=1. 3051. AFR(I,J,L) = AFR(I,J,L) + FV*(FRDN-FRUP) 3052. 350 FRUP=FRDN 3053. DO 370 L=1,N 3054. 370 AFR(I,J,L) = AFR(I,J,L)/(1.-AFB(I,J)) 3055. 400 CONTINUE 3056. C**** 3057. SDSTNC=100. 3058. PRINT *,'SDSTNC:',SDSTNC 3059. C1=90. 3060. PRINT *,'C1:',C1 3061. PRFR=.1 3062. PRINT *,'PRFR:',PRFR 3063. CALL HL0 3064. C**** 3065. RETURN 3066. END 3101. SUBROUTINE GHINIJ (I0,J0, WFC1) 3102. C**** INPUT: 3103. C**** AVH(I,J) - ARRAY OF VEGETATION HEIGHTS 3104. C**** SPGSN - SPECIFIC GRAVITY OF SNOW 3105. C**** OUTPUT: 3106. C**** VH - VEGETATION HEIGHT 3107. C**** SNOWM - SNOW MASKING DEPTH 3108. C**** WFC1 - WATER FIELD CAPACITY OF TOP SOIL LAYER, M 3109. C**** 3110. PARAMETER (IM=72,JM=46) 3111. INCLUDE 'SOILS035.COM' 3112. C**** SOILS28 Common block 9/25/90 3113. REAL*8 SDATA 3114. COMMON/SOILS2/SDATA(IM,JM,11*NGM+1),ALA(3,IM,JM),ACS(3,IM,JM) 3115. * ,AFB(IM,JM),AFR(IM,JM,NGM),AVH(IM,JM) 3116. DIMENSION SHCAP(IMT) 3117. DATA SHCAP/2.E6,2.E6,2.E6,2.5E6,2.4E6/ 3118. ONE=1. 3119. ID=I0+100*J0 3120. C**** SET UP LAYERS 3121. DO 10 N=1,11*NGM+1 3122. DZ(N)=SDATA(I0,J0,N) 3123. 10 CONTINUE 3124. DO 20 N=1,NGM 3125. IF(DZ(N).LE.0.) GO TO 21 3126. 20 CONTINUE 3127. 21 N=N-1 3128. IF(N.LE.0) THEN 3129. WRITE (*,*) 'GHINIJ: N <= 0: I,J,N=',I0,J0,N,(DZ(K),K=1,43) 3130. STOP 3131. END IF 3132. C**** CALCULATE THE BOUNDARIES, BASED ON THE THICKNESSES. 3133. ZB(1)=0. 3134. DO 30 L=1,N 3135. 30 ZB(L+1)=ZB(L)-DZ(L) 3136. C**** CALCULATE THE LAYER CENTERS, BASED ON THE BOUNDARIES. 3137. DO 40 L=1,N 3138. 40 ZC(L)=.5*(ZB(L)+ZB(L+1)) 3139. C**** FR: ROOT FRACTION IN LAYER L (1=FR(1)+FR(2)+...+FR(N)) 3140. DO 45 L=1,N 3141. FR(L)=AFR(I0,J0,L) 3142. 45 CONTINUE 3143. C**** VH: VEGETATION HEIGHT 3144. VH=AVH(I0,J0) 3145. SNOWM=VH*SPGSN 3146. C**** FB,FV: BARE, VEGETATED FRACTION (1=FB+FV) 3147. FB=AFB(I0,J0) 3148. FV=1.-FB 3149. C**** ALAI: LEAF AREA INDEX 3150. ALAI=ALA(1,I0,J0)+COST*ALA(2,I0,J0)+SINT*ALA(3,I0,J0) 3151. ALAI=MAX(ALAI,ONE) 3152. ALAIC=5.0 3153. ALAIE=ALAIC*(1.-EXP(-ALAI/ALAIC)) 3154. C**** RS: MINIMUM STOMATAL RESISTANCE 3155. RS=ALAI/(ACS(1,I0,J0)+COST*ACS(2,I0,J0)+SINT*ACS(3,I0,J0)) 3156. C??? CNC=ALAI/RS REDEFINED BEFORE BEING USED (QSBAL,COND) 3157. C 3158. CW WRITE(6,*)'N=',N,' R=',R 3159. CW WRITE(6,91) 3160. CW 91 FORMAT(1X,5X,'ZB',5X,'ZC',5X,'DZ'/1X,21('-')) 3161. CW DO 95 L=1,N 3162. CW 95 WRITE(6,100)ZB(L),ZC(L),DZ(L) 3163. CW WRITE(6,100)ZB(N+1) 3164. CW100 FORMAT(1X,3F7.3) 3165. CW WRITE(6,*) 3166. C**** 3167. DO 60 IBV=1,2 3168. DO 60 L=1,N 3169. THETS(L,IBV)=0. 3170. THETM(L,IBV)=0. 3171. DO 50 I=1,IMT-1 3172. THETS(L,IBV)=THETS(L,IBV)+Q(I,L)*THM(0,I) 3173. THETM(L,IBV)=THETM(L,IBV)+Q(I,L)*THM(NTH,I) 3174. 50 CONTINUE 3175. WS(L,IBV)=THETS(L,IBV)*DZ(L) 3176. 60 CONTINUE 3177. WS(0,2)=.0001*ALAI 3178. WFC1=FB*WS(1,1)+FV*(WS(0,2)+WS(1,2)) 3179. C**** 3180. CALL XKLH0 3181. C**** 3182. DO 90 IBV=1,2 3183. DO 90 L=1,N 3184. SHC(L,IBV)=0. 3185. DO 80 I=1,IMT 3186. SHC(L,IBV)=SHC(L,IBV)+Q(I,L)*SHCAP(I) 3187. 80 CONTINUE 3188. SHC(L,IBV)=(1.-THETS(L,IBV))*SHC(L,IBV)*DZ(L) 3189. 90 CONTINUE 3190. C**** 3191. C SHC(0,2) IS THE HEAT CAPACITY OF THE CANOPY 3192. AA=ALA(1,I0,J0) 3193. SHC(0,2)=(.010+.002*AA+.001*AA**2)*SHW 3194. C**** 3195. C HTPR IS THE HEAT OF PRECIPITATION. 3196. C SHTPR IS THE SPECIFIC HEAT OF PRECIPITATION. 3197. SHTPR=0. 3198. IF(PR.GT.0.)SHTPR=HTPR/PR 3199. C HTPRS IS THE HEAT OF LARGE SCALE PRECIPITATION 3200. HTPRS=SHTPR*PRS 3201. C**** 3202. RETURN 3203. END 3301. BLOCK DATA 3302. C**** INITIALIZES COEFFICIENTS FOR SOIL FUNCTIONS 3303. C**** A ARE THE MATRIC POTENTIAL COEFFICIENTS 3304. C**** B ARE THE CONDUCTIVITY COEFFICIENTS 3305. INCLUDE 'SOILS035.COM' 3306. C**** SOILS28 Common block 9/25/90 3307. DATA A/ 3308. & 0.2514, 0.0136, -2.8319, 0.5958, 3309. & 0.1481, 1.8726, 0.1025, -3.6416, 3310. & 0.2484, 2.4842, 0.4583, -3.9470, 3311. & 0.8781, -5.1816, 13.2385,-11.9501/ 3312. DATA B/ 3313. & -0.4910, -9.8945, 9.7976, -3.2211, 3314. & -0.3238,-12.9013, 3.4247, 4.4929, 3315. & -0.5187,-13.4246, 2.8899, 5.0642, 3316. & -3.0848, 9.5497,-26.2868, 16.6930/ 3317. DATA P/ 3318. & -0.1800, -7.9999, 5.5685, -1.8868, 3319. & -0.1000,-10.0085, 3.6752, 1.2304, 3320. & -0.1951, -9.7055, 2.7418, 2.0054, 3321. & -2.1220, 5.9983,-16.9824, 8.7615/ 3322. END 3501. SUBROUTINE GHINHT (SNOWDP,TG1,TG2,WTR1,WTR2,ACE1,ACE2) 3502. C**** INITIALIZES NEW GROUND (W,HT,SNW) FROM OLD (T,W,ICE,SNW) 3503. C**** EVALUATES THE HEAT IN THE SOIL LAYERS BASED ON THE 3504. C**** TEMPERATURES. 3505. C**** INPUT: 3506. C**** W - WATER IN SOIL LAYERS, M 3507. C**** TP - TEMPERATURE OF LAYERS, C 3508. C**** FICE - FRACTION OF ICE OF LAYERS 3509. C**** FSN - HEAT OF FUSION OF WATER 3510. C**** SHC - SPECIFIC HEAT CAPACITY OF SOIL 3511. C**** SHI - SPECIFIC HEAT CAPACITY OF ICE 3512. C**** SHW - SPECIFIC HEAT CAPCITY OF WATER 3513. C**** SNOWD - SNOW DEPTH, EQUIVALENT WATER M 3514. C**** OUTPUT: 3515. C**** HT - HEAT IN SOIL LAYERS 3516. INCLUDE 'SOILS035.COM' 3517. C**** SOILS28 Common block 9/25/90 3518. C**** ADD CALCULATION OF WFC2 3519. C**** BASED ON COMBINATION OF LAYERS 2-N, AS IN RETP2 3520. WFC1=FB*WS(1,1)+FV*(WS(0,2)+WS(1,2)) 3521. WFC2=0. 3522. FBV=FB 3523. DO 30 IBV=1,2 3524. DO 20 L=2,N 3525. WFC2=WFC2+FBV*WS(L,IBV) 3526. 20 CONTINUE 3527. FBV=FV 3528. 30 CONTINUE 3529. WFC1=1000.*WFC1 3530. WFC2=1000.*WFC2 3531. FICE(0,2)=1. 3532. FICE(1,1)=(ACE1+SNOWDP)/(WTR1+ACE1+SNOWDP+1.E-20) 3533. FICE(1,2)=FICE(1,1) 3534. TP(0,2)=TG1 3535. C**** W = SNOW(IF TOP LAYER) + WMIN + (WMAX-WMIN)*(WTR+ICE)/WFC 3536. W(0,2)=0. 3537. DO 1 IBV=1,2 3538. W(1,IBV)=SNOWDP/1000. 3539. WMIN=THETM(1,IBV)*DZ(1) 3540. WET1=(WTR1+ACE1)/WFC1 3541. IF(WET1.GT.1.) WET1=1. 3542. W(1,IBV)=W(1,IBV)+WMIN+(WS(1,IBV)-WMIN)*WET1 3543. SNOWD(IBV)=SNOWDP/1000. 3544. TP(1,IBV)=TG1 3545. DO 1 L=2,N 3546. FICE(L,IBV)=ACE2/(WTR2+ACE2+1.E-20) 3547. WMIN=THETM(L,IBV)*DZ(L) 3548. WET2=(WTR2+ACE2)/WFC2 3549. IF(WET2.GT.1.) WET2=1. 3550. W(L,IBV)=WMIN+(WS(L,IBV)-WMIN)*WET2 3551. TP(L,IBV)=TG2 3552. 1 CONTINUE 3553. C**** 3554. ENTRY GHEXHT 3555. C**** 3556. C**** COMPUTE HT (HEAT W/M+2) 3557. DO 10 IBV=1,2 3558. LL=2-IBV 3559. DO 10 L=LL,N 3560. IF(TP(L,IBV))2,4,6 3561. 2 HT(L,IBV)=TP(L,IBV)*(SHC(L,IBV)+W(L,IBV)*SHI)-W(L,IBV)*FSN 3562. GO TO 10 3563. 4 HT(L,IBV)=-FICE(L,IBV)*W(L,IBV)*FSN 3564. GO TO 10 3565. 6 HT(L,IBV)=TP(L,IBV)*(SHC(L,IBV)+W(L,IBV)*SHW) 3566. 10 CONTINUE 3567. IF(ID.EQ.0)THEN 3568. WRITE(*,*)'GHINHT ID CHECK',ID 3569. WRITE(*,*)'TG1,TG2',TG1,TG2 3570. WRITE(*,*)'TP',TP 3571. WRITE(*,*)'HT',HT 3572. WRITE(*,*)'W',W 3573. WRITE(*,*)'WTR1,WTR2',WTR1,WTR2 3574. WRITE(*,*)'ACE1,ACE2',ACE1,ACE2 3575. WRITE(*,*)'WFC1,WFC2',WFC1,WFC2 3576. WRITE(*,*)'SHC',SHC 3577. WRITE(*,*)'FICE',FICE 3578. ENDIF 3579. RETURN 3580. END 3701. SUBROUTINE RETP2 (TG2AV,WTR2AV,ACE2AV) 3702. C**** EVALUATES THE MEAN TEMPERATURE IN THE SOIL LAYERS 2-NGM 3703. C**** AS WELL AS THE WATER AND ICE CONTENT. 3704. C**** INPUT: 3705. C**** W - WATER IN SOIL LAYERS, M 3706. C**** HT - HEAT IN SOIL LAYERS 3707. C**** FSN - HEAT OF FUSION OF WATER 3708. C**** SHC - SPECIFIC HEAT CAPACITY OF SOIL 3709. C**** SHI - SPECIFIC HEAT CAPACITY OF ICE 3710. C**** SHW - SPECIFIC HEAT CAPCITY OF WATER 3711. C**** OUTPUT: 3712. C**** TG2AV - TEMPERATURE OF LAYERS 2 TO NGM, C 3713. C**** ICE2AV - ICE AMOUNT IN LAYERS 2 TO NGM, KG/M+2 3714. C**** WTR2AV - WATER IN LAYERS 2 TO NGM, KG/M+2 3715. INCLUDE 'SOILS035.COM' 3716. C**** SOILS28 Common block 9/25/90 3717. TG2AV=0. 3718. WTR2AV=0. 3719. ACE2AV=0. 3720. DO 3500 IBV=1,2 3721. WC=0. 3722. HTC=0. 3723. SHCC=0. 3724. DO 3420 L=2,N 3725. WC=WC+W(L,IBV) 3726. HTC=HTC+HT(L,IBV) 3727. 3420 SHCC=SHCC+SHC(L,IBV) 3728. TPC=0. 3729. FICEC=0. 3730. IF(WC.NE.0.) FICEC=-HTC/(FSN*WC) 3731. IF(FSN*WC+HTC.GE.0.)GO TO 3430 3732. TPC=(HTC+WC*FSN)/(SHCC+WC*SHI) 3733. FICEC=1. 3734. GO TO 3440 3735. 3430 IF(HTC.LE.0.) GO TO 3440 3736. TPC=HTC/(SHCC+WC*SHW) 3737. FICEC=0. 3738. 3440 CONTINUE 3739. FTP=FB 3740. IF(IBV.EQ.2) FTP=FV 3741. TG2AV=TG2AV+TPC*FTP 3742. WTR2AV=WTR2AV+WC*FTP*1000.*(1.-FICEC) 3743. ACE2AV=ACE2AV+WC*FTP*1000.*FICEC 3744. 3500 CONTINUE 3745. RETURN 3746. END 3801. C SUBROUTINE NEEDED LATER FOR EXTRA DIAGNOSTICS 3802. SUBROUTINE OUTW(I) 3803. C**** PRINTS THETA VALUES AT TIME STEP I 3804. INCLUDE 'SOILS035.COM' 3805. C**** SOILS28 Common block 9/25/90 3806. CALL WTAB 3807. ICHN=6 3808. SCNDS=I*DT 3809. DAY=INT(SCNDS/86400.) 3810. IDAY=DAY 3811. HOUR=INT((SCNDS-86400.*DAY)/86400.) 3812. IHOUR=HOUR 3813. C PRINT 1001 3814. WRITE(ICHN,1000) 3815. WRITE(ICHN,*)'GENERAL QUANTITIES (BARE SOIL OR VEGETATION)' 3816. WRITE(ICHN,*)'ID,DTS',ID,DTS 3817. CC WRITE(ICHN,1021) 3818. WRITE(ICHN,1045) 3819. WRITE(ICHN,1023)'DAY= ',IDAY,'PR= ',PR,'TS= ',TS-TFRZ, 3820. * 'Q1= ',Q1 3821. WRITE(ICHN,1023)'HOUR= ',IHOUR,'SNOWF= ',SNOWF,'TG= ',TG-TFRZ, 3822. * 'QS= ',QS 3823. WRITE(ICHN,1043)'T1= ',T1-TFRZ,'VG= ',VG,'CH= ',CH 3824. WRITE(ICHN,1044)'VSM= ',VSM 3825. WRITE(ICHN,1022) 3826. WRITE(ICHN,1021) 3827. WRITE(ICHN,1014)'BARE SOIL FB = ',FB 3828. 1014 FORMAT(1X,A17,F4.2) 3829. CC WRITE(ICHN,1021) 3830. WRITE(ICHN,1025) 3831. WRITE(ICHN,1026) 3832. WRITE(ICHN,1027) SNOWD(1),RNF(1),EVAP(1),XINFC(1),ZW(1),QB 3833. WRITE(ICHN,1021) 3834. WRITE(ICHN,1030) 3835. WRITE(ICHN,1031) 3836. DO 100 L=1,N 3837. WRITE(ICHN,1040)L,THETA(L,1),TP(L,1),FICE(L,1),RNFF(L,1),F(L,1), 3838. & H(L,1),XK(L,1),W(L,1),WS(L,1),SHC(L,1),FH(L,1),HT(L,1), 3839. & Q(1,L),Q(2,L),Q(3,L),Q(4,L) 3840. 100 CONTINUE 3841. CC WRITE(ICHN,1021) 3842. WRITE(ICHN,1022) 3843. WRITE(ICHN,1021) 3844. WRITE(ICHN,1014)'VEGETATION FV = ',FV 3845. CC WRITE(ICHN,1021) 3846. WRITE(ICHN,1035) 3847. WRITE(ICHN,1036) 3848. WRITE(ICHN,1037) SNOWD(2),RNF(2),EVAP(2),XINFC(2),ZW(2),QC,EVAPW, 3849. * EVAPD,DR,FW 3850. WRITE(ICHN,1021) 3851. WRITE(ICHN,1030) 3852. WRITE(ICHN,1031) 3853. L=0 3854. WRITE(ICHN,1049)L,THETA(L,2),TP(L,2),FICE(L,2),RNFF(L,2),F(L,2), 3855. & W(L,2),WS(L,2),SHC(L,2),FH(L,2),HT(L,2) 3856. DO 200 L=1,N 3857. WRITE(ICHN,1040)L,THETA(L,2),TP(L,2),FICE(L,2),RNFF(L,2),F(L,2), 3858. & H(L,2),XK(L,2),W(L,2),WS(L,2),SHC(L,2),FH(L,2),HT(L,2), 3859. & Q(1,L),Q(2,L),Q(3,L),Q(4,L) 3860. 200 CONTINUE 3861. CC WRITE(ICHN,1021) 3862. WRITE(ICHN,1022) 3863. WRITE(ICHN,1021) 3864. WRITE(ICHN,1055) 3865. 1055 FORMAT(1X,3X,9X,' KGM-2',2X,9X,' KGM-2',3X,9X,' KGM-2', 3866. * 4X,9X,'1E6JM-2',3X,9X,'1E6JM-2') 3867. WRITE(ICHN,1060) ARUNS,AEVAPW,AEPC,AFHG,AF0DT 3868. 1060 FORMAT(1X,3X,'ARUNS = ',F9.4,2X,'AEVAPW = ',0PF9.4,4X,'AEPC = ', 3869. * 0PF9.4,4X,'AFHG = ',-6PF9.4,2X,'AF0DT = ',-6PF9.4) 3870. WRITE(ICHN,1065) ARUNU,AEVAPD,AEPB,ATRG,HTPR*DTS 3871. 1065 FORMAT(1X,3X,'ARUNU = ',F9.4,2X,'AEVAPD = ',0PF9.4,4X,'AEPB = ', 3872. * 0PF9.4,4X,'ATRG = ',-6PF9.4,2X,'APHDT = ',-6PF9.4) 3873. WRITE(ICHN,1070) AEVAPB,ASHG,AERUNS 3874. 1070 FORMAT(1X,3X,' ',9X,2X,'AEVAPB = ',0PF9.4,4X,' ', 3875. * 9X,4X,'ASHG = ',-6PF9.4,2X,'AERNS = ',-6PF9.4) 3876. WRITE(ICHN,1073) ALHG,AF1DT,AEDIFS 3877. 1073 FORMAT(1X,3X,' ',9X,2X,' ',9X,4X,' ', 3878. * 9X,4X,'ALHG = ',-6PF9.4,2X,'AF1DT = ',-6PF9.4/ 3879. * 1X,3X,' ',9X,2X,' ',9X,4X,' ', 3880. * 9X,4X,' ',2X,9X,'AEDFS = ',-6PF9.4) 3881. RETURN 3882. 1000 FORMAT(1H ,121('=')) 3883. 1001 FORMAT(1H1) 3884. 1010 FORMAT(1X,A20,F10.0) 3885. 1020 FORMAT(1X,A20,1PE12.4) 3886. 1021 FORMAT(1H0) 3887. 1022 FORMAT(1H ,60('. '),'.') 3888. 1023 FORMAT(1X,A10,I10,A10,6PF10.2,2(A10,0PF8.2),A10,0PF8.4) 3889. 1043 FORMAT(1X,10X,10X,10X,10X,2(A10,0PF8.2),A10,0PF8.4) 3890. 1044 FORMAT(1X,10X,10X,38X,1(A10,F8.2)) 3891. 1045 FORMAT(1X,20X,10X,' 1E-6MS-1',10X,4X,'T(C)',10X,4X,'MS-1') 3892. 1024 FORMAT(1X,6(A10,E10.2)) 3893. 1015 FORMAT(1X,4(A8,F8.2)) 3894. 1019 FORMAT(1X,12X,4(A8,F8.2)) 3895. 1025 FORMAT(1H ,5X,'SNOWD',7X,'RNF',6X,'EVAP',6X,'XINFC', 3896. * 8X,'ZW',8X,'QB') 3897. 1026 FORMAT(1H ,' MH2O',2X,'1E-6MS-1',2X,'1E-6MS-1',3X, 3898. * '1E-6MS-1',3X,' M',5X,' ') 3899. 1027 FORMAT(1H ,0PF10.4,6PF10.4,6PF10.4,1X,6PF10.1,0PF10.4, 3900. * 0PF10.4) 3901. 1030 FORMAT(1H ,5X,'THETA',3X,'TP',2X,'FICE',4X,'RUNOFF' 3902. & ,8X,'FL',9X,'H',8X,'XK',6X,'W',5X,'WS',8X, 3903. & 'SHC',8X,'FH',8X,'HT',1X,'SAND',1X,'LOAM',1X,'CLAY',1X,'PEAT') 3904. 1031 FORMAT(1H ,5X,5X,2X,'(C)',2X,6X,'1E-6MS-1',2X,'1E-6MS-1', 3905. & 3X,' M',2X,'1E-6MS-1',1X,' M',1X,' M', 3906. & 1X,'1E6JM-3C-1',4X,' WM-2',3X,'1E6JM-2',4X,'%',4X,'%',4X,'%' 3907. & ,4X,'%'/1X,125('-')) 3908. 1035 FORMAT(1H ,5X,'SNOWD',7X,'RNF',6X,'EVAP',6X,'XINFC', 3909. * 8X,'ZW',8X,'QC',5X,'EVAPW',5X,'EVAPD',8X,'DR',8X,'FW') 3910. 1036 FORMAT(1H ,' MH2O',2X,'1E-6MS-1',2X,'1E-6MS-1',3X, 3911. * '1E-6MS-1',2X,' M',5X,' ',2X,'1E-6MS-1',2X, 3912. * '1E-6MS-1',2X,'1E-6MS-1',' ') 3913. 1037 FORMAT(1H ,0PF10.4,6PF10.4,6PF10.4,1X,6PF10.1,0PF10.4, 3914. * 0PF10.4,6PF10.4,6PF10.4,6PF10.4,0PF10.2) 3915. 1040 FORMAT(1X,I3,F7.3,F5.1,F6.3,1P,6PF10.4,6PF10.4,0PF10.3,6PF10.4, 3916. * 0PF7.4,0PF7.4,1X,-6PF10.4,0PF10.4,-6PF10.4,4(2PF5.1)) 3917. 1049 FORMAT(1X,I3,F7.3,F5.1,F6.3,1P,6PF10.4,6PF10.4,10X,10X, 3918. * 0PF7.4,0PF7.4,1X,-6PF10.4,0PF10.4,-6PF10.4,3(2PF5.1)) 3919. END 3920. C**** 4001. SUBROUTINE WTAB 4002. C**** RETURNS WATER TABLE ZW FOR IBV=1 AND 2. 4003. C**** INPUT: 4004. C**** ZB - LAYER BOUNDARIES, M 4005. C**** ZC - SOIL CENTERS, M 4006. C**** DZ - LAYER THICKNESSES, M 4007. C**** H - SOIL POTENTIAL OF LAYERS, M 4008. C**** F - FLUXES BETWEEN LAYERS, M S-1 4009. C**** XK - CONDUCTIVITIES OF LAYERS, M S-1 4010. C**** OUTPUT: 4011. C**** ZW(2) - WATER TABLE FOR IBV=1 AND 2, M 4012. INCLUDE 'SOILS035.COM' 4013. C**** SOILS28 Common block 9/25/90 4014. TOL=1.E-6 4015. DO 100 IBV=1,2 4016. C**** FIND NON-SATURATED LAYER 4017. DO 10 L=N,1,-1 4018. IF(W(L,IBV).LT.WS(L,IBV)*(1.-TOL))GO TO 20 4019. 10 CONTINUE 4020. L=1 4021. 20 CONTINUE 4022. C**** RETRIEVE MATRIC POTENTIAL 4023. C WRITE(6,*) 'ID,N,L,HMAT,IBV,XK(L,IBV)',ID,N,L,HMAT,IBV,XK(L,IBV) 4024. HMAT=H(L,IBV)-ZC(L) 4025. C**** CALCULATE DENOMINATOR, AND KEEP ZW ABOVE ZB(L+1) 4026. IF(XK(L,IBV).LE.1.E-20) THEN 4027. DENOM=-2.*HMAT/DZ(L) 4028. GO TO 90 4029. END IF 4030. DENOM=MAX(F(L,IBV)/XK(L,IBV)+1.,-2.*HMAT/DZ(L)) 4031. 90 CONTINUE 4032. C**** CALCULATE WATER TABLE 4033. C WRITE(6,*) 'DENOM',DENOM 4034. ZW(IBV)=ZB(L)-SQRT(-2.*HMAT*DZ(L)/(DENOM+1.E-20)) 4035. 100 CONTINUE 4036. RETURN 4037. END 4101. SUBROUTINE OUTGH 4102. C**** 4103. C**** CALLED AT END OF ADVNC. 4104. C**** USER WRITTEN SUBROUTINE. 4105. C**** CAN BE LOADED BEFORE SOILS MODULE TO BE USED INSTEAD OF 4106. C**** SOILS MODEL VERSION, WHICH IS BLANK. 4107. C**** 4108. C**** VALUES OF SOIL QUANTITIES CAN BE STORED AFTER EACH TIMESTEP 4109. C**** GCM TAU CAN BE EXAMINED, AS WELL AS GCM ID TO DETERMINE 4110. C**** WHETHER TO PRINT OUT. 4111. RETURN 4112. END 4201. SUBROUTINE GHIJ (I0,J0,TGGR,WEARTH) 4202. C**** 4203. C**** MAKES CALLS TO ADVANCE GROUND HYDROLOGY BY ONE TIME STEP 4204. C**** AT A GIVEN GRIDBOX. 4205. C**** 4206. C**** INPUT: 4207. C**** I0,J0 - THE GRIDBOX 4208. C**** 4209. C**** OUTPUT: 4210. C**** TGGR - THE LAND SURFACE TEMPERATURE, C 4211. C**** WEARTH - RELATIVE WETNESS OF FIRST LAYER 4212. C**** 4213. INCLUDE 'SOILS035.COM' 4214. REAL*8 AOUT(26) 4215. EQUIVALENCE (AOUT,ACNA) 4216. C**** Zero out output variables that will be summed over time steps 4217. DO 10 K=1,26 4218. 10 AOUT(K) = 0. 4219. C**** 4220. CALL GHINIJ (I0,J0,WFC1) 4223. CALL ADVNC 4226. TGGR = TBCS 4227. WBARE = 1. 4227.1 WVEG = 1. 4227.2 IF(FB.gt.0. .and. SNOWD(1).le.0.) WBARE = W(1,1)/WS(1,1) 4227.3 IF(FV.gt.0. .and. SNOWD(2).le.0.) WVEG = W(1,2)/WS(1,2) 4227.4 WEARTH = FB*WBARE + FV*WVEG 4228. RETURN 4229. END