/OPENFVS_DIR_RESTRUCT_TEST_BRANCH/src/variants/ca/src/htcalc.f

http://open-fvs.googlecode.com/
FORTRAN Legacy | 187 lines | 93 code | 0 blank | 94 comment | 0 complexity | 8057117b117f4cc5714ac23ff6d80812 MD5 | raw file
  1. SUBROUTINE HTCALC(JFOR,SINDX,KSPEC,AG,HGUESS,JOSTND,DEBUG)
  2. IMPLICIT NONE
  3. C----------
  4. C **HTCALC--CA DATE OF LAST REVISION: 01/10/11
  5. C----------
  6. C THIS ROUTINE CALCULATES A POTENTIAL HT GIVEN A SITE INDEX, SPECIES,
  7. C AND AGE. IT IS USED TO CALCULATE POTENTIAL HEIGHT GROWTH AND SITE
  8. C EQUIVALENTS.
  9. C----------
  10. COMMONS
  11. C
  12. C
  13. INCLUDE 'PRGPRM.F77'
  14. C
  15. C
  16. INCLUDE 'VARCOM.F77'
  17. C
  18. C
  19. COMMONS
  20. C----------
  21. LOGICAL DEBUG
  22. INTEGER INDEX(49),INDX,KSPEC,JFOR,JOSTND
  23. REAL HGUESS,AG,SINDX,B0,B1,B2,B3,B4,B5,B6,Z,A,B,C
  24. REAL X1,X2,TOPTRM,BOTTRM,TERM,TERM2,B50
  25. REAL DUNL1(6),DUNL2(6),DUNL3(6)
  26. C----------
  27. C SPECIES ORDER IN CA VARIANT:
  28. C 1=PC 2=IC 3=RC 4=WF 5=RF 6=SH 7=DF 8=WH 9=MH 10=WB
  29. C 11=KP 12=LP 13=CP 14=LM 15=JP 16=SP 17=WP 18=PP 19=MP 20=GP
  30. C 21=JU 22=BR 23=GS 24=PY 25=OS 26=LO 27=CY 28=BL 29=EO 30=WO
  31. C 31=BO 32=VO 33=IO 34=BM 35=BU 36=RA 37=MA 38=GC 39=DG 40=FL
  32. C 41=WN 42=TO 43=SY 44=AS 45=CW 46=WI 47=CN 48=CL 49=OH
  33. C----------
  34. DATA INDEX/ 3, 3, 3, 3, 5, 5, 3, 3, 5, 6,
  35. & 6, 6, 6, 6, 4, 3, 4, 4, 4, 4,
  36. & 6, 3, 3, 7, 3, 9, 9, 9, 9, 7,
  37. & 7, 7, 9, 10, 7, 10, 9, 9, 7, 7,
  38. & 10, 8, 10, 10, 10, 10, 10, 10, 10/
  39. C
  40. DATA DUNL1/ -88.9, -82.2, -78.3, -82.1, -56.0, -33.8 /
  41. DATA DUNL2/ 49.7067, 44.1147, 39.1441,
  42. & 35.4160, 26.7173, 18.6400 /
  43. DATA DUNL3/ 2.375, 2.025, 1.650, 1.225, 1.075, 0.875 /
  44. C
  45. IF(DEBUG)WRITE(JOSTND,15)
  46. 15 FORMAT(' ENTERING HTCALC')
  47. IF(DEBUG)WRITE(JOSTND,30)AG,SINDX,KSPEC
  48. 30 FORMAT(' IN HTCALC 30F AG,SINDX,KSPEC=',2F10.4,I10)
  49. C
  50. SELECT CASE (JFOR)
  51. C----------
  52. C REGION 5 LOGIC
  53. C R5 USE DUNNING/LEVITAN SITE CURVE.
  54. C----------
  55. CASE(1:5)
  56. IF(JFOR .LE. 5) THEN
  57. C----------
  58. C SET UP MAPPING TO THE CORRECT DUNNING-LEVITAN SITE CURVE
  59. C----------
  60. IF(SINDX .LE. 44.) THEN
  61. INDX=6
  62. ELSEIF (SINDX.GT.44. .AND. SINDX.LE.52.) THEN
  63. INDX=5
  64. ELSEIF (SINDX.GT.52. .AND. SINDX.LE.65.) THEN
  65. INDX=4
  66. ELSEIF (SINDX.GT.65. .AND. SINDX.LE.82.) THEN
  67. INDX=3
  68. ELSEIF (SINDX.GT.82. .AND. SINDX.LE.98.) THEN
  69. INDX=2
  70. ELSE
  71. INDX=1
  72. ENDIF
  73. IF(AG .LE. 40.) THEN
  74. HGUESS = DUNL3(INDX) * AG
  75. ELSE
  76. HGUESS = DUNL1(INDX) + DUNL2(INDX)*ALOG(AG)
  77. ENDIF
  78. ENDIF
  79. C
  80. CASE DEFAULT
  81. C----------
  82. C REGION 6 LOGIC
  83. C
  84. C
  85. C LOAD COEFFICIENTS AND GO TO A SITE CURVE BASED ON SPECIES.
  86. C COEFFS IN BB_ ARRAYS ARE LOADED IN BLKDAT:
  87. C 1 = KING DOUGLAS-FIR; WEYERH FOR PAP 8
  88. C 2 = DOLPH WHITE FIR PSW 185;
  89. C 3 = HANN-SCRIVANI DOUGLAS-FIR
  90. C 4 = HANN-SCRIVANI PONDEROSA PINE
  91. C 5 = DOLPH RED FIR PSW 206
  92. C 6 = DAHMS LODGEPOLE RP-PNW-8
  93. C 7 = POWERS BLACK OAK RES NOTE PSW-262
  94. C 8 = PORTER & WIANT TANOAK JOF 4/95 286-287
  95. C 9 = PORTER & WIANT MADRONE JOF 4/95 286-287
  96. C 10 = PORTER & WIANT RED ALDER JOF 4/95 286-287
  97. C
  98. C INDEX = ARRAY HOLDING THE INDEX FOR ACCESSING THE APPROPRIATE SITE
  99. C CURVE COEFFICIENTS IN THE BB_ ARRAYS.
  100. C----------
  101. INDX = INDEX(KSPEC)
  102. B0=BB0(INDX)
  103. B1=BB1(INDX)
  104. B2=BB2(INDX)
  105. B3=BB3(INDX)
  106. B4=BB4(INDX)
  107. B5=BB5(INDX)
  108. B6=BB6(INDX)
  109. IF(DEBUG) WRITE(JOSTND,*)' INDX,B0,B1,B2,B3,B4,B5,B6= ',
  110. & INDX,B0,B1,B2,B3,B4,B5,B6
  111. C----------
  112. C GO TO A DIFFERENT POTENTIAL HEIGHT CURVE DEPENDING ON THE INDEX
  113. C SCALE FACTORS ADDED TO SOME SPECIES SO THEY WOULD
  114. C HIT THE SITE HEIGHT.
  115. C----------
  116. SELECT CASE (INDX)
  117. C----------
  118. C KING DOUGLAS-FIR, WEYERHAUSER FOR PAP 8, 1966.
  119. C----------
  120. CASE(1)
  121. Z = B0/(SINDX - 4.5)
  122. A = B1 + B2*Z
  123. B = B3 + B4*Z
  124. C = B5 + B6*Z
  125. HGUESS = ((AG*AG)/(A + B*AG + C*AG*AG)) + 4.5
  126. C----------
  127. C DOLPH WHITE FIR, RES PAP PSW-185, FEB 1987
  128. C----------
  129. CASE(2)
  130. X1 = (B1*AG**B2*EXP(B3*AG))
  131. X2 = B4*(1.0 - EXP(B5*AG**B6))
  132. HGUESS = (SINDX - B0 + X1*X2)/X1 + 4.5
  133. IF(DEBUG)WRITE(JOSTND,*)' X1,X2,SINDX,AG,HGUESS= ',
  134. & X1,X2,SINDX,AG,HGUESS
  135. C----------
  136. C HANN-SCRIVANI OSU RES BULL 59, FEB 1987
  137. C----------
  138. CASE(3,4)
  139. TOPTRM = 1.0-EXP(-EXP(B0+B1*ALOG(SINDX-4.5)+B2*ALOG(AG)))
  140. BOTTRM = 1.0-EXP(-EXP(B0+B1*ALOG(SINDX-4.5)+B2*ALOG(50.)))
  141. HGUESS = ((SINDX-4.5)*TOPTRM/BOTTRM) + 4.5
  142. HGUESS = HGUESS*1.05
  143. IF(DEBUG)WRITE(JOSTND,*)' SINDX,AG,TOPTRM,BOTTRM,HGUESS= ',
  144. & SINDX,AG,TOPTRM,BOTTRM,HGUESS
  145. C----------
  146. C DOLPH RED FIR CURVES, RES PAP PSW 206
  147. C----------
  148. CASE(5)
  149. TERM=AG*EXP(AG*B3)*B2
  150. B = SINDX*TERM + B4*TERM*TERM + B5
  151. TERM2 = 50.0 * EXP(50.0*B3) * B2
  152. B50 = SINDX*TERM2 + B4*TERM2*TERM2 + B5
  153. HGUESS = ((SINDX-4.5) * (1.0-EXP(-B*(AG**B1)))) /
  154. & (1.0-EXP(-B50*(50.0**B1)))
  155. HGUESS = HGUESS + 4.5
  156. IF(DEBUG)WRITE(JOSTND,*)' SINDX,AG,TERM,B,TERM2,B50,HGUESS= ',
  157. & SINDX,AG,TERM,B,TERM2,B50,HGUESS
  158. C----------
  159. C DAHMS LODGEPOLE PINE, RES PAP PNW-8
  160. C----------
  161. CASE(6)
  162. HGUESS = SINDX * (B0 + B1*AG + B2*AG*AG)
  163. HGUESS = HGUESS*1.10
  164. IF(DEBUG)WRITE(JOSTND,*)' SINDX,AG,HGUESS= ',SINDX,AG,HGUESS
  165. C----------
  166. C POWERS BLACK OAK RES NOTE PSW-262
  167. C----------
  168. CASE(7)
  169. TERM = SQRT(AG)-SQRT(50.)
  170. HGUESS = (SINDX * (1 + B1*TERM)) - B0*TERM
  171. HGUESS = HGUESS * 0.70
  172. IF(DEBUG)WRITE(JOSTND,*)' SINDX,AG,TERM,HGUESS= ',
  173. & SINDX,AG,TERM,HGUESS
  174. C----------
  175. C PORTER & WIANT TANOAK, MADRONE, RED ALDER, JOF 4/1965 286-287p
  176. C----------
  177. CASE(8:10)
  178. HGUESS = SINDX / (B0 + B1/AG)
  179. HGUESS = HGUESS * 0.80
  180. IF(DEBUG)WRITE(JOSTND,*)' SINDX,AG,HGUESS= ',SINDX,AG,HGUESS
  181. C
  182. END SELECT
  183. C
  184. END SELECT
  185. C
  186. RETURN
  187. END