home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-387-Vol-3of3.iso / h / house_ii.zip / FOR / GROUND.FOR < prev    next >
Text File  |  1992-05-09  |  30KB  |  1,007 lines
  1.       SUBROUTINE GRNDEX(TODDB,WSPED,SOLRH,KM) 
  2. C 
  3. C TAPE60 WRITE.  ECHO OF NAMELIST INPUT 
  4. C TAPE18 READ.   NAMELIST INPUT 
  5. C TAPE41 READ.   SEQUENTIAL COMPRESSED WEATHER DATA  (WORKING FILE) 
  6. C TAPE43 WRITE/READ. GROUND TEMPERATURES(INTERNAL TO THIS ROUTINE)
  7. C 
  8. C  EXECUTIVE CODE FOR DETAILED GROUND HEAT TRANSFER 
  9. C   BY R D FISCHER, BATTELLE, WITH SUGGESTIONS FOR APPROACH FROM
  10. C   TOM BECKEY OF HONEYWELL 
  11. C  ICODE = 1  READ NAMELIST INPGND, INITIALIZE NODAL GEOM 
  12. C             AND HEAT TRANSFER COEFS (ICODE=2 IS 
  13. C             INCLUDED IN ICODE = 1 FOR NOW)
  14. C        =2   RUN PRECONDITIONING MONTHS
  15. C        =3   HOURLY GROUND HEAT TRANSFER CALC. 
  16. C             (USES MEAN HRLY VALUES FOR TODDB,WSPED,SOLRH) 
  17. C        =4   FINISH GROUND CALCULATIONS FOR THE MONTH
  18. C             FOR THE MONTH  INT(DAYMON/(NDAYS-1) TIMES TO
  19. C                 INITIALIZE FOR NEXT MONTH.
  20. C  TODDB  MEAN HRLY OUTDOOR DRYBULB AIR TEMP, F 
  21. C  WSPED  MEAN HRLY WIND SPEED,MPH
  22. C  SOLRH  MEAN HRLY TOTAL SOLAR RADIATION ON HORIZONTAL SURFACE 
  23. C  TIMEYR   TIME OF YEAR AT START OF GROUND PRECONDITIONING 
  24. C           (SET IN WWEATH) 
  25. C  NMNPRE  NO. OF MONTHS OF PRECONDITIONING(READ IN WWEATH) 
  26. C  NDHRS  NO. OF HOURS PER MONTH(INCL PRECONDITIONING DAY-- 
  27. C            SET IN WWEATH) 
  28. C  IGINIT =1 TO USE TSOILI FOR INITIAL GROUND TEMPS 
  29. C         =2 TO USE TG VALUES FROM TAPE43 FOR INITIALIZATION
  30. C        =-2 WRITE INITIAL TG'S AT START OF FIRST PRECOND DAY TO TAPE43 
  31. CMDK NKONST
  32. CMDK NWL
  33. CMDK NWN
  34. CMDK NZN
  35. CMDK NZW
  36. CMDK BLK53
  37. CMDK BLK55
  38. CMDK BLKBSF
  39. CMDK BLKHWG
  40. CMDK BLKQS
  41. CMDK CNSTRK
  42. CMDK ENCBK1
  43. CMDK ENCBK2
  44. CMDK ENCBLK
  45. CMDK IZWQ
  46. CMDK SOILB
  47. CMDK SURFAR
  48. CMDK TEMPB
  49. CMDK TYRGND
  50. CMDK UABAS
  51.       DIMENSION ODDB(24,4),WSPEED(24,4),SOLARH(24,4), 
  52.      +TDBT(24), WSPDT(24), SOLRT(24),TPR(28),IP(28) 
  53.      +,NDYM(12)
  54.       INTEGER*4 LAST1
  55.       NAMELIST/INPGND/ABSRPG,CPSOIL,IGINIT,NPRINT,PHASE,
  56.      +                RFIELD,RHSOIL,TEAMP,TEAV,TSOILI,
  57.      +                UBASF,UBWBG,XKSOIL,ZDEPTH
  58. C INPUTS IN NAMELIST INPGND:
  59. C  ABSRPG - ABSORPTIVITY OF GROUND TO SOLAR RADIATION
  60. C  CPSOIL - SPECIFIC HEAT OF SOIL AROUND BASEMENT, BTU/LBM-F
  61. C  IGINIT - =1 USE TSOILI FOR INITAIL GROUND TEMPERATURES
  62. C           =2 USE GROUND TEMP FROM TAPE43 FOR INITIALIZATION
  63. C              (SET IN A PREV RUN WITH IGINIT = -1 -- THIS OPTION
  64. C              WILL REDUCE EXECUTION TIME FOR HANDBOOK RUNS)
  65. C           = -2 INITIAL SOIL TEMPS WILL BE CALC FROM A KUSUDA
  66. C             RELATIONSHIP AND WRITTEN TO TAPE43
  67. C    NOTE: YOU MUST RUN ONCE WITH -2 AND THEN CHANGE TO 2!!!
  68. C  NPRINT - PRINT GROUND TEMP EVERY NPRINT HOURS
  69. C  PHASE  - PHASE ANGLE OF EARTH TEMP AT SURFACE(SEE KUSUDA REF.),
  70. C           RADIANS
  71. C  RFIELD - RADIUS FROM CENTER OF BASEMENT TO FAR-FIELD NODE, FT      
  72. C  RHSOIL - DENSITY OF SOIL, LBM/FT3
  73. C  TEAMP  - TEMPERATURE AMPLITUDE AT EARTHS SURFACE, F
  74. C  TEAV   - AVERAGE ANNUAL EARTH TEMPERATURE, F
  75. C  TSOILI - INITIAL VALUE FOR SOIL TEMPERATURE, F
  76. C  UBASF AND UBWBG  -  OVERALL HEAT FLOW TERMS AS USED IN
  77. C         THE SINGLE ZONE MODEL AND USED HERE FOR PRECONDITIONING
  78. C         THE GROUND NEAR THE BSMT WALLS AND FLOOR,(USE 0.5 FOR BOTH)     
  79. C         THESE VALUES ARE USED IN BSMTPRE.FOR
  80. C  XKSOIL - THERMAL CONDUCTIVITY OF SOIL, BTU/HR-FT-F
  81. C  ZDEPTH - DEPTH FROM GROUND SURFACE TO CONSTANT TEMP NODE, FT
  82.       DATA IP/1,3,5,6,7,2,4,8,10,12,13,14,
  83.      +        9,11,15,27,17,18,16,28,19,20, 
  84.      +        21,22,23,24,25,26/
  85.       DATA NDAYQ/1/,NDBG/0/ ,IPRBG/1/
  86.       DATA NDYM/31,28,31,30,31,30,31,31,30,31,30,31/ 
  87.       GO TO (100,200,300,400),ICODE 
  88.   100 CONTINUE
  89. C   
  90.       OPEN(18,FILE='TAPE18',STATUS='OLD',IOSTAT=IO18)
  91.       IF(IO18.NE.0)THEN
  92.         WRITE(60,*) ' GROUND: CANT OPEN TAPE18 WITH NAMELIST INPGND'
  93.         STOP ' GROUND: CANT OPEN TAPE18 WITH NAMELIST INPGND'
  94.         END IF
  95.       READ(18,INPGND,END=999)
  96.       CLOSE(18)
  97.       WRITE(60,INPGND)
  98. C     
  99. C  FIND BSMT FLOOR AND WALL DATA FROM ENCLOSURE DATA
  100.       KF=KONSTA(IZWD(1))
  101.       KW=KONSTA(IZWE(1))
  102.       NLMF=NLMP(KF) 
  103.       NLMW=NLMP(KW) 
  104.       THFLR=0.
  105.       THWAL=0.
  106.       DO 102 I=1,NLMF 
  107.   102 THFLR=THFLR+WX(I,KF)
  108.       DO 104 I=1,NLMW 
  109.   104 THWAL=THWAL+WX(I,KW)
  110.       ABASFT=ABASF(1)+ABASF(2)
  111.       HWBG=0. 
  112.       ABWBGT=0. 
  113.       DO 106 I=1,4
  114.       IZW=IZWE(I) 
  115.       ABWBGT=ABWBGT+AWALLA(IZW) 
  116.   106 HWBG=HWBG+HGTA(IZW) 
  117. C  NOTE: HEIGHT OF NODES IN BLW GRND BSMT WALL IS SET BY ENCL INPUT DATA
  118. C        SUGGESTION: WITH BASEMENT, SET CENTER OF NODE 8 TO 3 FT TO 
  119. C                    AGREE WITH MSMT OF SOIL TEMP IN HSE A AND B. 
  120. C                    WITH CRWSPC, SET NODE HEIGHTS = .1,.2,.3,.4 X
  121. C                    CRWSPC HEIGHT. 
  122.       DZN(1)=HGTA(IZWE(1))
  123.       DZN(2)=HGTA(IZWE(2))
  124.       DZN(8)=HGTA(IZWE(3))
  125.       DZN(9)=HGTA(IZWE(4))
  126.       IF(IGINIT.EQ.1)THEN 
  127.           DO 110 I=1,28 
  128.   110     TG(I)=TSOILI
  129.           ENDIF 
  130. C  SET UP NODAL GEOMETRY AND HEAT TRANSFER PARAMETERS 
  131.       CALL SETNOD(TIMEYR,HOGRND)
  132.       CALL SETND1(TIMEYR,HOGRND)
  133. C 
  134.       IF (IGINIT.EQ.2)THEN
  135. C  GET INITIAL GROUND TEMPS FROM TAPE43 
  136. C    TAPE41 TO FIRST PRECONDITIONING DAY AND ADVANCE POINTER ON TAPE41
  137. C    TO 1ST PRECONDITIONING DAY
  138.           OPEN(43,FILE='TAPE43',STATUS='UNKNOWN',IOSTAT=IO43) 
  139.           IF(IO43.NE.0)THEN
  140.             WRITE(60,*)' GRNDEX: CANT OPEN TAPE43 FILE'
  141.             STOP ' GRNDEX: CANT OPEN TAPE43 FILE'
  142.             ENDIF
  143.           READ(43,*,END=420)(TG(I),I=1,28)
  144.           CLOSE(43) 
  145.           NREAD=NDHRS*NMNPRE
  146.           DO 112 I=1,NREAD
  147. C              NMONTQ IS A DUMMY          
  148.           READ(41,503,END=430)NMONTQ 
  149.   112     CONTINUE
  150.           NMONTQ=NMNPRE 
  151.           NDAYT=NDHRS/24
  152.           NDAYTQ=NDAYT-1
  153.           ASSIGN 270 TO LAST1 
  154.           GO TO 260 
  155.           ENDIF 
  156. C 
  157. C  RUN GROUND PRECOND FOR NMPRE MONTHS
  158. C   FINAL NDAYS-1 WEATHER VALUES WILL BE STORED FOR REPEAT CALC 
  159. C     (NOTE:  ARRAYS ARE DIM FOR MAX OF NDAYS-1=4)
  160.   200 NMONTQ=1
  161.       ASSIGN 270 TO LAST1 
  162.       IF(NPRINT.GT.0)THEN 
  163.           DO 201 I=1,28 
  164.           IPQ=IP(I) 
  165.   201     TPR(I)=TG(IPQ)
  166.           WRITE(60,504)(TPR(I),I=1,28) 
  167.           ENDIF 
  168.       NDAYT=NDHRS/24
  169.       NDAYTQ=NDAYT-1
  170.   207 NDAY=0
  171.   208 NDAY=NDAY+1 
  172.       IF(NDAY.GT.NDAYT)GO TO 220
  173.       IF(NDAY.GT.1)NDY=NDAY-1 
  174. C  READ 24 HOURS OF DATA
  175.       DO 209 I=1,24 
  176. C         TPR(1) BELOW IS A DUMMY VAR      
  177.       READ(41,*)TDBT(I),TPR(1),WSPDT(I),SOLRT(I),TPR(1),MONLST 
  178. C  SAVE VALUES   FOR REPEAT CALC
  179.       IF(NDAY.GT.1.AND.NDY.LE.4)THEN 
  180.           ODDB(I,NDY)=TDBT(I) 
  181.           WSPEED(I,NDY)=WSPDT(I)
  182.           SOLARH(I,NDY)=SOLRT(I)
  183.           ENDIF 
  184.   209 CONTINUE
  185.       IHR=0 
  186.   210 IHR=IHR+1 
  187.       IF(IHR.EQ.25)THEN 
  188.           IF(NDBG.EQ.1)WRITE(60,505)NDAYQ
  189.           NDAYQ=NDAYQ+1 
  190.           GO TO 208 
  191.           ENDIF 
  192.       IF(IPRBG.EQ.1)THEN
  193.         CALL CURS(20)
  194.         WRITE(*,542)IHR,NDAY,NDAYQ,MONLST,KM
  195.   542 FORMAT(1X,'GRNDEX 1: INITIAL DAYS; IHR,NDAY,NDAYQ,MONLST,KM= ',
  196.      + 5I4)     
  197.       ENDIF
  198.       TODDB=TDBT(IHR) 
  199.       WSPED=WSPDT(IHR)
  200.       SOLRH=SOLRT(IHR)
  201. C  GET GROUND HEAT TRANSFER COEF
  202.       HOGRND=HGRN(WSPED)
  203.       CALL SETND2(TIMEYR,HOGRND)
  204.       CALL EXPLCT(TODDB,SOLRH,HOGRND,ABSRPG) 
  205.       GO TO 210 
  206. C 
  207. C  NOW REPEAT CALCULATIONS TO APPROX THE EFFECT ON GROUND 
  208. C   OF AN ENTIRE MONTH OF WEATHER 
  209. C  NOTE: ONLY ABOUT 20 DAYS ARE NEEDED, SO TOTAL DAYS/MON < MAX IS OK 
  210.   220 IEND=INT((NDYM(MONLST)-NDAYT)/NDAYTQ)
  211.       IRPEAT=0
  212.   230 IRPEAT=IRPEAT+1 
  213.       IF(IRPEAT.GT.IEND)GO TO 260 
  214.       NDAY=0
  215.   240 NDAY=NDAY+1 
  216.       IF(NDAY.GT.NDAYTQ)GO TO 230 
  217.       IHR=0 
  218.   250 IHR=IHR+1 
  219.       IF(IHR.EQ.25)THEN 
  220.           IF(NDBG.EQ.1)WRITE(60,505)NDAYQ
  221.           NDAYQ=NDAYQ+1 
  222.           GO TO 240 
  223.           ENDIF 
  224.       IF(IPRBG.EQ.1)THEN
  225.         CALL CURS(21)
  226.         IF(ICODE.EQ.4)THEN
  227.           WRITE(*,545)IHR,NDAY,NDAYQ,MONLST,KM
  228.         ELSE
  229.           WRITE(*,543)IHR,NDAY,NDAYQ,MONLST,KM
  230.         ENDIF
  231.   543 FORMAT(1X,'GRNDEX 2: REPEAT INITIAL DAYS; IHR,NDAY,NDAYQ,MONLST,',
  232.      +       'KM= ',5I4)     
  233.   545 FORMAT(1X,'GRNDEX 4: FINISH MONTH; IHR,NDAY,NDAYQ,MONLST,',
  234.      +       'KM= ',5I4)     
  235.       ENDIF
  236.       TODDB=ODDB(IHR,NDAY)
  237.       WSPED=WSPEED(IHR,NDAY)
  238.       SOLRH=SOLARH(IHR,NDAY)
  239.       HOGRND=HGRN(WSPED)
  240.       CALL SETND2(TIMEYR,HOGRND)
  241.       CALL EXPLCT(TODDB,SOLRH,HOGRND,ABSRPG) 
  242.       GO TO 250 
  243. C 
  244. C  NOW DO PRECONDITIONING FOR THE NEXT MONTH
  245.   260 IF(NPRINT.GT.0)THEN 
  246.           DO 262 I=1,28 
  247.           IPQ=IP(I) 
  248.   262     TPR(I)=TG(IPQ)
  249.           WRITE(60,504)(TPR(I),I=1,28)
  250.           ENDIF 
  251.       GO TO LAST1,(270,410) 
  252.   270 NMONTQ=NMONTQ+1 
  253.       IF(NMONTQ.LE.NMNPRE)GO TO 207 
  254.       IF(IGINIT.EQ.-2)THEN
  255.           OPEN(43,FILE='TAPE43',STATUS='UNKNOWN') 
  256.           WRITE(43,502)(TG(I),I=1,28)
  257.           CLOSE(43) 
  258.           ENDIF 
  259.       IHR=0 
  260.       NDY=0 
  261.       GO TO 500 
  262. C 
  263. C  HOURLY GROUND HEAT TRANSFER CALCULATIONS 
  264.   300 HOGRND=HGRN(WSPED)
  265.       CALL SETND2(TIMEYR,HOGRND)
  266.       CALL EXPLCT(TODDB,SOLRH,HOGRND,ABSRPG) 
  267. C 
  268. C  SAVE WEATHER DATA AFTER THE FIRST DAY
  269.   304 IHR=IHR+1 
  270.       IF(IHR.EQ.25)THEN 
  271.           IF(NDBG.EQ.1)WRITE(60,505)NDAYQ
  272.           NDAYQ=NDAYQ+1 
  273.           NDY=NDY+1 
  274.           IHR=0 
  275.           GO TO 304 
  276.         ENDIF 
  277.       IF(IPRBG.EQ.1)THEN
  278.         CALL CURS(22)
  279.         WRITE(*,544)IHR
  280.   544 FORMAT(16X,'GRNDEX 3: HOURLY CALCULATIONS FOR HR= ',I3)
  281.       ENDIF
  282.       IF(NDY.GE.1.AND.NDY.LE.4)THEN
  283.           ODDB(IHR,NDY)=TODDB 
  284.           WSPEED(IHR,NDY)=WSPED 
  285.           SOLARH(IHR,NDY)=SOLRH 
  286.           ENDIF 
  287.       GO TO 500 
  288. C 
  289. C  REPEAT CALC OF LAST FOUR DAYS TO FINISH MONTH
  290.   400 ASSIGN 410 TO LAST1 
  291. C  !!!
  292. C  USE HEAT FLOWS FROM WALLS AT LAST HOUR OF LAST DAY TO FINISH MONTH 
  293.       MONLST=KM 
  294. C         GET AVERAGES FOR THE PREVIOUS 24 HOURS
  295.       IF(NTBASM.EQ.0)STOP ' GRNDEX:NTBASM=0.! '      
  296.       DATA NTMQU/0/
  297.       IF(NTMQU.EQ.0)THEN
  298.         UBASFRF=UBASF
  299.         UBWBGRF=UBWBG
  300.         TBASMRF=TBASM
  301.         NTMQU=1
  302.       ENDIF
  303.       TBASM=TBASM/NTBASM
  304.       UBASF=0.
  305.       DO 402 I=1,2
  306.       IZW=IZWD(I)
  307.       IJ=15
  308.       IF(I.EQ.2)IJ=27
  309.       QCSURF(IZW)=QCSURFS(IZW)/NTBASM
  310.       UBASF=UBASF+(QCSURF(IZW)/(TBASM-TG(IJ)) - DDWKL(I+4))/AWALLA(IZW)
  311.       QCSURFS(IZW)=0.
  312.   402 CONTINUE
  313.       UBASF=0.5*UBASF
  314.       UBWBG=0.
  315.       DO 404 I=1,4
  316.       IZW=IZWE(I)
  317.       IJ=I
  318.       IF(I.EQ.3)IJ=8
  319.       IF(I.EQ.4)IJ=9
  320.       QCSURF(IZW)=QCSURFS(IZW)/NTBASM
  321.       UBWBG=UBWBG+(QCSURF(IZW)/(TBASM-TG(IJ)) - DDWKL(I))/AWALLA(IZW)
  322.       QCSURFS(IZW)=0.
  323.   404 CONTINUE
  324.       UBWBG=UBWBG/4.
  325.       WRITE(60,*)' GRNDEX 4: REPEAT CALC OF LAST 4 DAYS TO FINISH MONTH'
  326.       WRITE(60,*)' GRNDEX: ADJUSTING TBASM,UBASF,UBWBG VALUES..'
  327.       WRITE(60,*)' INPUT VALUES OF UBASF AND UBWBG= ',UBASFRF,UBWBGRF
  328.       WRITE(60,*)'  AT KM= ',KM,' NEW VALUES OF UBASF,UBWBG= ',UBASF,
  329.      +             UBWBG
  330.       WRITE(60,*)' INPUT VALUE OF BSMT AIR TEMP = ',TBASMRF
  331.       WRITE(60,*)'  MEAN VALUE OF BSMT AIR TEMP(LAST 24 HRS)= ',TBASM
  332.       WRITE(60,*)'  NTBASM= ',NTBASM
  333.       NTBASM=0
  334.       CALL CURS(23)
  335.       WRITE(*,*)' GRNDEX 4: REPEAT CALC OF LAST 4 DAYS TO FINISH MONTH'
  336.       GO TO 220 
  337.   410 CONTINUE
  338.       CALL CURS(21)
  339.       WRITE(*,*)' GRNDEX 4: COMPLETED FINISHING MONTH!                '
  340.       NDY=0 
  341.       IHR=0 
  342.       GO TO 500 
  343. C 
  344.   420 STOP 'GRNDEX:EOF READING TGS FROM TAPE43; SET IGINIT=-2 AND RERUN' 
  345.   430 STOP 'GRNDEX:EOF WHILE ADVANCING ON TAPE41' 
  346. C  ZERO HEAT FLOW ARRAYS
  347.   500 CALL ZERV(6,QBASF(1),QBASF(2),QBWBG(1),QBWBG(2),
  348.      +          QBWBG(3),QBWBG(4),0.,0.,0.,0.)
  349.       TBASM=0.
  350.       CALL CURS(20)
  351.       WRITE(*,506)
  352.       CALL CURS(21)
  353.       WRITE(*,506)
  354.       CALL CURS(23)
  355.       WRITE(*,506)
  356.       RETURN
  357.   999 STOP ' GROUND: EOF ON TAPE18 WITH NAMELIST INPGND'
  358. C 
  359.   502 FORMAT((1X,7(1X,E11.5)))
  360.   503 FORMAT(A10) 
  361.   504 FORMAT(1H0,'T S= ',8X,2G13.5/ 
  362.      +1X,39X,3G13.5/
  363.      +1X,13X,2G13.5/
  364.      +1X,13X,2G13.5/
  365.      +1X,39X,3G13.5/
  366.      +1X,13X,2G13.5/1X,2G13.5/1X,19X,G13.5,7X,G13.5/
  367.      +1X,2G13.5/
  368.      +1X,52X,2G13.5/1X,G13.5,13X,G13.5/ 
  369.      +1X,G13.5,19X,G13.5,20X,G13.5/1X,32X,G13.5)
  370.   505 FORMAT(1X,'GRNDEX:END OF DAY= ',I5) 
  371.   506 FORMAT(80(1H ))
  372.       END 
  373.       FUNCTION HGRN(WSPED)
  374. C 
  375. C  CONVECTIVE HEAT TRANSFER COEF AT GROUND DUE TO WIND VELOCITY 
  376. C    USE CARRIER CORRELATION FOR H AT GROUND SURF 
  377. C    (AN ALTERNATE WHICH GIVES SOMEWHAT LOWER VALUES IS THAT FOR
  378. C     A FLAT PLATE,IE: H=0.431*WSPED**.8 FOR L=10 FT. OR
  379. C     H*(1.0./L)**.2 FOR LOWER L) 
  380. C     UNITS: WSPED, MPH      HGRN, BTU/HR-FT2-F 
  381.       IF(WSPED.LT..4)THEN 
  382.           HGRN=1. 
  383.       ELSE
  384.           HGRN=SQRT(WSPED)/.6584
  385.           ENDIF 
  386.       RETURN
  387.       END 
  388.       SUBROUTINE SETNOD(TIMEYR,HOGRND)
  389. C 
  390. C  SET NODE VALUES FOR BASEMENT TO GROUND HEAT TRANSFER 
  391. C 
  392.       DIMENSION         DRN(28),RR(5),DN(28),DNO(28),VOL(28), 
  393.      +DDZ1(28),DDZ2(28),DDR1(28),DDR2(28),
  394.      +DDZ(6),IR(28),IZ(28),IA(28),A(16),DDWL(6),DRWL(6)
  395.      +,XK(28),I19(4),DDZK(7),I17(6)
  396. C 
  397. CMDK BLK52
  398. CMDK BLK53
  399. CMDK BLK55
  400. CMDK BLK57
  401. CMDK BLKHWG
  402. CMDK ENCBK2
  403. CMDK SOILB
  404. CMDK UABAS
  405.       LOGICAL NDEBG 
  406. C 
  407. C  F1F GIVES THE AVG RADIUS THAT HEAT MUST TRAVEL (ref: ANDREWS,
  408. C  A TRNSYS-COMPATIBLE MODEL OF GROUND-COUPLED STORAGE,SEPT 1979,
  409. C  BNL 51061 UC-59c)
  410.       F1F(Z1,Z2)=2./3.*(Z2**3-Z1**3)/(Z2**2-Z1**2)
  411.       F2F(Z3,Z4)=PI4*(Z4*Z4-Z3*Z3)
  412. C SETUP RADIAL CONNECTIONS (IR IS RADIAL CONNECTED NODE FOR NODE I
  413.       DATA IR/3,4,5,5,6,7,0,10,11,12,12,13,14,0,27, 
  414.      +28,18,19,20,0,22,19,24,25,0,0,17,17/
  415. C  SETUP AXIAL CONNECTIONS
  416.       DATA IZ/2,8,4,10,12,13,14,9,17,11,17,18,
  417.      +19,20,16,21,22,22,24,25,23,24,26,26,26,0,28,21/ 
  418. C  SET AREA INDICATOR 
  419.       DATA IA/1,1,2,2,3,4,6,1,1,
  420.      +2,2,3,4,6,14,14,7,3,4,6,5,
  421.      +8,5,9,6,0,15,15/
  422.       DATA I19/1,2,8,9/ 
  423.       DATA I17/1,3,5,6,0,7/ 
  424.       NDEBG=.FALSE. 
  425.       DATA PI/3.141596/ 
  426.       PI2=2.*PI 
  427.       PI4=PI/4. 
  428.       DO 10 I=1,30
  429.       DO 10 J=1,30
  430.    10 AK(I,J)=0.
  431.       REQBAS=SQRT(ABASFT/PI)
  432.       DEQBAS=2.*REQBAS
  433.       ABWBGC=PI*DEQBAS*HWBG 
  434. C   RATH IS AMPLIFICATION OF HEAT FLOW TO BASEMENT AIR AT WALL
  435. C     TO ACCOUNT  FOR SMALLER WALL AREA WITH CYL COORDINATES
  436.       RATH=ABWBGC/ABWBGT
  437.       REQ1=SQRT(REQBAS**2. - ABASF(2)/PI)
  438.       DEQ1=2.*REQ1
  439.       REQ2=REQBAS+THWAL 
  440.       A(14)=PI*REQ1*REQ1
  441.       A(15)=PI*(REQ2*REQ2-REQ1*REQ1)
  442.       A(16)=ABASFT-A(14)
  443.       DZN(5)=DZN(1)+DZN(2)
  444.       DZN(6)=DZN(5) 
  445.       DZN(7)=DZN(5) 
  446.       DZN(12)=DZN(8)+DZN(9) 
  447.       DZN(13)=DZN(12) 
  448.       DZN(14)=DZN(12) 
  449.       DZN(3)=DZN(1) 
  450.       DZN(4)=DZN(2) 
  451.       DZN(10)=DZN(8)
  452.       DZN(11)=DZN(9)
  453.       ZTF=ZDEPTH-HWBG-THFLR 
  454.       DZN(15)=ZTF/22. 
  455.       DZN(16)=3.*ZTF/22.
  456.       DZN(17)=DZN(15)+DZN(16) 
  457.       DZN(18)=DZN(17) 
  458.       DZN(21)=6.*ZTF/22.
  459.       DZN(22)=DZN(21) 
  460.       DZN(19)=DZN(18)+DZN(22) 
  461.       DZN(20)=DZN(19) 
  462.       DZN(23)=12.*ZTF/22. 
  463.       DZN(24)=DZN(23) 
  464.       DZN(25)=DZN(24) 
  465.       DZN(26)=2.*DZN(23)
  466.       DZN(27)=DZN(15) 
  467.       DZN(28)=DZN(16) 
  468.       IF(IGINIT.EQ.1)GO TO 16 
  469. C   SET INITIAL SOIL TEMPS USING KUSUDA RELATIONSHIP
  470. C     SET DEPTH TO CENTER OF NODES(ZN)
  471. C     (USED TO SET INITIAL GROUND TEMPS)
  472.       ZN(1)=0.5*DZN(1)
  473.       ZN(2)=DZN(1)+0.5*DZN(2) 
  474.       ZN(3)=ZN(1) 
  475.       ZN(4)=ZN(2) 
  476.       ZN(8)=DZN(1)+DZN(2)+0.5*DZN(8)
  477.       ZN(9)=DZN(1)+DZN(2)+DZN(8)+0.5*DZN(9) 
  478.       ZN(10)=ZN(8)
  479.       ZN(11)=ZN(9)
  480.       ZN(5)=0.5*DZN(5)
  481.       ZN(12)=DZN(5)+0.5*DZN(12) 
  482.       ZN(6)=ZN(5) 
  483.       ZN(13)=ZN(12) 
  484.       ZN(7)=ZN(5) 
  485.       ZN(14)=ZN(12) 
  486.       ZN(15)=HWBG+THFLR+0.5*DZN(15) 
  487.       ZN(16)=HWBG+THFLR+DZN(15)+0.5*DZN(16) 
  488.       ZN(17)=DZN(5)+DZN(12)+0.5*DZN(17) 
  489.       ZN(18)=ZN(17) 
  490.       ZN(19)=ZN(6)+ZN(13)+0.5*DZN(19) 
  491.       ZN(20)=ZN(19) 
  492.       ZN(21)=DZN(6)+DZN(13)+DZN(17)+0.5*DZN(21) 
  493.       ZN(22)=ZN(21) 
  494.       ZN(23)=DZN(6)+DZN(13)+DZN(19)+0.5*DZN(23) 
  495.       ZN(24)=ZN(23) 
  496.       ZN(25)=ZN(24) 
  497.       ZN(26)=DZN(6)+DZN(13)+DZN(19)+DZN(24)+0.5*DZN(24) 
  498.       ZN(27)=ZN(15) 
  499.       ZN(28)=ZN(16) 
  500.       ALFZ=XKSOIL/(RHSOIL*CPSOIL) 
  501.       IF(IGINIT.NE.2)THEN
  502.         DO 12 I=1,28
  503.         ZNODE=ZN(I) 
  504.         TG(I)=TGKUS(TIMEYR)
  505.    12   CONTINUE
  506.       ENDIF
  507. C   SET RADIAL WIDTHS OF NODES
  508.    16 DRN(1)=RFIELD/28. 
  509.       DRN(2)=DRN(1) 
  510.       DRN(8)=DRN(1) 
  511.       DRN(9)=DRN(1) 
  512.       DRN(3)=3.*RFIELD/28.
  513.       DRN(4)=DRN(3) 
  514.       DRN(10)=DRN(3)
  515.       DRN(11)=DRN(3)
  516.       DRN(5)=6.*RFIELD/28.
  517.       DRN(12)=DRN(5)
  518.       DRN(17)=DRN(1)+DRN(3) 
  519.       DRN(18)=DRN(5)
  520.       DRN(22)=DRN(17)+DRN(18) 
  521.       DRN(6)=18.*RFIELD/28. 
  522.       DRN(13)=DRN(6)
  523.       DRN(19)=DRN(6)
  524.       DRN(24)=DRN(22)+DRN(19) 
  525.       DRN(15)=REQ1
  526.       DRN(16)=DRN(15) 
  527.       DRN(21)=REQ2
  528.       DRN(23)=DRN(21) 
  529.       DRN(25)=0.
  530.       DRN(26)=0.
  531.       DRN(27)=REQ2-REQ1 
  532.       DRN(28)=DRN(27) 
  533.       RR(1)=REQ2
  534.       RR(2)=RR(1)+DRN(1)
  535.       RR(3)=RR(2)+DRN(3)
  536.       RR(4)=RR(3)+DRN(5)
  537.       RR(5)=RR(4)+DRN(6)
  538. C   DN VALUES ARE DIA TO NODE CENTROID
  539. C   DNO VALUES ARE DIA TO OUTER WALL OF NODE
  540.       DNO(1)=2.*RR(2) 
  541.       D=2.*RR(1)
  542.       DN(1)=F1F(D,DNO(1)) 
  543.       DN(2)=DN(1) 
  544.       DN(8)=DN(1) 
  545.       DN(9)=DN(1) 
  546.       DNO(2)=DNO(1) 
  547.       DNO(8)=DNO(1) 
  548.       DNO(9)=DNO(1) 
  549.       DNO(3)=DNO(1)+2.*DRN(3) 
  550.       DNO(4)=DNO(3) 
  551.       DNO(10)=DNO(3)
  552.       DNO(11)=DNO(3)
  553.       DN(3)=F1F(DNO(1),DNO(3))
  554.       DN(4)=DN(3) 
  555.       DN(10)=DN(3)
  556.       DN(11)=DN(3)
  557.       DNO(5)=DNO(3)+2.0*DRN(5)
  558.       DNO(12)=DNO(5)
  559.       DNO(18)=DNO(5)
  560.       DNO(22)=DNO(5)
  561.       DN(5)=F1F(DNO(3),DNO(5))
  562.       DN(12)=DN(5)
  563.       DN(18)=DN(5)
  564.       DNO(17)=DNO(3)
  565.       DN(17)=F1F(D,DNO(17)) 
  566.       DNO(6)=DNO(5)+2.*DRN(6) 
  567.       DNO(13)=DNO(6)
  568.       DNO(19)=DNO(6)
  569.       DNO(24)=DNO(6)
  570.       DN(6)=F1F(DNO(5),DNO(6))
  571.       DN(13)=DN(6)
  572.       DN(19)=DN(6)
  573.       DN(22)=F1F(D,DNO(22)) 
  574.       DN(15)=2./3.*DEQ1 
  575.       DN(16)=DN(15) 
  576.       DN(21)=2./3.*D
  577.       DN(23)=DN(21) 
  578.       DNO(15)=DEQ1
  579.       DNO(16)=DNO(15) 
  580.       DNO(21)=D 
  581.       DNO(23)=DNO(21) 
  582.       DN(24)=F1F(DNO(23),DNO(24)) 
  583.       DN(7)=DNO(6)+DRN(6) 
  584.       DN(14)=DN(7)
  585.       DN(20)=DN(7)
  586.       DN(25)=DN(7)
  587.       DNO(27)=DNO(21) 
  588.       DNO(28)=DNO(27) 
  589.       DN(27)=F1F(DNO(15),DNO(27)) 
  590.       DN(28)=DN(27) 
  591. C  CALC HEAT FLOW AREAS IN DEPTH DIRECTION
  592.       A(1)=F2F(D,DNO(1))
  593.       A(2)=F2F(DNO(1),DNO(3)) 
  594.       A(3)=F2F(DNO(3),DNO(5)) 
  595.       A(4)=F2F(DNO(5),DNO(6)) 
  596.       A(5)=PI4*D*D
  597.       A(6)=6.*A(4)
  598.       A(7)=A(1)+A(2)
  599.       A(8)=A(7)+A(3)
  600.       A(9)=A(8)+A(4)
  601.       A(10)=PI2*REQBAS*RATH*DZN(1)
  602.       A(11)=PI2*REQBAS*RATH*DZN(2)
  603.       A(12)=PI2*REQBAS*RATH*DZN(8)
  604.       A(13)=PI2*REQBAS*RATH*DZN(9)
  605. C  CALC NODE VOLUMES
  606.       DO 20 I=1,28
  607.       IF(I.EQ.26)GO TO 20 
  608.       L=IA(I) 
  609.       AZ=A(L) 
  610.       VOL(I)=AZ*DZN(I)
  611.    20 CONTINUE
  612.       DO 30 I=1,28
  613.       IF(I.EQ.26)GO TO 30 
  614.       K=IZ(I) 
  615.       L=IA(I) 
  616.       AZ=A(L) 
  617.       DDZ1(I)=0.5*DZN(I)/AZ 
  618.       DDZ2(I)=0.5*DZN(K)/AZ 
  619.    30 CONTINUE
  620. C  CALC DDZ'S FOR SURFACE CONNECTIONS 
  621.       DDZ(1)=0.5*DZN(1)/A(1)
  622.       DDZ(2)=0.5*DZN(3)/A(2)
  623.       DDZ(3)=0.5*DZN(5)/A(3)
  624.       DDZ(4)=0.5*DZN(6)/A(4)
  625.       DDZ(5)=0. 
  626.       DDZ(6)=0.5*DZN(7)/A(6)
  627. C   CALC DDR'S FOR RADIAL CONDUCTION
  628.       DO 40 I=1,28
  629.       IF(I.EQ.26)GO TO 40 
  630.       IF(I.EQ.7.OR.I.EQ.14)GO TO 40 
  631.       IF(I.EQ.20.OR.I.EQ.25)GO TO 40
  632.       J=IR(I) 
  633.       DENOM=PI2*DZN(I)
  634.       DDR1(I)=ALOG(DNO(I)/DN(I))/DENOM
  635.       DDR2(I)=ALOG(DN(J)/DNO(I))/DENOM
  636.    40 CONTINUE
  637. C 
  638. C  DDWL IS FROM OUTER SURF OF BASEMENT WALL TO ADJACENT GROUND CENTROID
  639. C 
  640.       DO 50 I=1,4 
  641.       L=I 
  642.       IF(I.EQ.3)L=8 
  643.       IF(I.EQ.4)L=9 
  644. C           DRWL IS THE RADIAL DISTANCE HEAT HAS TO FLOW FROM OUTER
  645. C           SURF OF BASEMENT WALL TO CENTROID OF ADJACENT GRND NODE      
  646.       DRWL(I)=0.5*(DN(L) - D)
  647.    50 DDWL(I)=ALOG(DN(L)/D)/(PI2*DZN(L))
  648.       DDWL(5)=0.5*DZN(15)/A(14)
  649.       DDWL(6)=0.5*DZN(27)/A(15) 
  650.       RETURN
  651. C 
  652. C * * * * 
  653.       ENTRY SETND1(TIMEYR,HOGRND) 
  654. C    CALC AK'S (CONDUCTANCES BETWEEN NODES) 
  655. C  SETUP FOR VARIABLE K, IF NEEDED LATER
  656.       DO 60 I=1,28
  657.    60 XK(I)=XKSOIL
  658.       DDWKL(1)=DDWL(1)/XK(1)
  659.       DDWKL(2)=DDWL(2)/XK(2)
  660.       DDWKL(3)=DDWL(3)/XK(8)
  661.       DDWKL(4)=DDWL(4)/XK(9)
  662. C         CONDWL(1-4) VALUES ARE USED IN OWALL.FOR      
  663.       CONDWL(1)=XK(1)/DRWL(1)
  664.       CONDWL(2)=XK(2)/DRWL(2)
  665.       CONDWL(3)=XK(8)/DRWL(3)
  666.       CONDWL(4)=XK(9)/DRWL(4)
  667.       DDWKL(5)=DDWL(5)/XK(15) 
  668.       DDWKL(6)=DDWL(6)/XK(27) 
  669. C         CONDWL(5,6) VALUES ARE USED IN FLOOR.FOR      
  670.       CONDWL(5)=XK(15)/(0.5*DZN(15))
  671.       CONDWL(6)=XK(27)/(0.5*DZN(27))
  672.       DO 62 I=1,6 
  673.       IF(I.EQ.5)GO TO 62
  674.       IQ=I17(I) 
  675.       DDZK(I)=DDZ(I)/XK(IQ) 
  676.    62 CONTINUE
  677. C  CALC INNER AK'S
  678.       DO 80 I=1,28
  679.       IF(I.EQ.26)GO TO 80 
  680.       J=IR(I) 
  681.       IF(J.EQ.0)GO TO 70
  682.       AK(I,J)=1./(DDR1(I)/XK(I)+DDR2(I)/XK(J))
  683.       AK(J,I)=AK(I,J) 
  684.    70 K=IZ(I) 
  685.       AK(I,K)=1./(DDZ1(I)/XK(I)+DDZ2(I)/XK(K))
  686.       AK(K,I)=AK(I,K) 
  687.    80 CONTINUE
  688. C   CALC BOUNDARY AK'S
  689. C  BOUNDARY AT BSMT FLR OR WALL IS SOIL TEMP OUTSIDE FLR OR WALL
  690.       AK(15,29)=1./DDWKL(5) 
  691.       AK(29,15)=AK(15,29) 
  692.       AK(27,29)=1./DDWKL(6) 
  693.       AK(29,27)=AK(27,29) 
  694.       DO 82 I=1,4 
  695.       IJ=I19(I) 
  696.       AK(IJ,29)=1./DDWKL(I) 
  697.       AK(29,IJ)=AK(IJ,29) 
  698.    82 CONTINUE
  699. C   CALC THERMAL CAPACITANCE VALUES (S'S) 
  700.       RCP=RHSOIL*CPSOIL 
  701.       DO 90 I=1,28
  702.       IF(I.EQ.26)GO TO 90 
  703.       RCPV=RCP*VOL(I) 
  704.       S(I)=1./RCPV
  705.    90 CONTINUE
  706. C 
  707.       IF(NDEBG)THEN 
  708.           WRITE(60,500)(DZN(I),I=1,28) 
  709.            WRITE(60,512)(ZN(I),I=1,28) 
  710.           WRITE(60,501)(DRN(I),I=1,28) 
  711.           WRITE(60,502)(RR(I),I=1,5) 
  712.           WRITE(60,503)(DN(I),DNO(I),I=1,28) 
  713.           WRITE(60,511)(A(I),I=1,15) 
  714.           WRITE(60,504)(VOL(I),I=1,28) 
  715.           WRITE(60,505)(DDZ1(I),DDZ2(I),I=1,28)
  716.           WRITE(60,506)(DDR1(I),DDR2(I),I=1,28)
  717.           WRITE(60,507)(DDZ(I),I=1,6)
  718.           WRITE(60,508)(DDWL(I),I=1,6) 
  719.           WRITE(60,510)(S(I),I    =1,28) 
  720.           ENDIF 
  721.       NDEBG=.FALSE. 
  722.       RETURN
  723. C 
  724.       ENTRY SETND2(TIMEYR,HOGRND) 
  725.       AK(1,30)=1./(DDZK(1)+1./(HOGRND*A(1)))
  726.       AK(30,1)=AK(1,30) 
  727.       AK(3,30)=1./(DDZK(2)+1./(HOGRND*A(2)))
  728.       AK(30,3)=AK(3,30) 
  729.       AK(5,30)=1./(DDZK(3)+1./(HOGRND*A(3)))
  730.       AK(30,5)=AK(5,30) 
  731.       AK(6,30)=1./(DDZK(4)+1./(HOGRND*A(4)))
  732.       AK(30,6)=AK(6,30) 
  733.       AK(7,30)=1./(DDZK(6)+1./(HOGRND*A(6)))
  734.       AK(30,7)=AK(7,30) 
  735.       IF(ICODE.EQ.2)THEN
  736. C         SET PARAMETERS FOR HEAT FLOW FROM BSMT AIR TO FLOOR AND WALLS      
  737.       ENDIF
  738.       RETURN
  739.   500 FORMAT(1X,'DZN S= ',/6(1X,5G13.5/)) 
  740.   501 FORMAT(1X,'DRN S= '/6(1X,5G13.5/))
  741.   502 FORMAT(1X,'RR S= '/5G13.5)
  742.   503 FORMAT(1X,'DN,DNO S= '/5(F7.2,1X,F7.2)/ 
  743.      +5(F7.2,1X,F7.2)/5(F7.2,1X,F7.2))
  744.   504 FORMAT(1X,'VOL S= '/6(1X,5G13.5/))
  745.   505 FORMAT(1X,'DDZ1,DDZ2= '/5(1X,G12.4,G12.4)/
  746.      +5(1X,G12.4,G12.4)/5(1X,G12.4,G12.4))
  747.   506 FORMAT(1X,'DDR1,DDR2= '/5(1X,G12.4,G12.4)/
  748.      +5(1X,G12.4,G12.4)/5(1X,G12.4,G12.4))
  749.   507 FORMAT(1X,'DDZ S= ',6G13.5) 
  750.   508 FORMAT(1X,'DDWL S= ',6G13.5)
  751.   509 FORMAT(11X,'I,J= ',2I5,'AK= ',G13.5)
  752.   510 FORMAT(1X,'S S= ',6(1X,5G13.5/))
  753.   511 FORMAT(1X,'A= ',5G13.5/4X,5G13.5/4X,6G13.5) 
  754.   512 FORMAT(1X,'ZN S= ',/6(1X,5G13.5/))
  755.       END 
  756.       SUBROUTINE EXPLCT(TODDB,SOLARH,HOGRND,ABSRPG)
  757. C 
  758. C  COMPUTE NEW TEMPS FOR GROUND AROUND BASEMENT USING 
  759. C          EXPLICIT,FORWARD,FINITE DIFFERENCE TECHNIQUE 
  760. C  --METHOD ORIGINATED WITH G. WHITACRE AT BATTELLE-- 
  761. C    THERMAL CONNECTIONS BETWEEN NODES ARE REFERENCED TO TOP,RIGHT
  762. C          BOTTOM,AND LEFT IN NODAL SKETCH
  763.       DIMENSION ITOP(28,2),ILFT(28,2),IRGT(28),IBOT(28) 
  764.      +,JTOP(28),JLFT(28),ISL(6)
  765. C      REAL STAB(28)
  766. CMDK NZN
  767. CMDK NZW
  768. CMDK IZWQ
  769. CMDK BLK52
  770. CMDK SOILB
  771. CMDK TEMPB
  772. CMDK UABAS
  773.       DATA ITOP/30, 1, 30,  3, 30, 30, 30,  2,  8,  4,
  774.      +          10,  5,  6,  7, 29, 15,  9, 12, 13, 14, 
  775.      +          16, 17, 21, 22, 20,  0, 29, 27, 16*0,11,3*0,28,18,0,19, 
  776.      +          4*0/
  777.       DATA IRGT/ 3,  4,  5,  5,  6,  7,  0, 10, 11, 12, 
  778.      +          12, 13, 14,  0, 27, 28, 18, 19, 20,  0, 
  779.      +          22, 19, 24, 25,  0,  0, 17, 17/ 
  780.       DATA IBOT/ 2,  8,  4,  10, 12, 13, 14,  9, 17, 11,
  781.      +          17, 18, 19, 20, 16, 21, 22, 22, 24, 25, 
  782.      +          23, 24, 26, 26, 26,  0, 28, 21/ 
  783.       DATA ILFT/29, 29,  1,  2,  3,  5,  6, 29, 29,  8, 
  784.      +           9, 10, 12, 13,  0,  0, 27, 17, 18, 19, 
  785.      +           0, 21,  0, 23, 24,  0, 15, 16, 4*0,4,6*0,11,4*0,28,0,
  786.      +            22,7*0,15,16/ 
  787.       DATA NU,NS1,NS2/1,1,0/
  788.       DATA JTOP/16*1,2,3*1,2,2,1,2,1,0,1,1/ 
  789.       DATA JLFT/4*1,2,6*1,2,2*1,2*0,2,1,2,1,0,1,0,2*1,0,1,1/
  790.       DATA SMAX/1./ 
  791.       DATA ISL/1,2,8,9,15,27/ 
  792.       DATA NV/0/,NSQ/0/ 
  793. C  NOTE: TIME BETWEEN CALLS TO THIS SUBPROGRAM SHOULD BE 1.0 HR (SET BELOW) 
  794.       DATA DTIMEG/1.0/ 
  795. C  SET SOIL TEMPS AT BACK-SIDE OF BSMT FLOOR AND BLW GRND BSMT WALL 
  796.       IF(ICODE.EQ.3)THEN
  797.         DO 2 I=1,6
  798.         IL=ISL(I) 
  799.         ISOIL=I 
  800.         IF(I.EQ.6)ISOIL=11
  801.         IF(I.LE.4)THEN
  802.           IZ=IZWE(I)
  803.         ELSE
  804.           IZ=IZWD(I-4)
  805.         ENDIF
  806.         TSOIL(ISOIL)=TMP(1,IZ)
  807.     2   CONTINUE      
  808.       ELSE
  809. C         SET TSOIL AT BSMT WALLS AND FLOOR AND UPDATE NODE TEMPS
  810. C         DURING PRECONDITIONING AND END-OF-MONTH CALCULATIONS.
  811.         CALL BSMTPRE
  812.       ENDIF
  813.       TG(30)=TODDB+ABSRPG*SOLARH/HOGRND 
  814. C   LOOP STARTS HERE
  815.       SMIN=20.
  816.       IQ=0
  817.    10 IQ=IQ+1 
  818.       IF(IQ.EQ.7)IQ=IQ+1
  819.       IF(IQ.EQ.14)IQ=IQ+1
  820.       IF(IQ.EQ.20)IQ=IQ+1
  821.       IF(IQ.EQ.25)IQ=IQ+2
  822.       IF(IQ.EQ.15)TG(29)=TSOIL(5) 
  823.       IF(IQ.EQ.27)TG(29)=TSOIL(11)
  824.       AKT=0.
  825.       SUM=0.
  826.       SUM1=0. 
  827.       JT=JTOP(IQ) 
  828.       DO 20 L=1,JT
  829.       IT=ITOP(IQ,L) 
  830.       IF(IT.EQ.0)GO TO 20 
  831.       AKT=AK(IQ,IT) 
  832.       IF(NU.EQ.0)WRITE(60,500)IQ,IT,AKT
  833.       SUM=SUM+AKT 
  834.       SUM1=SUM1+AKT*TG(IT)
  835.    20 CONTINUE
  836. C 
  837.       IR=IRGT(IQ) 
  838.       IF(IR.EQ.0)GO TO 30 
  839.       AKR=AK(IQ,IR) 
  840.       IF(NU.EQ.0)WRITE(60,500)IQ,IR,AKR
  841.       SUM=SUM+AKR 
  842.       SUM1=SUM1+AKR*TG(IR)
  843.    30 CONTINUE
  844. C 
  845.       IB=IBOT(IQ) 
  846.       IF(IB.EQ.0)GO TO 32 
  847.       AKB=AK(IQ,IB) 
  848.       IF(NU.EQ.0)WRITE(60,500)IQ,IB,AKB
  849.       SUM=SUM+AKB 
  850.       SUM1=SUM1+AKB*TG(IB)
  851.    32 CONTINUE
  852. C 
  853.       JR=JLFT(IQ) 
  854.       IF(JR.EQ.0)GO TO 50 
  855.       DO 40 L=1,JR
  856.       IL=ILFT(IQ,L) 
  857.       IF(IL.EQ.0)GO TO 40 
  858.       AKL=AK(IQ,IL) 
  859.       IF(NU.EQ.0)WRITE(60,500)IQ,IL,AKL
  860.       SUM=SUM+AKL 
  861.       IF(IQ.EQ.1)TG(29)=TSOIL(1)
  862.       IF(IQ.EQ.2)TG(29)=TSOIL(2)
  863.       IF(IQ.EQ.8)TG(29)=TSOIL(3)
  864.       IF(IQ.EQ.9)TG(29)=TSOIL(4)
  865.       SUM1=SUM1+AKL*TG(IL)
  866.    40 CONTINUE
  867.    50 TCN=DTIMEG*S(IQ) 
  868.       STB=1.-TCN*SUM
  869.       IF(STB.GT.SMIN)GO TO 60 
  870.       SMIN=STB
  871.       ISMIN=IQ
  872.    60 TG(IQ)=TG(IQ)*STB+TCN*SUM1
  873.       IF(IQ.LT.28)GO TO 10
  874.       TSOIL(6)=TG(7)
  875.       TSOIL(7)=TG(14) 
  876.       TSOIL(8)=TG(20) 
  877.       TSOIL(9)=TG(25) 
  878.       TSOIL(10)=TG(26)
  879.       NU=1
  880. C 
  881. C   END OF LOOP 
  882. C 
  883. C  CHECK TO SEE IF ALL GROUND TEMPS ARE STABLE
  884.       IF(SMIN.GT.SMAX)GO TO 70
  885.       IF(SMIN.LT.0.)GO TO 80
  886.       GO TO 90
  887.    70 IF(NS1.EQ.0)THEN
  888.       WRITE(60,501)ISMIN,SMIN,DTIMEG
  889.       NS1=1 
  890.       GO TO 90
  891.       ENDIF 
  892.    80 IF(NS2.EQ.0)THEN
  893.       WRITE(60,502)ISMIN,SMIN,DTIMEG
  894.       NS2=1 
  895.       NSQ=1 
  896.       GO TO 92
  897.       ENDIF 
  898.    90 CONTINUE
  899.       IF(NV.EQ.0.AND.NSQ.EQ.0)WRITE(60,503)ISMIN,SMIN
  900.       NV=1
  901.    92 CONTINUE
  902.       RETURN
  903.   500 FORMAT(1X,'I,J,AK= ',2I5,G13.5) 
  904.   501 FORMAT(1X,'OKAY TO INCR TIMSTP FOR GROUND T CALC'/
  905.      +1X,'ISMIN,SMAX,DTIME= ',I5,2G13.5)
  906.   502 FORMAT(1X,'--OOPS,TIMSTP FOR GROUND T CALC MUST BE DECR'/ 
  907.      +1X,'ISMIN,SMIN,DTIME= ',I5,2G13.5)
  908.   503 FORMAT(1H0,'*** STABLE GROUND-TEMP CALC---IQ,SMIN= ',I5,G13.5)
  909.       END 
  910.       FUNCTION TGKUS(TIME)
  911. C 
  912. C  KUSUDA RELATIONSHIP FOR FAR-FIELD GROUND TEMP
  913. C    ASHRAE TRANS VOL71(1) 1965,P61-75
  914. CMDK BLK55
  915. C 
  916.       DATA NTQ/0/ 
  917.       DATA PI/3.141596/ 
  918.       DATA HRYR/8766./
  919.       IF(NTQ.NE.0)GO TO 10
  920.       NTQ=1 
  921.       PI2H=2.*PI/HRYR 
  922.       C1=SQRT(PI/(ALFZ*HRYR)) 
  923.    10 C2=C1*ZNODE 
  924.       C3=PI2H*TIME-C2-PHASE 
  925.       TGKUS=TEAV-TEAMP/EXP(C2)*COS(C3)
  926.       RETURN
  927.       END 
  928.       SUBROUTINE BSMTPRE
  929. C
  930. C         CALC TEMPS OF BSMT WALLS AND FLOOR DURING PRECONDITIONING
  931. C         AND END OF MONTH CALCULATIONS.
  932. CMDK NWN
  933. CMDK NZN
  934. CMDK NZW
  935. CMDK IZWQ
  936. CMDK BLKBSF
  937. CMDK ENCBLK
  938. CMDK ENCBK2
  939. CMDK BLKQS
  940. CMDK SOILB
  941. CMDK SURFAR
  942. CMDK TEMP1
  943. CMDK TEMPB
  944. CMDK UABAS
  945.       NRM=NZNBAS
  946.       ID=9
  947. C              UPDATE BSMT WALL TEMPS BELOW GROUND      
  948. C              MUST DO THIS 3600/120 TIMES! (GROUND IS ON 1 HR 
  949. C              TIMESTEP, WHILE OWALLS IS ON A 2 MIN TIMESTEP
  950. C         NTGHR IS 3600/120      
  951.       DO 10 IU=1,NTGHR
  952.       DO 2 I = 1,4
  953. C         BSMT WALL NODES BELOW GROUND      
  954.       IZW=IZWE(I)
  955.       IJ=I
  956.       IF(I.GT.2)IJ=I+5
  957.       RB=1./(UBWBG * AWALLA(IZW) + DDWKL(I))
  958.       QCSURF(IZW)=1./RB * (TBASM - TG(IJ))
  959.       TSOIL(I)=TG(IJ) + QCSURF(IZW) * DDWKL(I)
  960.     2 CONTINUE
  961. C      DATA IQUU/0/
  962. C      IQUU=IQUU+1
  963.       DO 4 I=1,4
  964.       IZW=IZWE(I)
  965.       IJ=I
  966.       IF(I.GT.2)IJ=I+5
  967.       KON=KONSTA(IZW)
  968.       CALL OWALLS(NRM,IZW,ID,KON,DM,DM,DM,DM,IDM)
  969. C         UPDATE TMP VALUES (IS DONE IN ENDLP2 WHEN NOT PRECONDITIONING)      
  970.       DO 3 J=1,4
  971.     3 TMP(J,IZW)=TMPN(J,IZW)
  972.     4 CONTINUE
  973.    10 CONTINUE
  974. C   
  975.       ID=8
  976.       DO 20 IU=1,NTGHR
  977.       DO 12 I = 1, 2
  978. C         BASEMENT FLOOR NODES      
  979.       IZW=IZWD(I)
  980.       KON=KONSTA(IZW)
  981.       IF(I.EQ.1)THEN
  982.         IJ=15
  983.         ISOIL = 5
  984.       ELSEIF(I.EQ.2)THEN
  985.         IJ=27
  986.         ISOIL = 11
  987.       ENDIF
  988.       RB=1./(UBASF * AWALLA(IZW)  + DDWKL(I+4))
  989.       QCSURF(IZW) = 1./RB * (TBASM - TG(IJ))
  990.       TSOIL(ISOIL)=TG(IJ) + QCSURF(IZW) * DDWKL(I+4)
  991.    12 CONTINUE
  992.       DO 14 I = 1, 2
  993.       IZW=IZWD(I)
  994.       KON=KONSTA(IZW)
  995.       IF(I.EQ.1)THEN
  996.         ISOIL = 5
  997.       ELSEIF(I.EQ.2)THEN
  998.         ISOIL = 11
  999.       ENDIF
  1000.       CALL FLOOR(NRM,IZW,ID,KON,IDM,TSOIL(ISOIL),DM,IDM)
  1001.       DO 13 J=1,4
  1002.    13 TMP(J,IZW)=TMPN(J,IZW)
  1003.    14 CONTINUE
  1004.    20 CONTINUE
  1005.       RETURN
  1006.       END      
  1007.