1. C**** 2. C**** C089P.S Fortran source code for Physics routines 2003/11/28 3. C**** 4. C**** NASA/GISS Climate Model III Programmed by: 5. C**** Gary L. Russell 6. C**** NASA/Goddard Space Flight Center Gavin A. Schmidt 7. C**** Institute for Space Studies Reto A. Ruedy 8. C**** 2880 Broadway, New York, NY 10025 9. C**** U.S.A. 10. C**** 11. SUBROUTINE MSTCNV 12. C**** 13. C**** CONDSE loops over horizontal quarter boxes. For each quarter 14. C**** box, CONDSE loads one-dimesnional arrays of mass and heat and 15. C**** then calls MSTCNV and LSCOND. The returned arrays then 16. C**** reestablish the means and gradients of the prognostic variables. 17. C**** 18. INCLUDE 'C070.COM' 19. LOGICAL*4 QPOLE 20. PARAMETER (FPLUME=.25, TKF=273.16d0, TI=233.16d0) 21. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 22. COMMON /FLUXCB/ CLDMC(LMA,2*IM,2*JM-2),CLDSS(LMA,2*IM,2*JM-2) 23. COMMON /WORK00/ PKDN(IM,JM,LMA),PKUP(IM,JM,LMA),MSUM(IM,JM) 24. COMMON /WORK01/ VMF(LMA,LMA), LMOFL1(LMA),LMCMAX 25. COMMON /WORK02/ UAT(IM,JM,LMA),VAT(IM,JM,LMA), 26. * MG(LMA,2,2), PR(LMA,2,2),PK(LMA,2,2), 27. * H0MQ(LMA,2,2),HZMQ(LMA,2,2), Q0MQ(LMA,2,2),QZMQ(LMA,2,2), 28. * UV(IM+1,LMA),UM(IM+1,LMA),UMS(IM+1,LMA), 28.1 * RA(IM+1),RAQ(4), ID(IM+1) 29. C**** 30. C**** PRECCB PREC Precipitation from atmosphere (kg/m^2) 31. C**** EPRE Temperature of precipitation (C) 32. C**** 33. C**** FLUXCB CLDMC Moist convective cloud precipitation (kg/m^2) 34. C**** 35. C**** Calculate PK = P**RKAP, integrate pressures from the top down 36. IF(QPK) GO TO 40 37. DO 10 I=IM,IM*(JM-1)+1 38. MSUM(I,1) = MSTRAT 39. DO 10 L=LMA,1,-1 40. MSUM(I,1) = MSUM(I,1) + MA(I,1,L) 41. PKDN(I,1,L) = (GRAV*(MSUM(I,1)-MA(I,1,L)*.25))**RKAP 42. 10 PKUP(I,1,L) = (GRAV*(MSUM(I,1)-MA(I,1,L)*.75))**RKAP 43. DO 20 I=1,IM 44. MSUM(I, 1) = MSUM(IM,1) 45. MSUM(I,JM) = MSUM(1,JM) 46. DO 20 L=1,LMA 47. PKDN(I, 1,L) = PKDN(IM,1,L) 48. PKDN(I,JM,L) = PKDN(1,JM,L) 49. PKUP(I, 1,L) = PKUP(IM,1,L) 50. 20 PKUP(I,JM,L) = PKUP(1,JM,L) 51. QPK = .true. 52. C**** Save UAT and VAT, and zero out PREC, CLDSS and CLDMC 53. 40 DO 50 I=1,IM*JM*LMA 54. UAT(I,1,1) = UA(I,1,1) 55. 50 VAT(I,1,1) = VA(I,1,1) 55.5 DO 60 I=1,IM*JM 56. 60 PREC(I,1,1) = 0. 57. DO 70 L=1,LMA*IM*2*(JM*2-2) 58. CLDMC(L,1,1) = 0. 59. 70 CLDSS(L,1,1) = 0. 60. C**** Prevent negative water vapor mass in a quarter box 61. DO 80 L=1,LMA 62. DO 80 J=2,JM-1 63. DO 80 I=1,IM 64. QXY = RRT3*(ABS(QXM(I,J,L))+ABS(QYM(I,J,L))) + 1.D-20 65. IF(QXY.le.Q0M(I,J,L)) GO TO 80 66. QXM(I,J,L) = QXM(I,J,L)*(Q0M(I,J,L)/QXY) 67. QYM(I,J,L) = QYM(I,J,L)*(Q0M(I,J,L)/QXY) 68. 80 CONTINUE 69. C**** 70. C**** Loop over J 71. C**** 72. DO 850 J=1,JM 73. QPOLE = J.eq.1 .or. J.eq.JM 74. IF(.not.QPOLE) GO TO 200 75. C**** 76. C**** Define variables at the poles 77. C**** 78. I = 1 79. IQ = 1 80. I2 = 1 81. JQ = 2 82. IF(J.eq.JM) JQ = 1 83. J2 = 2*J+JQ-3 84. C**** Grid location parameters needed for momentum mixing 85. KMAX = IM+1 86. JVPO = J+JQ-2 87. IF(J.eq.JM) JVPO=JM-1 88. DO 110 K=1,IM 89. RA(K) = RAMVN(1) 90. 110 ID(K) = K + (JVPO-1)*IM + IM*JM*LMA 91. ID(IM+1) = 1 + (J-1)*IM 92. RA(IM+1) = 1. 93. C**** Gaseous air mass (kg/m^2), pressure (Pa) and exner function 94. MS = MSTRAT 95. DO 130 L=LMA,1,-1 96. MG(L,IQ,JQ) = MA(I,J,L) 97. PR(L,IQ,JQ) = (MS+.5*MA(I,J,L))*GRAV 98. PK(L,IQ,JQ) = (PKDN(I,J,L)+PKUP(I,J,L))*.5 99. MS = MS + MA(I,J,L) 100. C**** Potential enthalpy of a quarter box (J/m^2) 101. H0MQ(L,IQ,JQ) = H0M(I,J,L)*BYDXYP(J) 102. HZMQ(L,IQ,JQ) = HZM(I,J,L)*BYDXYP(J) 103. C**** Water vapor of a quarter box (kg/m^2) 104. Q0MQ(L,IQ,JQ) = Q0M(I,J,L)*BYDXYP(J) 105. QZMQ(L,IQ,JQ) = QZM(I,J,L)*BYDXYP(J) 106. C**** Momentum (kg/m*s) 107. VMF(L,L)= 0. 108. DO 120 K=1,KMAX 109. UV(K,L) = UAT(ID(K),1,L) 110. 120 UM(K,L) = UAT(ID(K),1,L)*MG(L,IQ,JQ) 111. 130 continue 112. EPRE(I,J,1) = H0M(I,J,1)*PK(1,IQ,JQ)/(SHCD*DXYP(J)*MA(I,J,1))-TKF 113. LMCMAX=0 114. GO TO 500 115. C**** 116. C**** Loop over I 117. C**** Define variables away from the poles 118. C**** 118.1 200 KMAX = 4 118.2 RA(1) = .5 118.3 RA(2) = .5 118.4 RA(3) = RAMVS(J) 118.5 RA(4) = RAMVN(J) 119. IM1 = IM 120. I = 1 122. 220 ID(1) = IM1 + (J-1)*IM 123. ID(2) = I + (J-1)*IM 124. ID(4) = ID(2) + IM*JM*LMA 125. ID(3) = ID(4) - IM 126. C**** Gaseous air mass (kg/m^2), pressure (Pa) and exner function 127. MS = MSTRAT 128. DO 270 L=LMA,1,-1 129. DO 260 IQJQ=1,4 130. MG(L,IQJQ,1) = MA(I,J,L) 131. PR(L,IQJQ,1) = (MS+.5*MA(I,J,L))*GRAV 132. 260 PK(L,IQJQ,1) = (PKDN(I,J,L)+PKUP(I,J,L))*.5 133. MS = MS + MA(I,J,L) 134. DO 270 K=1,KMAX 135. 270 UMS(K,L)=0. 136. EPRE(I,J,1) = H0M(I,J,1)*PK(1,1,1)/(SHCD*DXYP(J)*MA(I,J,1)) - TKF 137. LMCMAX=0 138. C**** 139. C**** Loop over horizontal quarter boxes in north-south direction 140. C**** 141. JQ = 1 142. 300 J2 = 2*J+JQ-3 143. RJ = 2.*RRT3*JQ - 3.*RRT3 144. RAQ(3) = .625 - .25*JQ 145. RAQ(4) = -.125 + .25*JQ 146. C**** 147. C**** Loop over horizontal quarter boxes in east-west direction 148. C**** 149. IQ = 1 150. 400 I2 = 2*I+IQ-2 151. RI = 2.*RRT3*IQ - 3.*RRT3 152. RAQ(1) = .625 - .25*IQ 153. RAQ(2) = -.125 + .25*IQ 154. C**** Potential enthalpy of a quarter box (J/m^2) 155. DO 420 L=LMA,1,-1 156. H0MQ(L,IQ,JQ) = BYDXYP(J)*(H0M(I,J,L)+RI*HXM(I,J,L)+RJ*HYM(I,J,L)) 157. HZMQ(L,IQ,JQ) = BYDXYP(J)* HZM(I,J,L) 158. C**** Water vapor of a quarter box (kg/m^2) 159. Q0MQ(L,IQ,JQ) = BYDXYP(J)*(Q0M(I,J,L)+RI*QXM(I,J,L)+RJ*QYM(I,J,L)) 160. QZMQ(L,IQ,JQ) = BYDXYP(J)* QZM(I,J,L) 161. C**** Momentum (kg/m*s) 162. VMF(L,L)= 0. 163. DO 410 K=1,KMAX 164. UV(K,L) = UAT(ID(K),1,L) 165. 410 UM(K,L) = UAT(ID(K),1,L)*MG(L,IQ,JQ) 166. 420 continue 167. C**** 168. C**** Call condensation subroutines with one-dimensional arrays 169. C**** 170. 500 CALL MSTCN2 (MG(1,IQ,JQ), PR(1,IQ,JQ),PK(1,IQ,JQ), 171. * H0MQ(1,IQ,JQ),HZMQ(1,IQ,JQ), Q0MQ(1,IQ,JQ),QZMQ(1,IQ,JQ), 172. * PREC(I,J,1), CLDMC(1,I2,J2), 173. * AIJ(I,J,71),AIJ(I,J,72),AJL(J,1,19)) 174. C CALL LSCOND 175. C**** 176. C**** Vertical mixing of momentum 177. C**** 178. DO 540 LMIN=1,LMCM 179. IF(VMF(LMIN,LMIN).le.0.) GO TO 540 180. LMAX=LMOFL1(LMIN) 181. DO 530 K=1,KMAX 182. UM(K,LMIN) = UM(K,LMIN) + (UV(K,LMIN+1)-UV(K,LMIN))*VMF(LMIN,LMIN) 183. DO 510 L=LMIN+1,LMAX-1 184. 510 UM(K,L) = UM(K,L) + UV(K,LMIN)*(VMF(LMIN,L-1)-VMF(LMIN,L)) 185. * - UV(K,L)*VMF(LMIN,L-1) + UV(K,L+1)*VMF(LMIN,L) 186. 520 UM(K,LMAX) = UM(K,LMAX) + (UV(K,LMIN)-UV(K,LMAX))*VMF(LMIN,LMAX-1) 187. DO 530 L=LMIN,LMAX 188. 530 UV(K,L) = UM(K,L)/MG(L,IQ,JQ) 189. 540 continue 190. C**** Accumulate momentum changes over quarter boxes 191. IF(QPOLE) GO TO 800 192. DO 550 L=1,LMA 193. DO 550 K=1,KMAX 194. 550 UMS(K,L) = UMS(K,L) + UM(K,L)*RAQ(K) 195. C**** 196. C**** End of loop over IQ 197. C**** 198. IQ = IQ + 1 199. IF(IQ.le.2) GO TO 400 200. C**** 201. C**** End of loop over JQ 202. C**** 203. JQ = JQ + 1 204. IF(JQ.le.2) GO TO 300 205. C**** 206. C**** Update model prognostic variables 207. C**** 208. IF(LMCMAX.eq.0) GO TO 750 209. PREC(I,J,1) = PREC(I,J,1)*.25 210. DO 720 L=1,LMCMAX 211. AJL(J,L,21) = AJL(J,L,21) - Q0M(I,J,L) + .25*DXYP(J)* 212. * (Q0MQ(L,2,2)+Q0MQ(L,1,2)+Q0MQ(L,2,1)+Q0MQ(L,1,1)) 213. c AJL(J,L,22) = AJL(J,L,22) - PLK*(H0M(I,J,L) - .25*DXYP(J)* 214. c * (H0MQ(L,2,2)+H0MQ(L,1,2)+H0MQ(L,2,1)+H0MQ(L,1,1))) 215. H0M(I,J,L) = .25*DXYP(J)* 216. * (H0MQ(L,2,2)+H0MQ(L,1,2)+H0MQ(L,2,1)+H0MQ(L,1,1)) 217. HXM(I,J,L) = (.25/RRT3)*DXYP(J)* 218. * (H0MQ(L,2,2)-H0MQ(L,1,2)+H0MQ(L,2,1)-H0MQ(L,1,1)) 219. HYM(I,J,L) = (.25/RRT3)*DXYP(J)* 220. * (H0MQ(L,2,2)+H0MQ(L,1,2)-H0MQ(L,2,1)-H0MQ(L,1,1)) 221. HZM(I,J,L) = .25*DXYP(J)* 222. * (HZMQ(L,2,2)+HZMQ(L,1,2)+HZMQ(L,2,1)+HZMQ(L,1,1)) 223. Q0M(I,J,L) = .25*DXYP(J)* 224. * (Q0MQ(L,2,2)+Q0MQ(L,1,2)+Q0MQ(L,2,1)+Q0MQ(L,1,1)) 225. QXM(I,J,L) = (.25/RRT3)*DXYP(J)* 226. * (Q0MQ(L,2,2)-Q0MQ(L,1,2)+Q0MQ(L,2,1)-Q0MQ(L,1,1)) 227. QYM(I,J,L) = (.25/RRT3)*DXYP(J)* 228. * (Q0MQ(L,2,2)+Q0MQ(L,1,2)-Q0MQ(L,2,1)-Q0MQ(L,1,1)) 229. QZM(I,J,L) = .25*DXYP(J)* 230. * (QZMQ(L,2,2)+QZMQ(L,1,2)+QZMQ(L,2,1)+QZMQ(L,1,1)) 231. DO 710 K=1,KMAX 232. 710 UA(ID(K),1,L) = UA(ID(K),1,L) + 233. * (UMS(K,L)/MA(I,J,L) - UAT(ID(K),1,L))*RA(K) 234. 720 continue 235. AIJ(I,J,15) = AIJ(I,J,15) + PREC(I,J,1) 236. C**** 237. C**** End of loop over I 238. C**** 239. 750 IM1 = I 240. I = I+1 241. IF(I.le.IM) GO TO 220 242. GO TO 850 243. C**** 244. C**** Update model prognostic variables at the poles 245. C**** 246. 800 IF(LMCMAX.eq.0) GO TO 850 247. DO 820 L=1,LMCMAX 248. AJL(J,L,21) = AJL(J,L,21) - Q0M(I,J,L) + Q0MQ(L,IQ,JQ)*DXYP(J) 249. c AJL(J,L,22) = AJL(J,L,22) - 250. c - PLK*(H0M(I,J,L) - H0MQ(L,IQ,JQ)*DXYP(J)) 251. H0M(I,J,L) = H0MQ(L,IQ,JQ)*DXYP(J) 252. HZM(I,J,L) = HZMQ(L,IQ,JQ)*DXYP(J) 253. Q0M(I,J,L) = Q0MQ(L,IQ,JQ)*DXYP(J) 254. QZM(I,J,L) = QZMQ(L,IQ,JQ)*DXYP(J) 255. DO 810 K=1,KMAX 256. 810 UA(ID(K),1,L) = UA(ID(K),1,L) + 257. * (UM(K,L)/MA(1,J,L) - UAT(ID(K),1,L))*RA(K) 258. 820 continue 259. AIJ(I,J,15) = AIJ(I,J,15) + PREC(I,J,1) 260. C**** 261. C**** End of loop over J 262. C**** 263. 850 continue 264. C**** Diagnostic for change of angular momentum by moist convection 265. DO 860 J=1,JM 266. IMAX=IM 267. IF(J.eq.1 .or. J.eq.JM) IMAX=1 268. I=IM 269. DO 860 L=1,LMA 270. DO 860 IP1=1,IMAX 271. AJL(J,L,24) = AJL(J,L,24) + 272. * (UA(I,J,L)-UAT(I,J,L))*(MA(I,J,L)+MA(IP1,J,L)) 273. 860 I=IP1 274. RETURN 275. END 500. 501. SUBROUTINE MSTCN2 (MAQ, PL,PLK, H0MQ,HZMQ, Q0MQ,QZMQ, 502. * PREC, CLDMC, AIJ71,AIJ72,AJL19) 503. C**** 504. C**** MSTCNV receives one-dimensional arrays of mass and momentum and 505. C**** modifies these arrays based on the Moist Convective algorithm. 506. C**** The linear upstream scheme is used for the downward subsidence. 507. C**** Stored vertical mass fluxes will be applied to momentum later. 508. C**** 509. IMPLICIT REAL*8 (A-H,M,O-Z) 510. PARAMETER (LMA=9,LMCM=7, 511. * RDRY=287.,RVAP=461.5, SHCD=1003.5, RKAP=RDRY/SHCD, 512. * grav=9.81d0, ELHE=2500000.,ELHS=2834000., 513. * FPLUME=.25, TKF=273.16d0, TI=233.16d0) 514. REAL*8 MAQ(LMA), PL(LMA),PLK(LMA), H0MQ(LMA),HZMQ(LMA), 515. * Q0MQ(LMA),QZMQ(LMA), CLDMC(LMA), AJL19(LMA) 517. COMMON /WORK01/ VMF(LMA,LMA), LMOFL1(LMA),LMCMAX 518. C**** Internal working arrays 519. REAL*8 MAP(LMA),H0MP(LMA),Q0MP(LMA), C(LMA),CM(LMA), 520. * CDHEAT(LMA),COND(LMA), HMOLD(LMA),QMOLD(LMA), 521. * FH0(0:LMA),FHZ(LMA), FQ0(0:LMA),FQZ(LMA) 522. C**** QSAT = (RAIR/RVAP)*610.71*EXP((L/RVAP)*(1/TKF-1/T)) / PRESS 523. QSAT(TM,PR) = 379.7915d0*EXP(ELHX*(1./TKF-1./TM)/RVAP)/PR 524. C**** 525. C**** Zero out CDHEAT,COND. Copy HMOLD,QMOLD from H0MQ,Q0MQ. 526. C**** 527. DO 230 L=1,LMA 528. CDHEAT(L) = 0. 529. COND(L) = 0. 530. HMOLD(L) = H0MQ(L) 531. 230 QMOLD(L) = Q0MQ(L) 532. C**** 533. C**** Loop over convective base layer LMIN 534. C**** 535. LMIN=0 536. 300 LMIN=LMIN+1 537. IF(LMIN.ge.LMCM) GO TO 600 538. C**** 539. C**** Test whether plume should be created from layer LMIN 540. C**** Assumes MAQ(LMIN)*FPLUME < MAQ(LMIN+1)*.95 541. C**** Assumes MAQ(LMIN)*FPLUME > MAQ(LMIN+1)*.001 542. C**** 543. C**** Test for moist static stability and energy 544. HDN = H0MQ(LMIN )/MAQ(LMIN ) 545. HUP = H0MQ(LMIN+1)/MAQ(LMIN+1) 546. QDN = Q0MQ(LMIN )/MAQ(LMIN ) 547. QUP = Q0MQ(LMIN+1)/MAQ(LMIN+1) 548. IF(PLK(LMIN)*(HUP-HDN)+ELHE*(QUP-QDN).ge.0.) GO TO 300 549. HEDGE= THBAR(HUP/SHCD,HDN/SHCD)*SHCD 550. ELHX = ELHE 551. DMSE = (HUP-HEDGE)*PLK(LMIN+1) + (HEDGE-HDN)*PLK(LMIN) + 552. * ELHE*(QSAT(HUP*PLK(LMIN+1)/SHCD,PL(LMIN+1))-QDN) 553. IF(DMSE.gt.-1.D-10*SHCD) GO TO 300 554. C**** Test for condensation 555. MAP(LMIN+1) = MAQ(LMIN)*FPLUME 556. H0MP(LMIN+1) = H0MQ(LMIN)*FPLUME 557. Q0MP(LMIN+1) = Q0MQ(LMIN)*FPLUME 558. C(LMIN) = -MAP(LMIN+1)/MAQ(LMIN+1) 559. FH0(LMIN) = C(LMIN)*(H0MQ(LMIN+1)-(1.+C(LMIN))*HZMQ(LMIN+1)) 560. WORK = (FH0(LMIN)+H0MP(LMIN+1))*(PLK(LMIN+1)-PLK(LMIN)) 561. H0MP(LMIN+1) = H0MP(LMIN+1) - WORK/PLK(LMIN+1) 562. TP = H0MP(LMIN+1)*PLK(LMIN+1)/(MAP(LMIN+1)*SHCD) 563. IF(TP.lt.TI) ELHX = ELHS 564. QSATMP = MAP(LMIN+1)*QSAT(TP,PL(LMIN+1)) 565. IF(Q0MP(LMIN+1).lt.QSATMP) GO TO 300 566. C**** Plume is established. Remove it from layer LMIN. 567. AJL19(LMIN) = AJL19(LMIN) + MAP(LMIN+1) 568. CM(LMIN) = -MAP(LMIN+1) 569. MAQ(LMIN) = MAQ(LMIN)*(1.-FPLUME) 570. H0MQ(LMIN) = H0MQ(LMIN)*(1.-FPLUME) 571. HZMQ(LMIN) = HZMQ(LMIN)*(1.-FPLUME) 572. Q0MQ(LMIN) = Q0MQ(LMIN)*(1.-FPLUME) 573. QZMQ(LMIN) = QZMQ(LMIN)*(1.-FPLUME) 574. C DO 310 K=1,KMAX 575. C UMP(K,LMIN+1) = UM(K,LMIN)*FPLUME 576. C 310 UM(K,LMIN) = UM(K,LMIN)*(1.-FPLUME) 577. VMF(LMIN,LMIN) = MAP(LMIN+1) 578. C**** Condense water vapor in the plume, releasing latent heat 579. IF(TP.ge.TKF .or. TP.lt.TI) GO TO 320 580. ELHX = ELHS 581. QSATMP = MAP(LMIN+1)*QSAT(TP,PL(LMIN+1)) 582. 320 GAMA = ELHX*ELHX*QSATMP / (SHCD*RVAP*TP*TP*MAP(LMIN+1)) 583. COND(LMIN+1) = (Q0MP(LMIN+1)-QSATMP) / (1.+GAMA) 584. CDHEAT(LMIN+1) = COND(LMIN+1)*ELHX 585. H0MP(LMIN+1) = H0MP(LMIN+1) + COND(LMIN+1)*ELHX/PLK(LMIN+1) 586. Q0MP(LMIN+1) = Q0MP(LMIN+1) - COND(LMIN+1) 587. C**** 588. C**** Test whether plume rises into successively higher layers 589. C**** 590. DO 380 L=LMIN+2,LMA 591. C**** Test whether plume has sufficient mass 592. IF(MAP(L-1).le..001d0*MAQ(L)) GO TO 400 593. C**** Test for static stability and energy 594. HDN = H0MP(L-1)/MAP(L-1) 595. HUP = H0MQ(L )/MAQ(L) 596. QDN = Q0MP(L-1)/MAP(L-1) 597. QUP = Q0MQ(L )/MAQ(L) 598. IF(PLK(L-1)*(HUP-HDN)+ELHE*(QUP-QDN).ge.0.) GO TO 400 599. HEDGE = THBAR(HUP/SHCD,HDN/SHCD)*SHCD 600. ELHX = ELHE 601. DMSE = (HUP-HEDGE)*PLK(L) + (HEDGE-HDN)*PLK(L-1) + 602. * ELHE*(QSAT(HUP*PLK(L)/SHCD,PL(L))-QDN) 603. IF(DMSE.gt.-1.D-10*SHCD) GO TO 400 604. C**** Test for condensation 605. FPRISE = 1. 606. IF(MAP(L-1).gt..95d0*MAQ(L)) FPRISE = .95d0*MAQ(L)/MAP(L-1) 607. MAP(L) = MAP(L-1)*FPRISE 608. H0MP(L) = H0MP(L-1)*FPRISE 609. Q0MP(L) = Q0MP(L-1)*FPRISE 610. C(L-1) = -MAP(L)/MAQ(L) 611. FH0(L-1)= C(L-1)*(H0MQ(L)-(1.+C(L-1))*HZMQ(L)) 612. WORK = (FH0(L-1)+H0MP(L))*(PLK(L)-PLK(L-1)) 613. H0MP(L) = H0MP(L) - WORK/PLK(L) 614. TP = H0MP(L)*PLK(L)/(MAP(L)*SHCD) 615. IF(TP.lt.TI) ELHX = ELHS 616. QSATMP = MAP(L)*QSAT(TP,PL(L)) 617. IF(Q0MP(L).lt.QSATMP) GO TO 400 618. C**** All or part of the plume will rise. Remove it from layer L-1. 619. AJL19(L-1) = AJL19(L-1) + MAP(L) 620. CM(L-1) = -MAP(L) 621. MAP(L-1) = MAP(L-1)*(1.-FPRISE) 622. H0MP(L-1) = H0MP(L-1)*(1.-FPRISE) 623. Q0MP(L-1) = Q0MP(L-1)*(1.-FPRISE) 624. C DO 360 K=1,KMAX 625. C UMP(K,L ) = UMP(K,L-1)*FPRISE 626. C 360 UMP(K,L-1) = UMP(K,L-1)*(1.-FPRISE) 627. VMF(LMIN,L-1) = MAP(L) 628. C**** Condense water vapor in the plume, releasing latent heat 629. IF(TP.ge.TKF .or. TP.lt.TI) GO TO 370 630. ELHX = ELHS 631. QSATMP = MAP(L)*QSAT(TP,PL(L)) 632. 370 GAMA = ELHX*ELHX*QSATMP / (SHCD*RVAP*TP*TP*MAP(L)) 633. COND(L) = (Q0MP(L)-QSATMP) / (1.+GAMA) 634. CDHEAT(L) = COND(L)*ELHX 635. H0MP(L) = H0MP(L) + COND(L)*ELHX/PLK(L) 636. 380 Q0MP(L) = Q0MP(L) - COND(L) 637. C**** 638. C**** Subsidense of momentum uses the upstream scheme 639. C**** 640. 400 LMAX = L-1 641. LMOFL1(LMIN) = LMAX 642. C DO 420 K=1,KMAX 643. C UM(K,LMIN) = UM(K,LMIN) - C(LMIN)*UM(K,LMIN+1) 644. C DO 410 L=LMIN+1,LMAX-1 645. C 410 UM(K,L) = UM(K,L) + C(L-1)*UM(K,L) - C(L)*UM(K,L+1) + UMP(K,L) 646. C 420 UM(K,LMAX) = UM(K,LMAX) + C(LMAX-1)*UM(K,LMAX) + UMP(K,LMAX) 647. C**** 648. C**** Subsidence of potential enthalpy and water vapor uses the linear 649. C**** upstream scheme. The plume is removed from layer LMAX before 650. C**** subsidence and the plume is added back to the higher layers after 651. C**** subsidence. Water vapor is restricted to be nonnegative. 652. C**** 653. DO 450 L=LMIN,LMAX-1 654. C**** Calculate FQ0 (kg) and FQZ (kg**2), air mass flux is negative 655. C FH0(L) = C(L)*(H0MQ(L+1)-(1.+C(L))*HZMQ(L+1)) 656. FHZ(L) = CM(L)*(C(L)*C(L)*HZMQ(L+1)-3.*FH0(L)) 657. FQ0(L) = C(L)*(Q0MQ(L+1)-(1.+C(L))*QZMQ(L+1)) 658. 450 FQZ(L) = CM(L)*(C(L)*C(L)*QZMQ(L+1)-3.*FQ0(L)) 659. C**** 660. C**** Modify the specific humidity moments so that the specific 661. C**** humidity mass in each division is non-negative 662. C**** 663. DO 470 L=LMIN+1,LMAX 664. C**** Air is leaving through the bottom edge: 2 or 3 divisions 665. IF(FQ0(L-1).le.0.) GO TO 460 666. C**** Bottom most division is negative, QMB = -FQ0(L-1) < 0: Case 2 or 4 667. QZMQ(L) = Q0MQ(L)/(1.+C(L-1)) 668. FQ0(L-1) = 0. 669. FQZ(L-1) = CM(L-1)*C(L-1)*C(L-1)*QZMQ(L) 670. GO TO 470 671. C**** No air is leaving through the top edge: 2 divisions 672. 460 IF(Q0MQ(L)+FQ0(L-1).ge.0.) GO TO 470 673. C**** Top most division is negative, QMT = Q0MQ(L)+FQ0(L-1) < 0: Case 3 674. QZMQ(L) = Q0MQ(L)/C(L-1) 675. FQ0(L-1) = -Q0MQ(L) 676. FQZ(L-1) = CM(L-1)*(C(L-1)+3.)*Q0MQ(L) 677. 470 CONTINUE 678. C**** 679. C**** Calculate new tracer mass and first moments of tracer mass 680. C**** 681. C**** Calculation in the first layer 682. H0MQ(LMIN) = H0MQ(LMIN) - FH0(LMIN) 683. HZMQ(LMIN) = (HZMQ(LMIN)*MAQ(LMIN) - FHZ(LMIN) 684. * + 3.*(CM(LMIN)*H0MQ(LMIN) - MAQ(LMIN)*FH0(LMIN))) 685. * / (MAQ(LMIN)-CM(LMIN)) 686. Q0MQ(LMIN) = Q0MQ(LMIN) - FQ0(LMIN) 687. QZMQ(LMIN) = (QZMQ(LMIN)*MAQ(LMIN) - FQZ(LMIN) 688. * + 3.*(CM(LMIN)*Q0MQ(LMIN) - MAQ(LMIN)*FQ0(LMIN))) 689. * / (MAQ(LMIN)-CM(LMIN)) 690. MAQ(LMIN) = MAQ(LMIN) - CM(LMIN) 691. C MAQ(LMIN) = MA(I,J,LMIN) 692. C**** Calculation in the interior layers 693. DO 480 L=LMIN+1,LMAX-1 694. H0MQ(L) = H0MQ(L) + (FH0(L-1)-FH0(L)) 695. HZMQ(L) = (HZMQ(L)*MAQ(L) + (FHZ(L-1)-FHZ(L)) 696. * + 3.*((CM(L-1)+CM(L))*H0MQ(L) 697. * - MAQ(L)*(FH0(L-1)+FH0(L)))) / (MAQ(L)+(CM(L-1)-CM(L))) 698. Q0MQ(L) = Q0MQ(L) + (FQ0(L-1)-FQ0(L)) 699. QZMQ(L) = (QZMQ(L)*MAQ(L) + (FQZ(L-1)-FQZ(L)) 700. * + 3.*((CM(L-1)+CM(L))*Q0MQ(L) 701. * - MAQ(L)*(FQ0(L-1)+FQ0(L)))) / (MAQ(L)+(CM(L-1)-CM(L))) 702. MAQ(L) = MAQ(L) + (CM(L-1)-CM(L)) 703. 480 CONTINUE 704. C**** Calculation in the top layer 705. H0MQ(LMAX) = H0MQ(LMAX) + FH0(LMAX-1) 706. HZMQ(LMAX) = (HZMQ(LMAX)*MAQ(LMAX) + FHZ(LMAX-1) 707. * + 3.*(CM(LMAX-1)*H0MQ(LMAX) - MAQ(LMAX)*FH0(LMAX-1))) 708. * / (MAQ(LMAX)+CM(LMAX-1)) 709. Q0MQ(LMAX) = Q0MQ(LMAX) + FQ0(LMAX-1) 710. QZMQ(LMAX) = (QZMQ(LMAX)*MAQ(LMAX) + FQZ(LMAX-1) 711. * + 3.*(CM(LMAX-1)*Q0MQ(LMAX) - MAQ(LMAX)*FQ0(LMAX-1))) 712. * / (MAQ(LMAX)+CM(LMAX-1)) 713. MAQ(LMAX) = MAQ(LMAX) + CM(LMAX-1) 714. C**** Add changes by plume in upper layers after subsidence 715. DO 490 L=LMIN+1,LMAX 716. H0MQ(L) = H0MQ(L) + H0MP(L) 717. Q0MQ(L) = Q0MQ(L) + Q0MP(L) 718. 490 MAQ(L) = MAQ(L) + MAP(L) 719. C 490 MAQ(L) = MA(I,J,L) 720. FH0(LMIN-1) = 0. 721. FH0(LMAX) = 0. 722. FQ0(LMIN-1) = 0. 723. FQ0(LMAX) = 0. 724. C**** 725. C**** REEVAPORATION AND PRECIPITATION 726. C**** 727. PREQ = COND(LMAX) 728. PRHEAT = CDHEAT(LMAX) 729. DO 540 L=LMAX-1,1,-1 730. IF(PREQ.le.0.) GO TO 530 731. CLDMC(L+1) = CLDMC(L+1) + PREQ 732. C**** REEVAPORATE ALL PRECIPITATION FROM ABOVE 733. EVAP = PREQ 734. PREQ = 0. 735. C**** FORWARD STEP COMPUTES HUMIDITY CHANGE BY RECONDENSATION 736. C**** Q = Q + F(TOLD,PRHEAT,QOLD+EVAP) 737. MCLOUD = .25*MAQ(L) 738. IF(L.le.LMIN) MCLOUD = .5*MAQ(L) 739. TOLD = HMOLD(L)*PLK(L)/(MAQ(L)*SHCD) 740. TN = TOLD - PRHEAT/(MCLOUD*SHCD) 741. QN = QMOLD(L)/MAQ(L)+EVAP/MCLOUD 742. ELHX = ELHE 743. IF(TOLD.lt.TKF) ELHX = ELHS 744. QSATC = QSAT(TN,PL(L)) 745. IF(QN-QSATC.le.0.) GO TO 520 746. DO 510 N=1,3 747. IF(N.ne.1) QSATC = QSAT(TN,PL(L)) 748. DQ = (QN-QSATC)/(1.+ELHX*ELHX*QSATC/(SHCD*RVAP*TN*TN)) 749. TN = TN + ELHX*DQ/SHCD 750. QN = QN - DQ 751. 510 PREQ = PREQ + DQ 752. PREQ = MCLOUD*PREQ 753. IF(PREQ.gt.EVAP) PREQ = EVAP 754. C**** UPDATE TEMPERATURE AND HUMIDITY DUE TO NET REVAPORATION IN CLOUD 755. 520 H0MQ(L) = H0MQ(L) + (ELHX*PREQ-PRHEAT)/PLK(L) 756. Q0MQ(L) = Q0MQ(L) + EVAP-PREQ 757. C**** ADD PRECIPITATION AND LATENT HEAT BELOW 758. 530 PRHEAT = CDHEAT(L) + ELHX*PREQ 759. PREQ = PREQ + COND(L) 760. C**** 761. 540 CONTINUE 762. C**** 763. CLDMC(1) = CLDMC(1) + PREQ 764. PREC = PREC + PREQ 765. C**** END OF CONVECTION CYCLE. SET SOME VARIABLES TO DEFAULTS. 766. DO 560 L=1,LMAX 767. COND(L) = 0. 768. CDHEAT(L) = 0. 769. 560 CONTINUE 770. DO 570 L=1,LMAX-1 771. HMOLD(L) = H0MQ(L) 772. 570 QMOLD(L) = Q0MQ(L) 773. C**** 774. C**** END OF OUTER LOOP OVER CLOUD BASE 775. C**** 776. IF(PL(LMAX).ge.70000.) AIJ71 = AIJ71 + 1. 777. IF(PL(LMIN)-PL(LMAX)+.5*(MAQ(LMIN)+MAQ(LMAX))*GRAV .ge. 45000.) 778. * AIJ72 = AIJ72 + 1. ! GRAV not defined 779. IF(LMCMAX.lt.LMAX) LMCMAX=LMAX 780. GO TO 300 781. 600 RETURN 782. END 1000. 1001. SUBROUTINE CONDSE 1002. C**** 1003. C**** CONDSE calculates large scale condensation 1004. C**** 1005. INCLUDE 'C070.COM' 1006. LOGICAL*4 QPOLE 1007. PARAMETER (AGEMAX=50.) 1008. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 1009. COMMON /FLUXCB/ CLDMC(LMA,2*IM,2*JM-2),CLDSS(LMA,2*IM,2*JM-2) 1010. COMMON /WORK00/ PKDN(IM,JM,LMA),PKUP(IM,JM,LMA) 1011. COMMON /WORK01/ LSCMAX 1012. COMMON /WORK02/ MG(LMA,2,2),PR(2*LMA,2,2), 1013. * TK(2*LMA,2,2),QQ(2*LMA,2,2) 1014. C**** 1015. C**** PRECCB PREC Precipitation from atmosphere (kg/m^2) 1016. C**** EPRE Energy of precipitation (J/m^2) 1017. C**** 1018. C**** FLUXCB CLDSS Super saturation cloud precipitation (kg/m^2) 1019. C**** 1020. C**** SURFCB BLDATA(4) Age of snow (days) 1021. C**** 1022. DO 10 I=1,IM*JM 1023. 10 PREC(I,1,2) = 0. 1024. IF(.not.QRAD) GO TO 900 1025. C**** Zero out large scale condensation clouds 1026. DO 20 L=1,LMA*2*IM*(2*JM-2) 1027. 20 CLDSS(L,1,1) = 0. 1028. C**** 1029. C**** Loop over J 1030. C**** 1031. DO 850 J=1,JM 1032. QPOLE = J.eq.1 .or. J.eq.JM 1033. IF(.not.QPOLE) GO TO 200 1034. C**** 1035. C**** Define variables at the poles 1036. C**** 1037. I = 1 1038. IQ = 1 1039. I2 = 1 1040. JQ = 2 1041. IF(J.eq.JM) JQ = 1 1042. J2 = 2*J+JQ-3 1043. C**** Gaseous air mass (kg/m^2), pressure (Pa) and exner function 1044. MS = MSTRAT 1045. DO 110 L=LMA,1,-1 1046. MG(L,IQ,JQ) = MA(I,J,L) 1047. PR(2*L ,IQ,JQ) = (MS+.25*MA(I,J,L))*GRAV 1048. PR(2*L-1,IQ,JQ) = (MS+.75*MA(I,J,L))*GRAV 1049. C PK(2*L ,IQ,JQ) = PKUP(I,J,L) 1050. C PK(2*L-1,IQ,JQ) = PKDN(I,J,L) 1051. MS = MS + MA(I,J,L) 1052. C**** Temperature (K) at poles 1053. TK(2*L ,IQ,JQ) = (H0M(I,J,L)+.5*HZM(I,J,L))* 1054. * PKUP(I,J,L) / (MA(I,J,L)*DXYP(J)*SHCD) 1055. TK(2*L-1,IQ,JQ) = (H0M(I,J,L)-.5*HZM(I,J,L))* 1056. * PKDN(I,J,L) / (MA(I,J,L)*DXYP(J)*SHCD) 1057. C**** Specific humidity at poles 1058. QQ(2*L ,IQ,JQ) = (Q0M(I,J,L)+.5*QZM(I,J,L)) / (MA(I,J,L)*DXYP(J)) 1059. 110 QQ(2*L-1,IQ,JQ) = (Q0M(I,J,L)-.5*QZM(I,J,L)) / (MA(I,J,L)*DXYP(J)) 1060. C LSCMAX = 0 1061. GO TO 500 1062. C**** 1063. C**** Loop over I 1064. C**** Define variables away from the poles 1065. C**** 1066. 200 I = 1 1067. C**** Gaseous air mass (kg/m^2), pressure (Pa) and exner function 1068. 250 MS = MSTRAT 1069. DO 270 L=LMA,1,-1 1070. DO 260 IQJQ=1,4 1071. MG(L,IQJQ,1) = MA(I,J,L) 1072. PR(2*L ,IQJQ,1) = (MS+.25*MA(I,J,L))*GRAV 1073. 260 PR(2*L-1,IQJQ,1) = (MS+.75*MA(I,J,L))*GRAV 1074. C PK(2*L ,IQJQ,1) = PKUP(I,J,L) 1075. C 260 PK(2*L-1,IQJQ,1) = PKDN(I,J,L) 1076. 270 MS = MS + MA(I,J,L) 1077. C LSCMAX = 0 1078. C**** 1079. C**** Loop over horizontal quarter boxes in north-south direction 1080. C**** 1081. JQ = 1 1082. 300 J2 = 2*J+JQ-3 1083. RJ = 2.*RRT3*JQ - 3.*RRT3 1084. C**** 1085. C**** Loop over horizontal quarter boxes in east-west direction 1086. C**** 1087. IQ = 1 1088. 400 I2 = 2*I+IQ-2 1089. RI = 2.*RRT3*IQ - 3.*RRT3 1090. C**** Temperature (K) of an eighth box 1091. DO 410 L=1,LMA 1092. TK(2*L,IQ,JQ) = 1093. * (H0M(I,J,L)+RI*HXM(I,J,L)+RJ*HYM(I,J,L)+.5*HZM(I,J,L))* 1094. * PKUP(I,J,L) / (MA(I,J,L)*DXYP(J)*SHCD) 1095. TK(2*L-1,IQ,JQ) = 1096. * (H0M(I,J,L)+RI*HXM(I,J,L)+RJ*HYM(I,J,L)-.5*HZM(I,J,L))* 1097. * PKDN(I,J,L) / (MA(I,J,L)*DXYP(J)*SHCD) 1098. C**** Specific humidity of an eighth box 1099. QQ(2*L,IQ,JQ) = 1100. * (Q0M(I,J,L)+RI*QXM(I,J,L)+RJ*QYM(I,J,L)+.5*QZM(I,J,L)) / 1101. * (MA(I,J,L)*DXYP(J)) 1102. 410 QQ(2*L-1,IQ,JQ) = 1103. * (Q0M(I,J,L)+RI*QXM(I,J,L)+RJ*QYM(I,J,L)-.5*QZM(I,J,L)) / 1104. * (MA(I,J,L)*DXYP(J)) 1105. C**** 1106. C**** Large scale condensation 1107. C**** 1108. 500 CALL LSCOND (MG(1,IQ,JQ), PR(1,IQ,JQ), TK(1,IQ,JQ), QQ(1,IQ,JQ), 1109. * PREC(I,J,2), CLDSS(1,I2,J2)) 1110. IF(QPOLE) GO TO 800 1111. C**** 1112. C**** End of loop over IQ 1113. C**** 1114. IQ = IQ + 1 1115. IF(IQ.le.2) GO TO 400 1116. C**** 1117. C**** End of loop over JQ 1118. C**** 1119. JQ = JQ + 1 1120. IF(JQ.le.2) GO TO 300 1121. C**** 1122. C**** Update model prognostic variables away from poles 1123. C**** 1124. C IF(LSCMAX.le.0) GO TO 750 1125. PREC(I,J,2) = PREC(I,J,2)*.25 1126. DO 710 L=1,LMA ! LSCMAX 1127. L2=2*L 1128. L1=L2-1 1129. c AJL(J,L,26) = AJL(J,L,26) - PLK*(H0M(I,J,L) - 1130. c + .125*SHCD*DXYP(J)*MA(I,J,L)* 1131. c * ((TK(L2,2,2)+TK(L2,1,2)+TK(L2,2,1)+TK(L2,1,1))/PKUP(I,J,L) + 1132. c + (TK(L1,2,2)+TK(L1,1,2)+TK(L1,2,1)+TK(L1,1,1))/PKDN(I,J,L))) 1133. H0M(I,J,L) = .125*SHCD*DXYP(J)*MA(I,J,L)* 1134. * ((TK(L2,2,2)+TK(L2,1,2)+TK(L2,2,1)+TK(L2,1,1))/PKUP(I,J,L) + 1135. + (TK(L1,2,2)+TK(L1,1,2)+TK(L1,2,1)+TK(L1,1,1))/PKDN(I,J,L)) 1136. HXM(I,J,L) = (.125/RRT3)*SHCD*DXYP(J)*MA(I,J,L)* 1137. * ((TK(L2,2,2)-TK(L2,1,2)+TK(L2,2,1)-TK(L2,1,1))/PKUP(I,J,L) + 1138. + (TK(L1,2,2)-TK(L1,1,2)+TK(L1,2,1)-TK(L1,1,1))/PKDN(I,J,L)) 1139. HYM(I,J,L) = (.125/RRT3)*SHCD*DXYP(J)*MA(I,J,L)* 1140. * ((TK(L2,2,2)+TK(L2,1,2)-TK(L2,2,1)-TK(L2,1,1))/PKUP(I,J,L) + 1141. + (TK(L1,2,2)+TK(L1,1,2)-TK(L1,2,1)-TK(L1,1,1))/PKDN(I,J,L)) 1142. HZM(I,J,L) = .25*SHCD*DXYP(J)*MA(I,J,L)* 1143. * ((TK(L2,2,2)+TK(L2,1,2)+TK(L2,2,1)+TK(L2,1,1))/PKUP(I,J,L) - 1144. + (TK(L1,2,2)+TK(L1,1,2)+TK(L1,2,1)+TK(L1,1,1))/PKDN(I,J,L)) 1145. Q0M(I,J,L) = .125*DXYP(J)*MA(I,J,L)* 1146. * ((QQ(L2,2,2)+QQ(L2,1,2)+QQ(L2,2,1)+QQ(L2,1,1)) + 1147. + (QQ(L1,2,2)+QQ(L1,1,2)+QQ(L1,2,1)+QQ(L1,1,1))) 1148. QXM(I,J,L) = (.125/RRT3)*DXYP(J)*MA(I,J,L)* 1149. * ((QQ(L2,2,2)-QQ(L2,1,2)+QQ(L2,2,1)-QQ(L2,1,1)) + 1150. + (QQ(L1,2,2)-QQ(L1,1,2)+QQ(L1,2,1)-QQ(L1,1,1))) 1151. QYM(I,J,L) = (.125/RRT3)*DXYP(J)*MA(I,J,L)* 1152. * ((QQ(L2,2,2)+QQ(L2,1,2)-QQ(L2,2,1)-QQ(L2,1,1)) + 1153. + (QQ(L1,2,2)+QQ(L1,1,2)-QQ(L1,2,1)-QQ(L1,1,1))) 1154. QZM(I,J,L) = .25*DXYP(J)*MA(I,J,L)* 1155. * ((QQ(L2,2,2)+QQ(L2,1,2)+QQ(L2,2,1)+QQ(L2,1,1)) - 1156. + (QQ(L1,2,2)+QQ(L1,1,2)+QQ(L1,2,1)+QQ(L1,1,1))) 1157. IF(Q0M(I,J,L).ge.0.) GO TO 710 1158. WRITE (6,*) 'CONDSE: Negative humidity set to 0. I,J,L,Q0,DQ0=', 1159. * I,J,L,Q0M(I,J,L)/(DXYP(J)*MA(I,J,L)) 1160. Q0M(I,J,L) = 0. 1161. 710 continue 1162. C**** 1163. C**** End of loop over I 1164. C**** 1165. 750 I = I + 1 1166. IF(I.le.IM) GO TO 250 1167. GO TO 850 1168. C**** 1169. C**** Update model prognostic variables at the poles 1170. C**** 1171. C 800 IF(LSCMAX.le.0) GO TO 750 1172. 800 DO 810 L=1,LMA ! LSCMAX 1173. L2=2*L 1174. L1=L2-1 1175. H0M(I,J,L) = .5*SHCD*DXYP(J)*MA(I,J,L)* 1176. * (TK(L2,IQ,JQ)/PKUP(I,J,L) + TK(L1,IQ,JQ)/PKDN(I,J,L)) 1177. HZM(I,J,L) = SHCD*DXYP(J)*MA(I,J,L)* 1178. * (TK(L2,IQ,JQ)/PKUP(I,J,L) - TK(L1,IQ,JQ)/PKDN(I,J,L)) 1179. Q0M(I,J,L) = .5*DXYP(J)*MA(I,J,L)*(QQ(L2,IQ,JQ) + QQ(L1,IQ,JQ)) 1180. QZM(I,J,L) = DXYP(J)*MA(I,J,L)*(QQ(L2,IQ,JQ) - QQ(L1,IQ,JQ)) 1181. IF(Q0M(I,J,L).ge.0.) GO TO 810 1182. WRITE (6,*) 'CONDSE: Negative humidity set to 0. I,J,L,Q0,DQ0=', 1183. * I,J,L,Q0M(I,J,L)/(DXYP(J)*MA(I,J,L)) 1184. Q0M(I,J,L) = 0. 1185. 810 continue 1186. C**** 1187. C**** End of loop over J 1188. C**** 1189. 850 continue 1190. C**** 1191. C**** Convert precipitation temperature into precipitation energy 1192. C**** 1193. 900 DO 930 J=1,JM 1194. IMAX=IM 1195. IF(J.eq.1 .or. J.eq.JM) IMAX=1 1196. DO 930 I=1,IMAX 1197. PREC(I,J,1) = PREC(I,J,1) + PREC(I,J,2) 1198. TPRE = EPRE(I,J,1) 1199. IF(TPRE.ge.0.) then 1200. C EPRE(I,J,1) = PREC(I,J,1)*TPRE*SHCI 1201. C EPRE(I,J,2) = PREC(I,J,2)*TPRE*SHCI 1202. EPRE(I,J,1) = 0. 1203. EPRE(I,J,2) = 0. 1204. else 1205. C EPRE(I,J,1) = PREC(I,J,1)*(TPRE*SHCI-ELHM) 1206. C EPRE(I,J,2) = PREC(I,J,2)*(TPRE*SHCI-ELHM) 1207. EPRE(I,J,1) = -PREC(I,J,1)*ELHM 1208. EPRE(I,J,2) = -PREC(I,J,2)*ELHM 1209. endif 1210. IF(PREC(I,J,1).le.0.) GO TO 920 1211. AIJ(I,J, 5) = AIJ(I,J, 5) + PREC(I,J,1) 1212. AIJ(I,J,23) = AIJ(I,J,23) + EPRE(I,J,1) 1213. AIJ(I,J,52) = AIJ(I,J,52) + EPRE(I,J,1) 1213.5 TMNMX(I,J,5) = TMNMX(I,J,5) + PREC(I,J,1) 1214. PRECZ = PREC(I,J,1)*ZATMO(I,J) 1215. AJ(J,20) = AJ(J,20) + PREC(I,J,1)* FGRND(I,J) 1216. BJ(J,20) = BJ(J,20) + PREC(I,J,1)* FGICE(I,J) 1217. CJ(J,20) = CJ(J,20) + PREC(I,J,1)* FLAKE(I,J)*(1.-RSI(I,J)) 1218. DJ(J,20) = DJ(J,20) + PREC(I,J,1)* FLAKE(I,J)* RSI(I,J) 1219. EJ(J,20) = EJ(J,20) + PREC(I,J,1)*FOCEAN(I,J)*(1.-RSI(I,J)) 1220. FJ(J,20) = FJ(J,20) + PREC(I,J,1)*FOCEAN(I,J)* RSI(I,J) 1221. AJ(J,21) = AJ(J,21) + PRECZ * FGRND(I,J) 1222. BJ(J,21) = BJ(J,21) + PRECZ * FGICE(I,J) 1223. CJ(J,21) = CJ(J,21) + PRECZ * FLAKE(I,J)*(1.-RSI(I,J)) 1224. DJ(J,21) = DJ(J,21) + PRECZ * FLAKE(I,J)* RSI(I,J) 1225. EJ(J,21) = EJ(J,21) + PRECZ *FOCEAN(I,J)*(1.-RSI(I,J)) 1226. FJ(J,21) = FJ(J,21) + PRECZ *FOCEAN(I,J)* RSI(I,J) 1227. AJ(J,39) = AJ(J,39) + EPRE(I,J,1)* FGRND(I,J) 1228. BJ(J,39) = BJ(J,39) + EPRE(I,J,1)* FGICE(I,J) 1229. CJ(J,39) = CJ(J,39) + EPRE(I,J,1)* FLAKE(I,J)*(1.-RSI(I,J)) 1230. DJ(J,39) = DJ(J,39) + EPRE(I,J,1)* FLAKE(I,J)* RSI(I,J) 1231. EJ(J,39) = EJ(J,39) + EPRE(I,J,1)*FOCEAN(I,J)*(1.-RSI(I,J)) 1232. FJ(J,39) = FJ(J,39) + EPRE(I,J,1)*FOCEAN(I,J)* RSI(I,J) 1233. IF(PREC(I,J,2).le.0.) GO TO 920 1234. AJ(J,61) = AJ(J,61) + PREC(I,J,2)* FGRND(I,J) 1235. BJ(J,61) = BJ(J,61) + PREC(I,J,2)* FGICE(I,J) 1236. CJ(J,61) = CJ(J,61) + PREC(I,J,2)* FLAKE(I,J)*(1.-RSI(I,J)) 1237. DJ(J,61) = DJ(J,61) + PREC(I,J,2)* FLAKE(I,J)* RSI(I,J) 1238. EJ(J,61) = EJ(J,61) + PREC(I,J,2)*FOCEAN(I,J)*(1.-RSI(I,J)) 1239. FJ(J,61) = FJ(J,61) + PREC(I,J,2)*FOCEAN(I,J)* RSI(I,J) 1240. C**** 1241. C**** Update age of snow 1242. C**** 1243. 920 EXPSNO = 1. 1244. IF(TPRE.lt.0.) EXPSNO = EXP(-PREC(I,J,1)) 1245. IF(TPRE.gt.0.) TPRE = 0. 1246. DAGE = (1./NDAY - BLDATA(I,J,4)/(NDAY*AGEMAX)) / 1247. * (1.+TPRE*TPRE/100.) 1248. 930 BLDATA(I,J,4) = (BLDATA(I,J,4) + DAGE)*EXPSNO 1249. C**** 1250. C**** Accumulate ADAILY diagnostics 1251. C**** 1252. DO 960 KR=1,4 1253. I=IJDD(1,KR) 1254. J=IJDD(2,KR) 1255. ADAILY(JHOUR, 5,KR) = ADAILY(JHOUR, 5,KR) + EPRE(I,J,1) 1256. 960 ADAILY(JHOUR,49,KR) = ADAILY(JHOUR,49,KR) + PREC(I,J,1) 1257. RETURN 1258. END 1500. 1501. SUBROUTINE LSCOND (MAQ, P, TK, QQ, PREC, CLDSS) 1502. C**** 1503. C**** LSCOND receives one-dimensional arrays of mass and energy and 1504. C**** modifies these arrays based on the large scale condensation 1505. C**** algorithm. 1506. C**** 1507. IMPLICIT REAL*8 (A-H,M,O-Z) 1508. PARAMETER (LMA=9, 1509. * RDRY=287.,RVAP=461.5, SHCD=1003.5, RKAP=RDRY/SHCD, 1510. * ELHE=2500000.,ELHS=2834000., TKF=273.16d0) 1511. REAL*8 MAQ(LMA), P(2*LMA), TK(2*LMA), QQ(2*LMA), CLDSS(LMA) 1512. COMMON /WORK01/ LSCMAX 1513. C**** QSAT = (RAIR/RVAP)*610.71*EXP((L/RVAP)*(1/TKF-1/T)) / PRESS 1514. QSAT(TM,PR) = 379.7915d0*EXP(ELBYRV*(1./TKF-1./TM))/PR 1515. C**** 1516. C**** Large scale condensation 1517. C**** 1518. COND = 0. 1519. EL = ELHS 1520. ELBYCP = EL/SHCD 1521. DO 320 L2=2*LMA,1,-1 1522. L=(L2+1)/2 1523. IF(TK(L2).ge.TKF) EL = ELHE 1524. C IF(TK(L2).le.TIC) EL = ELHS (define TIC = 233.16) 1525. TK(L2) = TK(L2) - (COND/MAQ(L))*ELBYCP 1526. QQ(L2) = QQ(L2) + COND/MAQ(L) 1527. COND = 0. 1528. ELBYCP = EL/SHCD 1529. ELBYRV = EL/RVAP 1530. QS = QSAT(TK(L2),P(L2)) 1531. IF(QQ(L2).le.QS) GO TO 320 1532. DQSDT = QS*ELBYRV/(TK(L2)*TK(L2)) 1533. DQSUM = (QQ(L2)-QS)/(DQSDT*ELBYCP+1.) 1534. TK(L2) = TK(L2) + DQSUM*ELBYCP 1535. QQ(L2) = QQ(L2) - DQSUM 1536. DO 310 N=1,2 1537. QS = QSAT(TK(L2),P(L2)) 1538. DQSDT = QS*ELBYRV/(TK(L2)*TK(L2)) 1539. DQ = (QQ(L2)-QS)/(DQSDT*ELBYCP+1.) 1540. TK(L2) = TK(L2) + DQ*ELBYCP 1541. QQ(L2) = QQ(L2) - DQ 1542. 310 DQSUM = DQSUM + DQ 1543. IF(L2.ne.2*L-1 .or. QQ(L2)+QQ(L2+1).ge.0.) GO TO 315 1544. C**** Reduce condensation so that sum of vapor in two half boxes 1545. C**** is non-negative 1546. DQ = - (QQ(L2)+QQ(L2+1)) 1547. TK(L2) = TK(L2) - DQ*ELBYCP 1548. QQ(L2) = QQ(L2) + DQ 1549. DQSUM = DQSUM - DQ 1550. 315 COND = DQSUM*MAQ(L) 1551. CLDSS(L) = CLDSS(L) + COND 1552. 320 continue 1553. C**** End of inside L loop 1554. PREC = PREC + .5*COND 1555. RETURN 1556. END 2000. 2001. SUBROUTINE ORBIT (OBLIQ,ECCEN,OMEGVP,DAY,DISTQ,SIND,COSD,LAMBDA) 2002. C**** 2003. C**** ORBIT receives the orbital parameters and time of year, and 2004. C**** returns the distance from the sun and its declination angle. 2005. C**** The reference for the following caculations is: V.M.Blanco 2006. C**** and S.W.McCuskey, 1961, "Basic Physics of the Solar System", 2007. C**** pages 135 - 151. 2008. C**** 2011. C**** Input: OBLIQ = latitude of tropics in degrees 2012. C**** ECCEN = eccentricity of the orbital ellipse 2013. C**** OMEGVP = spatial angle from vernal equinox to perihelion 2013.1 C**** in radians with sun as angle vertex 2014. C**** DAY = day of the year in days; 0 = Jan 1, hour 0 2015. C**** 2016. C**** Constants: EDAYPY = Earth days per year = 365 2017. C**** VEREQX = occurence of vernal equinox = day 79 = Mar 21 2018. C**** 2019. C**** Intermediate quantities: 2019.1 C**** BSEMI = semi minor axis in units of the semi major axis 2019.2 C**** TAofVE = TA(VE) = true anomaly of the vernal equinox = -OMEGVP 2019.3 C**** EAofVE = EA(VE) = eccentric anomaly of the vernal equinox 2019.4 C**** MAofVE = MA(VE) = mean anomaly of the vernal equinox 2020. C**** PERIHE = perihelion in days 2021. C**** MA = mean anomaly in temporal radians from perihelion 2022. C**** EA = eccentric anomaly = spatial angle from perihelion 2022.1 C**** to current location with orbit center as angle vertex 2023. C**** TA = true anomaly = spatial angle from perihelion to 2024. C**** current location with sun as angle vertex 2025. C**** GREENW = longitude of Greenwich in the Earth's reference frame 2026. C**** in radians 2027. C**** DIST = distance to the sun in units of the semi major axis 2027.5 C**** 2028. C**** Output: DISTQ = square of the distance of the sun in units of 2028.1 C**** the semi major axis 2029. C**** SIND = sine of the declination angle 2030. C**** COSD = cosine of the declination angle 2031. C**** LAMBDA = sun longitude in Earth's rotating reference frame 2031.1 C*** in radians 2032. C**** 2033. IMPLICIT REAL*8 (A-Z) 2037. PARAMETER (TWOPI=6.283185307179586477d0, EDAYPY=365., VEREQX=79.) 2041. OBLIQR = OBLIQ*(TWOPI/360.) 2043. C**** 2044. C**** Using Kepler's equation: MA = EA - ECCEN sin(EA) 2045. C**** calculate the mean anomaly of the current day 2046. C**** 2047. BSEMI = SQRT(1.-ECCEN*ECCEN) 2048. TAofVE = - OMEGVP 2048.1 EAofVE = ATAN2(SIN(TAofVE)*BSEMI,COS(TAofVE)+ECCEN) 2048.2 MAofVE = EAofVE - ECCEN*SIN(EAofVE) 2049. C PERIHE = VEREQX - MAofVE*EDAYPY/TWOPI 2050. MA = (DAY-VEREQX)*TWOPI/EDAYPY + MAofVE 2051. C**** 2052. C**** Solve for EA using Kepler's equation: MA = EA - ECCEN sin(EA) 2053. C**** 2054. EA = MA + ECCEN*(SIN(MA)+ECCEN*SIN(2.*MA)/2.) 2055. 110 DEA = (MA-EA+ECCEN*SIN(EA)) / (1.-ECCEN*COS(EA)) 2056. EA = EA + DEA 2057. IF(DABS(DEA).gt.1.d-8) GO TO 110 2058. C**** 2059. C**** Calculate the distance to the sun and the true anomaly 2060. C**** 2062. COSEA = COS(EA) 2063. SINEA = SIN(EA) 2064. DISTQ = (1.-ECCEN*COSEA)*(1.-ECCEN*COSEA) 2065. TA = ATAN2(SINEA*BSEMI,COSEA-ECCEN) 2066. C**** 2067. C**** Change the reference frame to be a nonrotating reference frame 2068. C**** with the Earth at the center, the positive x axis parallel to 2069. C**** the ray from the sun to the Earth were it at vernal equinox, and 2069.1 C**** the x-y plane to be the Earth's equatorial plane. 2070. C**** The distance from the current Earth to that ray (or x axis) is: 2071. C**** DIST sin(TA-TAofVE). The sun is located at: 2072. C**** 2073. C**** SUN = (-DIST cos(TA-OMEGVP), 2074. C**** -DIST sin(TA-OMEGVP) cos(OBLIQ), 2075. C**** DIST sin(TA-OMEGVP) sin(OBLIQ)) 2076. C**** 2077. C**** At vernal equinox: mod(TA-TAofVE,2*pi) = 0, SUN = (-DIST,0,0) 2080. C**** 2081. SIND = SIN(TA-TAofVE)*SIN(OBLIQR) 2082. COSD = SQRT(1.-SIND*SIND) 2084. C GREENW = (DAY-VEREQX)*TWOPI*(EDAYPY+1.)/EDAYPY**2 2084.1 C = (MA-MAofVE)*(EDAYPY+1.)/EDAYPY 2085. C SUNX = -COS(TA-TAofVE) 2086. C SUNY = -SIN(TA-TAofVE)*COS(OBLIQR) 2087. C LAMBDA = ATAN2(SUNY,SUNX) - GREENW 2088. C LAMBDA = MOD(LAMBDA,TWOPI) 2089. C**** 2090. RETURN 2091. END 2500. 2501. SUBROUTINE COSZ0 2502. C**** 2503. C**** COSZ0 calculates the Earth's zenith angle, weighted either 2504. C**** by time or by sun light. 2505. C**** 2505.5 IMPLICIT REAL*8 (A-H,M,O-Z) 2506. PARAMETER (IM=72,JM=46, TWOPI=6.283185307179586477d0) 2506.5 REAL*4 VADATA,DGLAT,DGLON 2507. REAL*8 COSZ(IM,JM),COSZA(IM,JM), 2508. * SINJ(JM),COSJ(JM),RI(IM),SINI(IM),COSI(IM),LT1,LT2 2509. COMMON /RADNIO/ VADATA(11,4,3),DGLAT(JM),DGLON(IM) 2510. COMMON /WORK02/ LT1(IM),SLT1(IM),S2LT1(IM), 2511. * LT2(IM),SLT2(IM),S2LT2(IM) 2512. C**** ZERO1 must equal the cut-off value for COSZ used in SOLAR 2514. PARAMETER (ZERO1=1.d-2) 2515. C**** Compute area weighted latitudes and their sines and cosines 2516. DLAT = TWOPI*NINT(360./(JM-1))/720. 2517. FJEQ = (1+JM)/2. 2518. RLATS = -.25*TWOPI 2519. SLATS = -1. 2520. CLATS = 0. 2521. DO 20 J=1,JM 2522. RLATN = DLAT*(J+.5-FJEQ) 2523. IF(J.eq.JM) RLATN = .25*TWOPI 2524. SLATN = SIN(RLATN) 2525. CLATN = COS(RLATN) 2526. RLATM = (RLATN*SLATN+CLATN-RLATS*SLATS-CLATS)/(SLATN-SLATS) 2527. DGLAT(J)= RLATM*(360./TWOPI) 2528. SINJ(J) = SIN(RLATM) 2529. COSJ(J) = COS(RLATM) 2530. RLATS = RLATN 2531. SLATS = SLATN 2532. 20 CLATS = CLATN 2533. C**** Compute sines and cosines of longitude. The left edge of 2534. C**** grid box I=1 is assumed to on the International Date Line. 2535. DO 40 I=1,IM 2536. RI(I) = (TWOPI/IM)*(I-.5) 2537. DGLON(I)= (I-.5)*(360./IM) 2538. SINI(I) = SIN(RI(I)) 2539. 40 COSI(I) = COS(RI(I)) 2540. RETURN 2541. C**** 2542. C**** 2543. ENTRY COSZT (SIND,COSD,ROT1,ROT2,COSZ) 2544. C**** 2545. C**** COSZT computes the zenith angle weighted by time from ROT1 2546. C**** to ROT2, Greenwich Mean Time in radians. 2547. C**** Input: SIND,COSD = sine and cosine of the declination angle 2548. C**** ROT1 = initial time, 0 < ROT1 < 2*PI 2549. C**** ROT2 = final time, ROT1 < ROT2 < ROT1+2*PI 2550. C**** Output: COSZ = S(cos(Z)*dT) / S(dT) 2551. C**** 2552. DROT = ROT2-ROT1 2553. C**** Compute sines and cosines of initial and final times 2554. 100 SR1 = SIN(ROT1) 2555. CR1 = COS(ROT1) 2556. SR2 = SIN(ROT2) 2557. CR2 = COS(ROT2) 2558. C**** Compute initial and final local times (measured from noon) 2559. C**** and their sines and cosines 2560. DO 120 I=1,IM 2561. LT1(I) = ROT1+RI(I) 2562. SLT1(I) = SR1*COSI(I)+CR1*SINI(I) 2563. LT2(I) = ROT2+RI(I) 2564. 120 SLT2(I) = SR2*COSI(I)+CR2*SINI(I) 2565. C**** 2566. C**** CALCULATION FOR POLAR GRID BOXES 2567. C**** 2568. DO 200 J=1,JM,JM-1 2569. SJSD = SINJ(J)*SIND 2570. CJCD = COSJ(J)*COSD 2571. IF(SJSD+CJCD.LE.ZERO1) GO TO 180 2572. IF(SJSD-CJCD.ge.0.) GO TO 160 2573. C**** AVERAGE COSZ FROM DAWN TO DUSK NEAR THE POLES 2574. DUSK = ACOS(-SJSD/CJCD) 2575. SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD 2576. DAWN = -DUSK 2577. SDAWN = -SDUSK 2578. COSZ(1,J) = (SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN))/TWOPI 2579. GO TO 200 2580. C**** CONSTANT DAYLIGHT NEAR THE POLES 2581. 160 COSZ(1,J) = SJSD 2582. GO TO 200 2583. C**** CONSTANT NIGHTIME NEAR THE POLES 2584. 180 COSZ(1,J) = 0. 2585. 200 CONTINUE 2586. C**** 2587. C**** LOOP OVER NON-POLAR LATITUDES 2588. C**** 2589. DO 500 J=2,JM-1 2590. SJSD = SINJ(J)*SIND 2591. CJCD = COSJ(J)*COSD 2592. IF(SJSD+CJCD.LE.ZERO1) GO TO 460 2593. IF(SJSD-CJCD.ge.0.) GO TO 420 2594. C**** COMPUTE DAWN AND DUSK (AT LOCAL TIME) AND THEIR SINES 2595. DUSK = ACOS(-SJSD/CJCD) 2596. SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD 2597. DAWN = -DUSK 2598. SDAWN = -SDUSK 2599. C**** NEITHER CONSTANT DAYTIME NOR CONSTANT NIGHTIME AT THIS LATITUDE, 2600. C**** LOOP OVER LONGITUDES 2601. ZERO2 = ZERO1/CJCD 2602. DO 400 I=1,IM 2603. C**** FORCE DUSK TO LIE BETWEEN LT1 AND LT1+2*PI 2604. IF(DUSK.gt.LT1(I)+ZERO2) GO TO 220 2605. DUSK = DUSK+TWOPI 2606. DAWN = DAWN+TWOPI 2607. 220 IF(DAWN.LT.LT2(I)-ZERO2) GO TO 240 2608. C**** CONTINUOUS NIGHTIME FROM INITIAL TO FINAL TIME 2609. COSZ(I,J) = 0. 2610. GO TO 400 2611. 240 IF(DAWN.ge.LT1(I)) GO TO 300 2612. IF(DUSK.LT.LT2(I)) GO TO 260 2613. C**** CONTINUOUS DAYLIGHT FROM INITIAL TIME TO FINAL TIME 2614. COSZ(I,J) = SJSD+CJCD*(SLT2(I)-SLT1(I))/DROT 2615. GO TO 400 2616. 260 IF(DAWN+TWOPI.LT.LT2(I)-ZERO2) GO TO 280 2617. C**** DAYLIGHT AT INITIAL TIME AND NIGHT AT FINAL TIME 2618. COSZ(I,J) = (SJSD*(DUSK-LT1(I))+CJCD*(SDUSK-SLT1(I)))/DROT 2619. GO TO 400 2620. C**** DAYLIGHT AT INITIAL AND FINAL TIMES WITH NIGHTIME IN BETWEEN 2621. 280 COSZ(I,J) = (SJSD*(LT2(I)-DAWN-TWOPI+DUSK-LT1(I))+CJCD* 2622. * (SLT2(I)-SDAWN+SDUSK-SLT1(I)))/DROT 2623. GO TO 400 2624. 300 IF(DUSK.LT.LT2(I)) GO TO 320 2625. C**** NIGHT AT INITIAL TIME AND DAYLIGHT AT FINAL TIME 2626. COSZ(I,J) = (SJSD*(LT2(I)-DAWN)+CJCD*(SLT2(I)-SDAWN))/DROT 2627. GO TO 400 2628. C**** NIGHTIME AT INITIAL AND FINAL TIMES WITH DAYLIGHT IN BETWEEN 2629. 320 COSZ(I,J) = (SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN))/DROT 2630. 400 CONTINUE 2631. GO TO 500 2632. C**** CONSTANT DAYLIGHT AT THIS LATITUDE 2633. 420 DO 440 I=1,IM 2634. 440 COSZ(I,J) = SJSD+CJCD*(SLT2(I)-SLT1(I))/DROT 2635. GO TO 500 2636. C**** CONSTANT NIGHTIME AT THIS LATITUDE 2637. 460 DO 480 I=1,IM 2638. 480 COSZ(I,J) = 0. 2639. 500 CONTINUE 2640. RETURN 2641. C**** 2642. C**** 2643. ENTRY COSZS (SIND,COSD,ROT1,ROT2,COSZ,COSZA) 2644. C**** 2645. C**** COSZS computes the zenith angle from ROT1 to ROT2 twice, 2646. C**** weighted by time and weighted by sun light. 2647. C**** Input: SIND,COSD = sine and cosine of the declination angle 2648. C**** ROT1 = initial time, 0 < ROT1 < 2*PI 2649. C**** ROT2 = final time, ROT1 < ROT2 < ROT1+2*PI 2650. C**** Output: COSZ = S(cos(Z)*dT) / S(dT) 2651. C**** COSZA = S(cos(Z)*cos(Z)*dT) / S(cos(Z)*dT) 2652. C**** 2653. DROT = ROT2-ROT1 2654. C**** COMPUTE THE SINES AND COSINES OF THE INITIAL AND FINAL GMT'S 2655. SR1 = SIN(ROT1) 2656. CR1 = COS(ROT1) 2657. SR2 = SIN(ROT2) 2658. CR2 = COS(ROT2) 2659. C**** COMPUTE THE INITIAL AND FINAL LOCAL TIMES (MEASURED FROM NOON TO 2660. C**** NOON) AND THEIR SINES AND COSINES 2661. DO 520 I=1,IM 2662. LT1(I) = ROT1+RI(I) 2663. SLT1(I) = SR1*COSI(I)+CR1*SINI(I) 2664. CLT1 = CR1*COSI(I)-SR1*SINI(I) 2665. S2LT1(I)= 2.*SLT1(I)*CLT1 2666. LT2(I) = ROT2+RI(I) 2667. SLT2(I) = SR2*COSI(I)+CR2*SINI(I) 2668. CLT2 = CR2*COSI(I)-SR2*SINI(I) 2669. 520 S2LT2(I)= 2.*SLT2(I)*CLT2 2670. C**** 2671. C**** CALCULATION FOR POLAR GRID BOXES 2672. C**** 2673. DO 600 J=1,JM,JM-1 2674. SJSD = SINJ(J)*SIND 2675. CJCD = COSJ(J)*COSD 2676. IF(SJSD+CJCD.LE.ZERO1) GO TO 580 2677. IF(SJSD-CJCD.ge.0.) GO TO 560 2678. C**** AVERAGE COSZ FROM DAWN TO DUSK NEAR THE POLES 2679. CDUSK = -SJSD/CJCD 2680. DUSK = ACOS(CDUSK) 2681. SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD 2682. S2DUSK= 2.*SDUSK*CDUSK 2683. DAWN = -DUSK 2684. SDAWN = -SDUSK 2685. S2DAWN= -S2DUSK 2686. ECOSZ = SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN) 2687. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SDAWN)+ 2688. * .5*CJCD*(DUSK-DAWN+.5*(S2DUSK-S2DAWN))) 2689. COSZ(1,J) = ECOSZ/TWOPI 2690. COSZA(1,J) = QCOSZ/ECOSZ 2691. GO TO 600 2692. C**** CONSTANT DAYLIGHT NEAR THE POLES 2693. 560 ECOSZ = SJSD*TWOPI 2694. QCOSZ = SJSD*ECOSZ+.5*CJCD*CJCD*TWOPI 2695. COSZ(1,J) = ECOSZ/TWOPI 2696. COSZA(1,J) = QCOSZ/ECOSZ 2697. GO TO 600 2698. C**** CONSTANT NIGHTIME NEAR THE POLES 2699. 580 COSZ(1,J) = 0. 2700. COSZA(1,J) = 0. 2701. 600 CONTINUE 2702. C**** 2703. C**** LOOP OVER NON-POLAR LATITUDES 2704. C**** 2705. DO 900 J=2,JM-1 2706. SJSD = SINJ(J)*SIND 2707. CJCD = COSJ(J)*COSD 2708. IF(SJSD+CJCD.LE.ZERO1) GO TO 860 2709. IF(SJSD-CJCD.ge.0.) GO TO 820 2710. C**** COMPUTE DAWN AND DUSK (AT LOCAL TIME) AND THEIR SINES 2711. CDUSK = -SJSD/CJCD 2712. DUSK = ACOS(CDUSK) 2713. SDUSK = SQRT(CJCD*CJCD-SJSD*SJSD)/CJCD 2714. S2DUSK= 2.*SDUSK*CDUSK 2715. DAWN = -DUSK 2716. SDAWN = -SDUSK 2717. S2DAWN= -S2DUSK 2718. C**** NEITHER CONSTANT DAYTIME NOR CONSTANT NIGHTIME AT THIS LATITUDE, 2719. C**** LOOP OVER LONGITUDES 2720. ZERO2 = ZERO1/CJCD 2721. DO 800 I=1,IM 2722. C**** FORCE DUSK TO LIE BETWEEN LT1 AND LT1+2*PI 2723. IF(DUSK.gt.LT1(I)+ZERO2) GO TO 620 2724. DUSK = DUSK+TWOPI 2725. DAWN = DAWN+TWOPI 2726. 620 IF(DAWN.LT.LT2(I)-ZERO2) GO TO 640 2727. C**** CONTINUOUS NIGHTIME FROM INITIAL TO FINAL TIME 2728. COSZ(I,J) = 0. 2729. COSZA(I,J) = 0. 2730. GO TO 800 2731. 640 IF(DAWN.ge.LT1(I)) GO TO 700 2732. IF(DUSK.LT.LT2(I)) GO TO 660 2733. C**** CONTINUOUS DAYLIGHT FROM INITIAL TIME TO FINAL TIME 2734. ECOSZ = SJSD*DROT+CJCD*(SLT2(I)-SLT1(I)) 2735. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SLT2(I)-SLT1(I))+ 2736. * .5*CJCD*(DROT+.5*(S2LT2(I)-S2LT1(I)))) 2737. COSZ(I,J) = ECOSZ/DROT 2738. COSZA(I,J) = QCOSZ/ECOSZ 2739. GO TO 800 2740. 660 IF(DAWN+TWOPI.LT.LT2(I)-ZERO2) GO TO 680 2741. C**** DAYLIGHT AT INITIAL TIME AND NIGHT AT FINAL TIME 2742. ECOSZ = SJSD*(DUSK-LT1(I))+CJCD*(SDUSK-SLT1(I)) 2743. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SLT1(I))+ 2744. * .5*CJCD*(DUSK-LT1(I)+.5*(S2DUSK-S2LT1(I)))) 2745. COSZ(I,J) = ECOSZ/DROT 2746. COSZA(I,J) = QCOSZ/ECOSZ 2747. GO TO 800 2748. C**** DAYLIGHT AT INITIAL AND FINAL TIMES WITH NIGHTIME IN BETWEEN 2749. 680 ECOSZ = SJSD*(DROT-DAWN-TWOPI+DUSK)+ 2750. * CJCD*(SLT2(I)-SDAWN+SDUSK-SLT1(I)) 2751. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SLT1(I)+SLT2(I)-SDAWN)+ 2752. * .5*CJCD*(DUSK+DROT-DAWN-TWOPI+ 2753. * .5*(S2DUSK-S2LT1(I)+S2LT2(I)-S2DAWN))) 2754. COSZ(I,J) = ECOSZ/DROT 2755. COSZA(I,J) = QCOSZ/ECOSZ 2756. GO TO 800 2757. 700 IF(DUSK.LT.LT2(I)) GO TO 720 2758. C**** NIGHT AT INITIAL TIME AND DAYLIGHT AT FINAL TIME 2759. ECOSZ = SJSD*(LT2(I)-DAWN)+CJCD*(SLT2(I)-SDAWN) 2760. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SLT2(I)-SDAWN)+ 2761. * .5*CJCD*(LT2(I)-DAWN+.5*(S2LT2(I)-S2DAWN))) 2762. COSZ(I,J) = ECOSZ/DROT 2763. COSZA(I,J) = QCOSZ/ECOSZ 2764. GO TO 800 2765. C**** NIGHTIME AT INITIAL AND FINAL TIMES WITH DAYLIGHT IN BETWEEN 2766. 720 ECOSZ = SJSD*(DUSK-DAWN)+CJCD*(SDUSK-SDAWN) 2767. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SDUSK-SDAWN)+ 2768. * .5*CJCD*(DUSK-DAWN+.5*(S2DUSK-S2DAWN))) 2769. COSZ(I,J) = ECOSZ/DROT 2770. COSZA(I,J) = QCOSZ/ECOSZ 2771. 800 CONTINUE 2772. GO TO 900 2773. C**** CONSTANT DAYLIGHT AT THIS LATITUDE 2774. 820 DO 840 I=1,IM 2775. ECOSZ = SJSD*DROT+CJCD*(SLT2(I)-SLT1(I)) 2776. QCOSZ = SJSD*ECOSZ+CJCD*(SJSD*(SLT2(I)-SLT1(I))+ 2777. * .5*CJCD*(DROT+.5*(S2LT2(I)-S2LT1(I)))) 2778. COSZ(I,J) = ECOSZ/DROT 2779. 840 COSZA(I,J) = QCOSZ/ECOSZ 2780. GO TO 900 2781. C**** CONSTANT NIGHTIME AT THIS LATITUDE 2782. 860 DO 880 I=1,IM 2783. COSZ(I,J) = 0. 2784. 880 COSZA(I,J) = 0. 2785. 900 CONTINUE 2786. RETURN 2787. END 3000. 3001. SUBROUTINE RADIA 3002. C**** 3003. C**** RADIA calculates the solar and thermal radiation heating 3004. C**** rates and applies that heating to the atmospheric energy. 3005. C**** 3006. INCLUDE 'C070.COM' 3007. REAL*8 COE(3),PT(6) 3008. INTEGER*4 MODRT,IQUART(0:3),JQUART(0:3) 3009. LOGICAL*4 QPOLE 3010. PARAMETER (OBLIQ=23.44d0, ECCEN=.0167d0, OMEGVP=4.938d0, 3011. * TAUMCW=.075d0, TAUMCI=.025d0, TAUSSW=.3d0, TAUSSI=.01d0, 3012. * TAUMIN=0., TAUCLD=1., TKF=273.16d0, TCIR=243.16d0, 3013. * RHOW=1000., S1BYG1=.57735d0, STBO=5.67032d-8) 3014. COMMON /FLUXCB/ CLDMC(LMA,2*IM,2*JM-2),CLDSS(LMA,2*IM,2*JM-2) 3015. COMMON /WORK00/ PKDN(IM,JM,LMA),PKUP(IM,JM,LMA) 3016. COMMON /WORK01/ COSZ1(IM,JM),COSZ2(IM,JM),COSZA(IM,JM),TOTCLD(LMA) 3017. C COMMON /WORK04/ is being used by the radiation routines 3017.5 REAL*8 AF(JM,KAJ,6) 3017.6 EQUIVALENCE (AF,AJ) 3018. C**** 3019. C**** RA85CM9.COM Radiation Common Block 3/29/91 3020. C**** 3021. PARAMETER (NKSR=6,NKTR=25) 3022. REAL*4 VADATA,DGLAT,DGLON, FGOLDH,FRACSL,TKCICE, 3023. * PLE,HLE,TLB,TLT,TL,QL, U0GAS,ULGAS,TRACER,RTAU, 3024. * POWATR,PGRND,PWICE,PGICE,AGESN,SNOWGR,SNOWWI,SNOWGI, 3025. * TGOW,TGGR,TGWI,TGGI,TS,WS,WEARTH, S0,COSZ,PVT,BXA,SRBXAL, PPMV, 3026. * TRDFLB,TRUFLB,TRNFLB,TRFCRL,TRSLCR, 3027. * SRDFLB,SRUFLB,SRNFLB,SRFHRL,SRSLHR, 3028. * SRIVIS,SROVIS,PLAVIS,SRINIR,SRONIR,PLANIR,SRXATM, 3029. * SRDVIS,SRUVIS,ALBVIS,SRDNIR,SRUNIR,ALBNIR,FSRNFG, 3030. * SRTVIS,SRRVIS,SRAVIS,SRTNIR,SRRNIR,SRANIR,FTRUFG, 3031. * TRDFGW,TRUFGW,TRUFTW,BTEMPW,TRDFSL,TRUFSL,DTRUFG, 3032. * TRSLTS,TRSLTG,TRSLWV,TRSLBS,TTRUFG 3033. C**** Control and input parameters 3034. COMMON /RADNIO/ VADATA(11,4,3),DGLAT(JM),DGLON(IM), 3035. A FGOLDH(18),ID5(5),ITR(4),FRACSL,TKCICE,KFORCE, 3036. C**** Basic RADNIO input data 3037. B PLE(40),HLE(40),TLB(40),TLT(40),TL(40),QL(40), 3038. C U0GAS(40,9),ULGAS(40,9),TRACER(40,4),RTAU(40), 3039. D POWATR,PGRND,PWICE,PGICE,AGESN,SNOWGR,SNOWWI,SNOWGI, 3040. E TGOW,TGGR,TGWI,TGGI,TS,WS,WEARTH, 3041. F S0,COSZ,PVT(11),BXA(153),SRBXAL(15,2), 3042. G PPMV(9),JYEARR,JDAYR,JLAT,ILON, 3043. C**** Basic RADNIO ouput data 3044. H TRDFLB(40),TRUFLB(40),TRNFLB(40),TRFCRL(40),TRSLCR, 3045. I SRDFLB(40),SRUFLB(40),SRNFLB(40),SRFHRL(40),SRSLHR, 3046. J SRIVIS,SROVIS,PLAVIS,SRINIR,SRONIR,PLANIR,SRXATM(4), 3047. K SRDVIS,SRUVIS,ALBVIS,SRDNIR,SRUNIR,ALBNIR,FSRNFG(4), 3048. L SRTVIS,SRRVIS,SRAVIS,SRTNIR,SRRNIR,SRANIR,FTRUFG(4), 3049. M TRDFGW,TRUFGW,TRUFTW,BTEMPW,TRDFSL,TRUFSL,DTRUFG(4), 3050. N TRSLTS,TRSLTG,TRSLWV,TRSLBS,TTRUFG,LBOTCL,LTOPCL 3051. C**** 3052. DATA IFIRST/1/, JYLAST/0/, JDLAST/0/, 3052.1 * IQUART/-1,0,0,-1/, JQUART/-2,-1,-2,-1/ 3053. C**** 3054. C**** FIXDCB FGRND Ground fraction (1) 3055. C**** FGICE Glacial ice fraction (1) 3056. C**** FWATER Water fraction (1) 3057. C**** 3057.1 C**** OCENCB MO Ocean mass (kg/m^2) 3057.2 C**** G0M Mean potential enthalpy of ocean (J) 3057.3 C**** GZM Vertical gradient of ocean potential enthalpy (J) 3057.4 C**** S0M Mean salt of ocean (J) 3057.5 C**** SZM Vertical gradient of ocean salt (kg) 3057.6 C**** GZM(1) Lake surface temperature ( C) 3057.7 C**** 3058. C**** GRNDCB SBG Compressed snow depth of bare ground (m) 3059. C**** SVG Compressed snow depth of vegetated ground (m) 3064. C**** FVEG Ground ratios for 8 vegetation types (1) 3065. C**** 3066. C**** GICECB MGI Mass of glacial ice (kg/m^2) 3067. C**** HGI Enthalpy minus latent heat of glacial ice (J/m^2) 3068. C**** 3069. C**** SEAICB RSI Horizontal ratio of sea ice to ocean (1) 3070. C**** MSI Mass of sea ice (kg/m^2) 3071. C**** HSI Enthalpy minus latent heat of sea ice (J/m^2) 3072. C**** 3073. C**** SURFCB BLDATA(1) Ground temperature over ground fraction (C) 3073.5 C**** BLDATA(2) Ground wetness over ground fraction 3074. C**** BLDATA(3) Composite surface wind speed (m/s) 3075. C**** BLDATA(4) Age of snow (days) 3076. C**** 3078. C**** FLUXCB CLDMC Moist convective cloud precipitation (kg/m^2) 3079. C**** CLDSS Super saturation cloud precipitation (kg/m^2) 3080. C**** 3081. IF(QRAD) IDACC(2)=IDACC(2)+1 3082. IF(IFIRST.ne.1) GO TO 40 3083. IFIRST=0 3084. CALL COSZ0 3085. C**** Set unchanging parameters used in this routine 3086. LMR=LMA+3 3087. DO 10 L=1,LMA 3088. 10 PLE(L) = (SIGEA(L-1)*(MDRYA-MSTRAT)+MSTRAT)*(GRAV/100.) 3089. PLE(LMA+1) = MSTRAT*( GRAV/100.) 3090. PLE(LMA+2) = MSTRAT*(.5*GRAV/100.) 3091. PLE(LMR) = MSTRAT*(.2d0*GRAV/100.) 3092. PLE(LMR+1) = 1.d-5 3093. DO 20 LR=LMA+1,LMR 3094. COE(LR-LMA) = NRAD*DTS*GRAV/(100.*(PLE(LR)-PLE(LR+1))) 3095. QL(LR) = .3d-5 3096. 20 RTAU(LR)= 0. 3096.5 CALL RCOMP0 3097. C**** Set cloud layers used by diagnostics 3098. DO 31 LLOW=1,LTM 3099. IF(SIGA(LLOW+1)*(MDRYA-MSTRAT)+MSTRAT.lt.8012.) GO TO 32 3100. 31 CONTINUE 3101. 32 DO 33 LMID=LLOW+1,LTM 3102. IF(SIGA(LMID+1)*(MDRYA-MSTRAT)+MSTRAT.lt.4383.) GO TO 34 3103. 33 CONTINUE 3104. 34 LHI=LTM 3105. WRITE (6,903) LLOW,LLOW+1,LMID,LMID+1,LHI 3106. C**** Change greenhouse gases and sulphate aerosols each year 3107. 40 IF(JYEAR.eq.JYLAST) GO TO 50 3108. JYLAST=JYEAR 3109. JYEARR=JYEAR 3110. CALL RCOMPY 3115. C**** Compute declination angle, distance to the sun, earth 3116. C**** albedoes, and other parameters at beginning of each new day 3117. 50 IF(JDAY.eq.JDLAST) GO TO 60 3118. JDLAST=JDAY 3119. JDAYR=JDAY 3120. CALL ORBIT (OBLIQ,ECCEN,OMEGVP,JDAY-.5d0,DISTQ,SIND,COSD,SLONG) 3121. S0 = 1367./DISTQ 3123. CALL RCOMPD 3124. C**** Calculate average cosine of zenith angle for current source 3125. C**** term time step and for radiation period 3126. 60 ROT1 = MOD(IDT,NDAY)*TWOPI/NDAY 3127. ROT2 = ROT1 + TWOPI/NDAY 3128. CALL COSZT (SIND,COSD,ROT1,ROT2,COSZ1) 3129. IF(.NOT.QRAD) GO TO 800 3130. ROT2 = ROT1 + TWOPI*NRAD/NDAY 3131. CALL COSZS (SIND,COSD,ROT1,ROT2,COSZ2,COSZA) 3131.5 MODRT = MOD(IDT-IDT0,NRAD*4)/NRAD 3132. C**** 3133. C**** MAIN J LOOP 3134. C**** 3135. DO 700 J=1,JM 3136. JLAT=J 3137. XFRADJ = .2d0 + 1.2d0*COSP(J)*COSP(J) 3138. CALL RCOMPJ 3139. IMAX=IM 3140. QPOLE = J.eq.1 .or. J.eq.JM 3141. IF(.NOT.QPOLE) GO TO 100 3142. C**** Conditions at the poles 3143. QPOLE = .TRUE. 3144. IMAX=1 3145. IQ =1 3146. JQ =1 3147. IF(J.eq.JM) JQ=JM*2-2 3148. C**** 3149. C**** MAIN I LOOP 3150. C**** 3151. 100 DO 600 I=1,IMAX 3152. ILON=I 3153. IF(QPOLE) GO TO 110 3154. IQ=I*2+IQUART(MODRT) 3155. JQ=J*2+JQUART(MODRT) 3157. C**** Determine fractions for surface types 3158. 110 POWATR= FWATER(I,J)*(1.-RSI(I,J)) 3159. PWICE = FWATER(I,J)*RSI(I,J) 3160. PGICE = FGICE(I,J) 3161. PGRND = FGRND(I,J) 3162. PT(1) = FGRND(I,J) 3162.1 PT(2) = FGICE(I,J) 3162.2 PT(3) = FLAKE(I,J)*(1.-RSI(I,J)) 3162.3 PT(4) = FLAKE(I,J)* RSI(I,J) 3162.4 PT(5) = FOCEAN(I,J)*(1.-RSI(I,J)) 3162.5 PT(6) = FOCEAN(I,J)* RSI(I,J) 3163. C**** 3164. C**** Determine cloud optical depths, accumulate cloud diagnostics 3165. C**** 3166. CCOVER = 0. 3167. TAUSUM = 0. 3168. DO 250 L=LMA,1,-1 3169. C**** Load vertical arrays for pressure (mb), temperature (K), and 3170. C**** water vapor mixing ration (1) 3171. PLE(L) = PLE(L+1) + MA(I,J,L)*GRAV/100. 3172. TL(L) = ((H0M(I,J,L)-.5*HZM(I,J,L))*PKDN(I,J,L) + 3172.1 * (H0M(I,J,L)+.5*HZM(I,J,L))*PKUP(I,J,L))*.5 / 3172.2 C * (MA(I,J,L)*DXYP(J)*SHCD + Q0M(I,J,L)*(SHCV-SHCD)) 3172.3 * (MA(I,J,L)*DXYP(J)*SHCD) 3173. QL(L) = Q0M(I,J,L) / (MA(I,J,L)*DXYP(J)) 3174. RTAU(L)= 0. 3175. C**** Calculate super saturation and/or moist convective optical depths 3176. IF(CLDMC(L,IQ,JQ)+CLDSS(L,IQ,JQ).le.0.) GO TO 250 3177. IF(TL(L).lt.TKF) GO TO 210 3178. RTAU(L) = SQRT(MA(I,J,L)*(CLDMC(L,IQ,JQ)*TAUMCW*TAUMCW + 3179. * CLDSS(L,IQ,JQ)*TAUSSW*TAUSSW)) 3180. GO TO 230 3181. 210 IF(TL(L).le.TCIR) GO TO 220 3182. RTAU(L) = SQRT(MA(I,J,L)*(CLDMC(L,IQ,JQ)* 3183. * ((TAUMCW*(TL(L)-TCIR)+TAUMCI*(TKF-TL(L)))/(TKF-TCIR))**2 + 3184. * CLDSS(L,IQ,JQ)* 3185. * ((TAUSSW*(TL(L)-TCIR)+TAUSSI*(TKF-TL(L)))/(TKF-TCIR))**2)) 3186. GO TO 230 3187. 220 RTAU(L) = SQRT(MA(I,J,L)*(CLDMC(L,IQ,JQ)*TAUMCI*TAUMCI + 3188. * CLDSS(L,IQ,JQ)*TAUSSI*TAUSSI)) 3189. C**** Accumulate cloud diagnostics 3190. C**** Determine cloud top pressure (mb) and cloud top temperature (K) 3191. 230 IF(TAUSUM.ge.TAUCLD .or. TAUSUM+RTAU(L).lt.TAUCLD) GO TO 240 3192. CTPRES = PLE(L+1) + (PLE(L)-PLE(L+1))*(TAUCLD-TAUSUM)/RTAU(L) 3193. CTTEMP = TL(L) 3194. CCOVER = 1. 3195. 240 TAUSUM = TAUSUM + RTAU(L) 3196. IF(CLDMC(L,IQ,JQ).ge.CLDSS(L,IQ,JQ)) then 3197. AJL(J,L,15) = AJL(J,L,15) + RTAU(L) 3198. else 3199. AJL(J,L,16) = AJL(J,L,16) + RTAU(L) 3200. endif 3201. 250 CONTINUE 3202. IF(TAUSUM.le.TAUMIN) GO TO 300 3203. DO 251 IT=1,6 3204. 251 AF(J,68,IT) = AF(J,68,IT) + PT(IT)*TAUSUM 3207. AIJ(I,J,18) = AIJ(I,J,18) + TAUSUM 3208. IF(CCOVER.eq.0.) GO TO 280 3209. DO 252 IT=1,6 3210. AF(J,57,IT) = AF(J,57,IT) + PT(IT)*CTPRES 3211. AF(J,59,IT) = AF(J,59,IT) + PT(IT)*CTTEMP 3213. 252 AF(J,69,IT) = AF(J,69,IT) + PT(IT) 3221. AIJ(I,J,19) = AIJ(I,J,19) + 1. 3222. AIJ(I,J,39) = AIJ(I,J,39) + CTPRES 3223. AIJ(I,J,40) = AIJ(I,J,40) + CTTEMP 3224. IF(RTAU(6)+RTAU(7)+RTAU(8)+RTAU(9) .lt. TAUCLD) GO TO 260 3225. AIJ(I,J,43) = AIJ(I,J,43) + 1. 3226. GO TO 280 3227. 260 IF(TAUSUM-RTAU(1)-RTAU(2)-RTAU(3) .lt. TAUCLD) GO TO 270 3228. AIJ(I,J,42) = AIJ(I,J,42) + 1. 3229. GO TO 280 3230. 270 AIJ(I,J,41) = AIJ(I,J,41) + 1. 3231. 280 DO 281 KR=1,4 3232. 281 IF(I.eq.IJDD(1,KR) .and. J.eq.IJDD(2,KR)) GO TO 282 3233. GO TO 400 3234. 282 JH=JHOUR 3235. DO 283 NR=1,NRAD 3236. IF(JH.ge.24) JH=JH-24 3237. ADAILY(JH,21,KR) = ADAILY(JH,21,KR) + RTAU(6) 3238. ADAILY(JH,22,KR) = ADAILY(JH,22,KR) + RTAU(5) 3239. ADAILY(JH,23,KR) = ADAILY(JH,23,KR) + RTAU(4) 3240. ADAILY(JH,24,KR) = ADAILY(JH,24,KR) + RTAU(3) 3241. ADAILY(JH,25,KR) = ADAILY(JH,25,KR) + RTAU(2) 3242. ADAILY(JH,26,KR) = ADAILY(JH,26,KR) + RTAU(1) 3243. ADAILY(JH,27,KR) = ADAILY(JH,27,KR) + CCOVER 3244. 283 JH=JH+1 3245. GO TO 400 3246. C**** 3247. C**** Zero out all optical depths if TAUSUM < TAUMIN 3248. C**** 3249. 300 IF(TAUSUM.le.0.) GO TO 400 3250. DO 310 L=1,LMA 3251. 310 RTAU(L) = 0. 3252. C**** 3253. C**** RADIATION, SOLAR AND THERMAL 3254. C**** 3254.5 400 COSZ = COSZA(I,J) 3254.7 WS = BLDATA(I,J,3) 3254.8 AGESN = BLDATA(I,J,4) 3255. DO 420 K=1,3 3256. 420 TL(LMA+K) = RQT(I,J,K)/SHCD 3261. IF(FOCEAN(I,J).gt.0.) then 3262. GO1 = (G0M(I,J,1)-GZM(I,J,1))/(MO(I,J,1)*DXYP(J)) 3263. SO1 = (S0M(I,J,1)-SZM(I,J,1))/(MO(I,J,1)*DXYP(J)) 3264. TGOW = TEMGS(GO1,SO1) + TKF 3265. endif 3265.5 IF(FLAKE(I,J).gt.0.) then 3266. TGOW = GZM(I,J,1) + TKF 3267. endif 3268. IF(PWICE.gt.0.) then 3268.1 TGWI = (HSI(I,J,1)/(XSI1*MSI(I,J,1))+ELHM) / SHCI + TKF 3268.2 SNOWWI = MSI(I,J,1) - ACE1I 3268.3 endif 3269. IF(PGICE.gt.0.) then 3269.1 TGGI = (HGI(I,J,1)/(XGI1*MGI(I,J))+ELHM) / SHCI + TKF 3269.2 SNOWGI = MGI(I,J) - ACE1I 3269.3 endif 3269.5 IF(PGRND.gt.0.) then 3270. TGGR = BLDATA(I,J,1) + TKF 3270.1 WEARTH= BLDATA(I,J,2) 3270.2 FBARE = FVEG(I,J,1) + FVEG(I,J,10) 3270.3 SNOWGR= (SBG(I,J)*FBARE + SVG(I,J)*(1.-FBARE)) * RHOW 3270.4 DO 430 K=1,10 3270.5 430 PVT(K) = FVEG(I,J,K) 3275.1 endif 3275.9 PGK= (PLE(1)*100.)**RKAP 3276. TS = (H0M(I,J,1)-HZM(I,J,1)*S1BYG1)*PGK / (SHCD*DXYP(J)*MA(I,J,1)) 3277. FGOLDH(2) = XFRADJ*(1.-PGRND) 3278. FGOLDH(3) = XFRADJ*PGRND 3278.1 FGOLDH(5) = .005d0*(WS+.75) * PT(5) 3279. CALL RCOMP 3281. C**** Store solar and thermal heating rates for subsequent 4 hours 3282. SRHR(I,J,0) = SRNFLB(1) 3283. TRHR(I,J,0) = STBO*(POWATR*TGOW**4 + PWICE*TGWI**4 + 3284. + PGICE*TGGI**4 + PGRND*TGGR**4) - TRNFLB(1) 3285. DO 440 L=1,LMA 3286. SRHR(I,J,L) = SRFHRL(L) 3287. 440 TRHR(I,J,L) = -TRFCRL(L) 3287.1 DO 450 K=1,4 3287.2 450 SRFRAC(I,J,K) = FSRNFG(K) 3288. C**** Update stratospheric temperatures above the dynamical top 3289. COSZ = COSZ2(I,J) 3290. DO 460 LR=1,3 3291. 460 RQT(I,J,LR) = RQT(I,J,LR) + 3292. * (SRFHRL(LMA+LR)*COSZ-TRFCRL(LMA+LR))*COE(LR) 3293. C**** 3294. C**** Accumulate diagnostics 3295. C**** 3296. C**** Longitudinal averages for budget pages 3297. DO 510 IT=1,6 3298. AF(J, 1,IT) = AF(J, 1,IT) + PT(IT)*COSZ*S0 3301. AF(J, 4,IT) = AF(J, 4,IT)+PT(IT)*COSZ*(SRNFLB(LMA+4)-SRNFLB(1)) 3305. AF(J,56,IT) = AF(J,56,IT)+PT(IT)*COSZ*(SRNFLB(LMA+1)-SRNFLB(1)) 3309. AF(J, 5,IT) = AF(J, 5,IT) + PT(IT)*COSZ*SRDFLB(1) 3317. AF(J,55,IT) = AF(J,55,IT) + PT(IT)*(BTEMPW-TKF) 3321. AF(J,67,IT) = AF(J,67,IT) + PT(IT)*TRDFLB(1) 3325. AF(J,70,IT) = AF(J,70,IT) - PT(IT)*(TRNFLB(LMA+4)-TRNFLB(1)) 3329. AF(J,71,IT) = AF(J,71,IT) - PT(IT)*(TRNFLB(LMA+1)-TRNFLB(1)) 3333. AF(J,72,IT) = AF(J,72,IT) + PT(IT)*COSZ*S0*PLAVIS 3337. AF(J,73,IT) = AF(J,73,IT) + PT(IT)*COSZ*S0*PLANIR 3341. AF(J,74,IT) = AF(J,74,IT) + PT(IT)*COSZ*S0*ALBVIS 3345. AF(J,75,IT) = AF(J,75,IT) + PT(IT)*COSZ*S0*ALBNIR 3349. AF(J,76,IT) = AF(J,76,IT) + PT(IT)*COSZ*S0*SRRVIS 3353. AF(J,77,IT) = AF(J,77,IT) + PT(IT)*COSZ*S0*SRRNIR 3357. AF(J,78,IT) = AF(J,78,IT) + PT(IT)*COSZ*S0*SRAVIS 3361. 510 AF(J,79,IT) = AF(J,79,IT) + PT(IT)*COSZ*S0*SRANIR 3365. C**** Latitude - layer diagnostics 3366. DO 520 L=1,LMA 3367. AJL(J,L,13) = AJL(J,L,13) + COSZ*SRFHRL(L) 3368. 520 AJL(J,L,14) = AJL(J,L,14) - TRFCRL(L) 3369. DO 530 LR=1,3 3370. AQJL(J,LR,3) = AQJL(J,LR,3) + COSZ*SRFHRL(LMA+LR) 3371. 530 AQJL(J,LR,4) = AQJL(J,LR,4) - TRFCRL(LMA+LR) 3372. C**** Longitude - latitude diagnostics 3373. AIJ(I,J,21) = AIJ(I,J,21) - (TRNFLB(LMA+4)-TRNFLB(1)) 3374. AIJ(I,J,24) = AIJ(I,J,24) + COSZ*SRNFLB(LMA+4) 3375. AIJ(I,J,25) = AIJ(I,J,25) + COSZ*S0 3376. AIJ(I,J,26) = AIJ(I,J,26) + COSZ*SRNFLB(1) 3377. AIJ(I,J,27) = AIJ(I,J,27) + COSZ*SRDFLB(1) 3380. C**** Diurnal cycle diagnostics 3381. DO 540 KR=1,4 3382. 540 IF(I.eq.IJDD(1,KR) .and. J.eq.IJDD(2,KR)) GO TO 550 3383. GO TO 600 3384. 550 JH=JHOUR 3385. DO 560 NR=1,NRAD 3386. IF(JH.ge.24) JH=JH-24 3387. ADAILY(JH,2,KR) = ADAILY(JH,2,KR) + (1.-SRNFLB(LMA+4)/S0) 3388. ADAILY(JH,3,KR) = ADAILY(JH,3,KR) + 3389. * (1.-SRNFLB(1)/(SRDFLB(1)+1.d-20)) 3390. ADAILY(JH,4,KR) = ADAILY(JH,4,KR) + (COSZ* 3391. * (SRNFLB(LMA+4)-SRNFLB(LMA+1))-TRNFLB(LMA+4)+TRNFLB(LMA+1)) 3392. 560 JH=JH+1 3393. C**** 3394. C**** End of outside DO loops over I and J 3395. C**** 3396. 600 CONTINUE 3397. 700 CONTINUE 3398. C**** 3399. C**** UPDATE THE TEMPERATURES BY RADIATION 3400. C**** 3401. 800 DO 820 J=1,JM 3402. IMAX=IM 3403. IF(J.eq.1 .or. J.eq.JM) IMAX=1 3404. DO 820 I=1,IMAX 3404.5 ATMABS = 0. 3405. DO 810 L=1,LMA 3406. H0M(I,J,L) = H0M(I,J,L) + (SRHR(I,J,L)*COSZ1(I,J)+TRHR(I,J,L))* 3407. * DTS*DXYP(J)*(1./PKDN(I,J,L)+1./PKUP(I,J,L))*.5 3407.1 HZM(I,J,L) = HZM(I,J,L) + (SRHR(I,J,L)*COSZ1(I,J)+TRHR(I,J,L))* 3407.2 * DTS*DXYP(J)*(1./PKUP(I,J,L)-1./PKDN(I,J,L)) 3407.3 810 ATMABS = ATMABS + (SRHR(I,J,L)*COSZ1(I,J)+TRHR(I,J,L)) 3408. C**** Accumulate diagnostics 3408.1 AJ(J,32) = AJ(J,32) + ATMABS* FGRND(I,J) 3408.2 BJ(J,32) = BJ(J,32) + ATMABS* FGICE(I,J) 3408.3 CJ(J,32) = CJ(J,32) + ATMABS* FLAKE(I,J)*(1.-RSI(I,J)) 3408.4 DJ(J,32) = DJ(J,32) + ATMABS* FLAKE(I,J)* RSI(I,J) 3408.5 EJ(J,32) = EJ(J,32) + ATMABS*FOCEAN(I,J)*(1.-RSI(I,J)) 3408.6 820 FJ(J,32) = FJ(J,32) + ATMABS*FOCEAN(I,J)* RSI(I,J) 3409. DO 830 KR=1,4 3410. I=IJDD(1,KR) 3410.1 J=IJDD(2,KR) 3411. 830 ADAILY(JHOUR,1,KR) = ADAILY(JHOUR,1,KR) + S0*COSZ1(I,J) 3413. RETURN 3414. C**** 3415. 903 FORMAT (' Low clouds in layers 1-',I2,' Mid level clouds in', 3416. * ' LAYERS',I3,'-',I2,' High clouds in layers',I3,'-',I2) 3417. END 4000. 4001. SUBROUTINE SURFCE 4002. C**** 4003. C**** SURFCE calculates the surface fluxes which include sensible 4004. C**** heat, evaporation, thermal radiation, and momentum drag. 4005. C**** These fluxes are applied to the atmospheric prognostic variables. 4006. C**** The composite surface temperature, surface specific humidity, 4007. C**** and surface wind components are saved in BLDATA. 4008. C**** 4009. INCLUDE 'C070.COM' 4010. C**** Parameters used by surface interaction 4011. C**** S1BYG1 = fraction of surface level from layer 1 center to ground 4012. C**** WSBYW1 = ratio of surface wind speed to layer 1 speed 4013. C**** CIAX = surface wind factor used to calculate cross isobar angle 4014. C**** WSADD = addition to surface wind speed used by surface fluxes 4015. C**** ALPHA = explicit to implicit fraction for surface fluxes 4016. C**** 4017. PARAMETER (S1BYG1=.57735d0, WSBYW1=.75, CIAX=.3d0, WSADD=2., 4018. * LBLM=3, ALPHA=1., 4019. * RHOW=1000., RHOI=910., RHOS=100., ALAMI=2.1762d0, ALAMS=.35d0, 4020. * BYRLI=1./(RHOI*ALAMI), BYRLS=1./(RHOS*ALAMS), 4023. * TKF=273.16d0, STBO=5.67032d-8) 4024. REAL*8 CDNL(IM,JM) 4025. LOGICAL*4 QPOLE,QDASRF,QDIAGD 4026. C**** 4027. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 4028. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 4029. * W0(IM,JM,4), E0(IM,JM,4),E1(IM,JM,4), UAS(IM,JM),VAS(IM,JM) 4031. COMMON /WORK00/ PKDN(IM,JM,LMA),PKUP(IM,JM,LMA),MSUM(IM,JM) 4032. COMMON /WORK01/ UAT(IM,JM,LMA), VAT(IM,JM,LMA) 4033. COMMON /WORK02/ is being used by subroutine COSZT 4033.5 COMMON /WORK03/ TGRND(IM,JM,4),TGRN2(IM,JM,4), COSZ1(IM,JM), 4034. * UVMS(IM+1), FTYPE(4),TCG(4),DEW(4),HDT(4),DMU(4),DMV(4), 4035. * H0MT(LBLM),HXMT(LBLM),HYMT(LBLM),HZMT(LBLM), RA(IM+1), 4036. * Q0MT(LBLM),QXMT(LBLM),QYMT(LBLM),QZMT(LBLM), ID(IM+1) 4037. C**** 4038. C**** FIXDCB FGRND Ground fraction (1) 4039. C**** FGICE Glacial ice fraction (1) 4040. C**** FWATER Water fraction (1) 4041. C**** 4041.1 C**** OCENCB G0M Mean enthalpy of ocean layers (J) 4041.2 C**** GZM Vertical gradient of ocean enthalpy (J) 4041.3 C**** 4042. C**** GRNDCB SBG Compressed snow depth of bare ground (m) 4043. C**** SVG Compressed snow depth of vegetated ground (m) 4043.9 C**** 4044. C**** GICECB MGI Mass of glacial ice (kg/m^2) 4045. C**** HGI Enthalpy minus latent heat of glacial ice (J/m^2) 4046. C**** 4050. C**** SEAICB RSI Horizontal ratio of sea ice to ocean (1) 4051. C**** MSI Mass of sea ice (kg/m^2) 4052. C**** HSI Enthalpy minus latent heat of sea ice (J/m^2) 4053. C**** 4054. C**** SURFCB BLDATA(1) Ground temperature over ground fraction (C) 4055. C**** BLDATA(2) Ground wetness over ground fraction 4056. C**** BLDATA(3) Composite surface wind speed (m/s) 4057. C**** 4058. C**** FLUXCB COSZ1 Average zenith angle during source time step 4059. C**** SROX Solar radiation absorbed by open water (J/m^2) 4060. C**** DMUA U momentum downward from atmosphere (kg/m*s) 4061. C**** DMVA V momentum downward from atmosphere (kg/m*s) 4062. C**** W0 Water flux downward into surface (kg/m^2) 4063. C**** E0 Energy downward into surface (J/m^2) 4064. C**** E1 Energy from ground layer 1 to layer 2 (J/m^2) 4064.5 C**** UAS U atmospheric surface wind (m/s) on A grid 4064.6 C**** VAS V atmospheric surface wind (m/s) on A grid 4065. C**** 4066. COMMON /SOILS2/ SDATA(IM,JM,LMG*11+1) 4067. COMMON /SOILPR/ SNOWBG,SNOWVG, WATRBG(0:LMG),WATRVG(0:LMG), 4068. * HEATBG(0:LMG),HEATVG(0:LMG) 4069. COMMON /SOILIN/ PRECTO,PRECLS,EPRETO,EPRELS, PG100,TKS,QS, 4070. * CDM,WSWSAD,RHOA, SRH,TRHDN, COSYR,SINYR, SDUMMY(5), IGCM 4071. COMMON /SOILDO/ ACNA,ACNC, 4072. * ABETA,ABETAB,ABETAV,ABETAP,ABETAT,ABETAD, 4073. * AEPB,AEPC,AEPP, AEVAPB,AEVAPD,AEVAPW, ADIFS,ARUNS,ARUNU, 4074. * ATRG,ASHG,ALHG, AEDIFS,AF0DT,AF1DT,AFHG, AERUNS,AERUNU 4075. C**** 4076. C**** ACNA AERODYNAMIC CONDUCTANCE (m/s) 4077. C**** ACNC CANOPY CONDUCTANCE (m/s) 4078. C**** ABETA OVERALL EVAPORATION EFFICIENCY 4079. C**** ABETAB BARE SOIL EVAPORATION EFFICIENCY 4080. C**** ABETAV CANOPY EVAPORATION EFFICIENCY 4081. C**** ABETAP PENMAN EVAPORATION EFFICIENCY 4082. C**** ABETAT TRANSPIRATION EFFICIENCY 4083. C**** ABETAD ROOT BETA FOR TRANSPIRATION 4084. C**** AEPB POTENTIAL EVAP FROM BARE SOIL (kg/m^2) 4085. C**** AEPC POTENTIAL EVAP FROM CANOPY (kg/m^2) 4086. C**** AEPP PENMAN POTENTIAL EVAP (kg/m^2) 4087. C**** AEVAPB EVAPORATION FROM BARE SOIL (kg/m^2) 4088. C**** AEVAPD EVAPORATION FROM DRY CANOPY (kg/m^2) 4089. C**** AEVAPW EVAPORATION FROM WET CANOPY (kg/m^2) 4090. C**** ADIFS WATER FLUX FROM TOP LAYER TO SECOND (kg/m^2) 4091. C**** ARUNS SURFACE RUNOFF (kg/m^2) 4092. C**** ARUNU UNDERGROUND RUNOFF (kg/m^2) 4093. C**** ATRG THERMAL RADIATION FROM THE LAND SURFACE (J M-2) 4094. C**** ASHG SENSIBLE HEAT FROM THE LAND SURFACE (J M-2) 4095. C**** ALHG LATENT HEAT FROM THE LAND SURFACE (J M-2) 4096. C**** AEDIFS HEAT CARRIED BY WATER FROM TOP LAYER (J M-2) 4097. C**** AF0DT HEAT FLUX INTO LAND SURFACE (J M-2) 4098. C**** AF1DT HEAT FLUX INTO GROUND (J M-2) 4099. C**** AFHG HEAT FLUX FROM GROUND TO CANOPY (J M-2) 4100. C**** AERUNS HEAT CARRIED BY SURFACE RUNOFF (J M-2) 4101. C**** AERUNU HEAT CARRIED BY UNDERGROUND RUNOFF (J M-2) 4102. C**** 4103. DATA IFIRST/1/, JDLAST/-9/ 4104. C**** QSAT = (RAIR/RVAP)*610.71*EXP((L/RVAP)*(1/TKF-1/T)) / PRESS 4105. QSAT(TM,PR) = 379.7915d0*EXP(ELHX*(1./TKF-1./TM)/RVAP)/PR 4106. IF(IFIRST.ne.1) GO TO 10 4107. IFIRST=0 4108. DTSURF = DTS/NSURF 4109. BYDTS = 1./DTS 4110. C**** Read in neutral drag coefficient for ground and glacial ice 4111. CALL READR4 (53,IM*JM,CDNL,CDNL) 4112. CLOSE (53) 4113. C**** Read in soils data for ground hydrology 4113.5 DO 5 L=1,LMG*11+1 4114. 5 CALL READR4 (52,IM*JM,SDATA(1,1,L),SDATA(1,1,L)) 4115. CLOSE (52) 4116. C**** Initialize ground hydrology 4117. IGCM = -1 4118. CALL GHINIT (FLAKE,FVEG,DTSURF,DUMMY) 4119. C**** Calculate sine and cosine of time of year foe each new day 4120. 10 IF(JDAY.eq.JDLAST) GO TO 20 4121. JDLAST = JDAY 4122. COSYR = COS (TWOPI*(JDAY-.5)/365.) 4123. SINYR = SIN (TWOPI*(JDAY-.5)/365.) 4124. C**** Zero out surface energy and evaporation and initialize TGRND 4125. 20 DO 40 I=1,IM*JM 4125.5 IF(RSI(I,1).gt.0.) then 4126. TGRND(I,1,2) = (HSI(I,1,1)/(XSI1*MSI(I,1,1)) + ELHM) / SHCI 4127.1 TGRN2(I,1,2) = (HSI(I,1,2)/(XSI2*MSI(I,1,1)) + ELHM) / SHCI 4127.2 endif 4127.5 IF(FGICE(I,1).gt.0.) then 4127.6 TGRND(I,1,3) = (HGI(I,1,1)/(XGI1*MGI(I,1)) + ELHM) / SHCI 4127.7 TGRN2(I,1,3) = (HGI(I,1,2)/(XGI2*MGI(I,1)) + ELHM) / SHCI 4127.8 endif 4127.9 40 continue 4128. DO 50 I=1,IM*JM*17 4129. 50 SROX(I,1) = 0. 4130. C**** 4131. C**** Outside loop over surface time steps 4132. C**** 4133. DO 920 NS=0,NSURF-1 4134. QDASRF = MOD((IDT-IDT0)*NSURF+NS,NDIAGS).eq.0 4135. IF(QDASRF) IDACC(3) = IDACC(3) + 1 4136. QDIAGD = MOD(IDAY+NS,NSURF).eq.0 4137. C**** Copy UA,VA to UAT,VAT. The former arrays will be updated by 4138. C**** surface friction, while the latter arrays will be fixed. 4139. DO 60 I=1,IM*JM*(LBLM+1) 4140. UAT(I,1,1) = UA(I,1,1) 4141. 60 VAT(I,1,1) = VA(I,1,1) 4142. C**** Calculate cosine of zenith angle 4143. IF(NSURF.le.1) GO TO 100 4144. ROT1 = (NS + NSURF*MOD(IDT,NDAY))*TWOPI / (NDAY*NSURF) 4145. ROT2 = (NS+1 + NSURF*MOD(IDT,NDAY))*TWOPI / (NDAY*NSURF) 4146. CALL COSZT (SIND,COSD,ROT1,ROT2,COSZ1) 4147. C**** 4148. C**** Outside loop over J 4149. C**** 4150. 100 DO 910 J=1,JM 4151. HEMI = 1. 4152. IF(J.le.JM/2) HEMI = -1. 4153. QPOLE = J.eq.1 .OR. J.eq.JM 4154. IF(QPOLE) GO TO 110 4155. C**** Conditions at non-polar points 4156. IMAX=IM 4157. KMAX=4 4158. RA(1) = .5 4159. RA(2) = .5 4160. RA(3) = RAMVS(J) 4161. RA(4) = RAMVN(J) 4162. GO TO 150 4163. C**** Conditions at the poles 4164. 110 IMAX=1 4165. JVPO=1 4166. IF(J.eq.JM) JVPO=JM-1 4167. U1 = 0. 4168. V1 = 0. 4169. DO 120 I=1,IM 4170. U1 = U1 - VAT(I,JVPO,1)*SINI(I) 4171. 120 V1 = V1 + VAT(I,JVPO,1)*COSI(I) 4172. U1 = U1*HEMI*2./IM 4173. V1 = V1*2./IM 4174. DO 130 I=1,IM 4175. RA(I) = RAMVN(1) 4176. 130 ID(I) = I + (JVPO-1)*IM + IM*JM*LMA 4177. RA(IM+1) = 1. 4178. ID(IM+1) = 1 + (J-1)*IM 4179. C**** 4180. C**** Outside loop over I 4181. C**** 4182. 150 IM1=IM 4183. DO 900 I=1,IMAX 4184. IF(QPOLE) GO TO 200 4185. C**** Interpolate velocities to primary grid away from poles 4186. ID(1) = IM1+(J-1)*IM 4187. ID(2) = I +(J-1)*IM 4188. ID(4) = ID(2)+IM*JM*LMA 4189. ID(3) = ID(4)-IM 4190. U1 = .5*(UAT(IM1,J,1)+UAT(I,J,1)) 4191. V1 = .5*(VAT(I,J-1,1)+VAT(I,J,1)) 4192. C**** 4193. C**** Determine surface conditions 4194. C**** 4195. 200 FTYPE(1) = FWATER(I,J)*(1.-RSI(I,J)) 4196. FTYPE(2) = FWATER(I,J)*RSI(I,J) 4197. FTYPE(3) = FGICE(I,J) 4198. FTYPE(4) = FGRND(I,J) 4199. PS = GRAV*(MSUM(I,J)-MA(I,J,1)*.5*(1.-S1BYG1)) 4200. PG = GRAV* MSUM(I,J) 4202. PGK = PG**RKAP 4203. HS = (H0M(I,J,1)-HZM(I,J,1)*S1BYG1)*PGK/(DXYP(J)*MA(I,J,1)) 4204. Q1 = Q0M(I,J,1)/(DXYP(J)*MA(I,J,1)) 4205. QS = (Q0M(I,J,1)-QZM(I,J,1)*S1BYG1)/(DXYP(J)*MA(I,J,1)) 4206. BETAS = RKAP*HS/PG 4207. IF(QS.ge.0.) GO TO 210 4208. QZM(I,J,1) = Q0M(I,J,1)/S1BYG1 4209. QS = 0. 4210. 210 WSSQ = (U1*U1 + V1*V1)*WSBYW1*WSBYW1 4211. WS = SQRT(WSSQ) 4213.1 C**** Calculate atmospheric variables used for diagnostics 4213.2 QA1DN = (Q0M(I,J,1)-.5*QZM(I,J,1)) / (DXYP(J)*MA(I,J,1)) 4213.3 QA1UP = (Q0M(I,J,1)+.5*QZM(I,J,1)) / (DXYP(J)*MA(I,J,1)) 4213.4 TA1DN = (H0M(I,J,1)-.5*HZM(I,J,1))*PGK / 4213.5 * (DXYP(J)*MA(I,J,1)* SHCD ) - TKF 4213.6 TA1UP = (H0M(I,J,1)+.5*HZM(I,J,1))*PGK / 4213.7 * (DXYP(J)*MA(I,J,1)* SHCD ) - TKF 4214. C**** Zero out quantities to be summed over surface types 4214.5 DEWS = 0. 4215. USS = 0. 4216. VSS = 0. 4217. TCSS = 0. 4218. QSS = 0. 4219. DMUS = 0. 4220. DMVS = 0. 4221. TCGS = 0. 4222. QGS = 0. 4222.5 SRHDTS= 0. 4223. TRHDTS= 0. 4224. SNHDTS= 0. 4225. EVHDTS= 0. 4226. CIAS = 0. 4227. RIGSS = 0. 4228. CDMS = 0. 4229. C CDHS = 0. 4230. TKES = 0. 4231. C**** 4232. IF(FTYPE(1).le.0.) GO TO 220 4233. C**** 4234. C**** Open Ocean or Open Lake 4235. C**** 4236. ITYPE = 1 4237. PTYPE = FTYPE(1) 4238. CDN = .0009d0 + .0000025d0*WSSQ 4238.5 SRH = SRFRAC(I,J,1)*COSZ1(I,J) 4238.6 SRHDT = SRH*DTSURF 4239. IF(FOCEAN(I,J).gt.0.) then 4240. GO1 = G0M(I,J,1)/(MO(I,J,1)*DXYP(J)) 4241. SO1 = S0M(I,J,1)/(MO(I,J,1)*DXYP(J)) 4242. TCG(1) = TEMGS(GO1,SO1) 4243. else 4244. TCG(1) = GZM(I,J,1) 4245. endif 4246. SROX(I,J) = SROX(I,J) + SRHDT 4247. ELHX = ELHE 4248. GO TO 300 4249. C**** 4250. 220 IF(FTYPE(2).le.0.) GO TO 240 4251. C**** 4252. C**** Sea Ice or Lake Ice 4253. C**** 4254. ITYPE = 2 4255. PTYPE = FTYPE(2) 4256. SNOW = MSI(I,J,1) - ACE1I 4258. TCG(2)= TGRND(I,J,2) 4259. TCG2 = TGRN2(I,J,2) 4264. dF1dTG= 2. / (ACE1I*BYRLI + SNOW*BYRLS) 4265. CDN = .00095417d0 + .0000005d0*ACE1I + .0000005d0*MSI(I,J,2) 4266. HCG1 = SHCI*XSI1*MSI(I,J,1) 4266.1 HCG2 = SHCI*XSI2*MSI(I,J,1) 4267. ELHX = ELHS 4267.5 SRH = SRFRAC(I,J,2)*COSZ1(I,J) 4267.6 SRHDT = SRH*DTSURF 4268. GO TO 300 4269. C**** 4270. 240 IF(FTYPE(3)+FTYPE(4).le.0.) GO TO 600 4271. CDN = CDNL(I,J) 4272. IF(FTYPE(3).le.0.) GO TO 260 4273. C**** 4274. C**** Glacial Ice 4275. C**** 4276. ITYPE = 3 4277. PTYPE = FTYPE(3) 4278. SNOW = MGI(I,J) - ACE1I 4279. TCG(3)= TGRND(I,J,3) 4280. TCG2 = TGRN2(I,J,3) 4283. dF1dTG= 2. / (ACE1I*BYRLI + SNOW*BYRLS) 4284. HCG1 = SHCI*XGI1*MGI(I,J) 4284.1 HCG2 = SHCI*XGI2*MGI(I,J) 4285. ELHX = ELHS 4285.5 SRH = SRFRAC(I,J,3)*COSZ1(I,J) 4285.6 SRHDT = SRH*DTSURF 4286. GO TO 300 4287. C**** 4288. 260 IF(FTYPE(4).le.0.) GO TO 600 4289. C**** 4290. C**** Ground 4291. C**** 4292. ITYPE = 4 4293. PTYPE = FTYPE(4) 4294. FBARE = FVEG(I,J,1) + FVEG(I,J,10) 4295. TCG(4) = BLDATA(I,J,1) 4296. SNOWBG = SBG(I,J) 4297. SNOWVG = SVG(I,J) 4298. DO 270 L=1,LMG 4299. WATRBG(L) = WBG(I,J,L) 4300. WATRVG(L) = WVG(I,J,L) 4301. HEATBG(L) = HBG(I,J,L) 4302. 270 HEATVG(L) = HVG(I,J,L) 4303. WATRVG(0) = WVG(I,J,0) 4304. HEATVG(0) = HVG(I,J,0) 4305. PRECTO = PREC(I,J,1)*BYDTS/RHOW 4306. PRECLS = PREC(I,J,2)*BYDTS/RHOW 4307. EPRETO = EPRE(I,J,1)*BYDTS 4308. EPRELS = EPRE(I,J,2)*BYDTS 4309. PG100 = PG/100. 4310. TKS = HS/SHCD 4311. WSWSAD = WS + WSADD 4312. RHOA = 1./BETAS 4312.5 SRH = SRFRAC(I,J,4)*COSZ1(I,J) 4312.6 SRHDT = SRH*DTSURF 4313. TRHDN = TRHR(I,J,0) 4313.5 C AMASS1 = MA(I,J,1) 4314. C**** Accumulate ground diagnostics 4315. IF(NS.gt.0) GO TO 300 4316. XBARE = RHOW*PTYPE*FBARE 4317. XVEGE = RHOW*PTYPE*(1.-FBARE) 4318. AJ(J,18) = AJ(J,18) + PTYPE*TCG(4) 4319. AJ(J,49) = AJ(J,49) + XVEGE*WATRVG(0) 4320. AJ(J,50) = AJ(J,50) + XBARE*WATRBG(1) + XVEGE*WATRVG(1) 4321. AJ(J,51) = AJ(J,51) + XBARE*WATRBG(2) + XVEGE*WATRVG(2) 4323. AIJ(I,J,17) = AIJ(I,J,17) + FBARE* 4324. * (WATRBG(1)+WATRBG(2)+WATRBG(3)+WATRBG(4)+WATRBG(5)+WATRBG(6)) 4325. * + (1.-FBARE)* 4326. * (WATRVG(1)+WATRVG(2)+WATRVG(3)+WATRVG(4)+WATRVG(5)+WATRVG(6)) 4327. IF(XBARE*SNOWBG+XVEGE*SNOWVG .le. 0.) GO TO 300 4328. AJ(J,31) = AJ(J,31) + PTYPE 4329. AJ(J,49) = AJ(J,50) - XBARE*SNOWBG - XVEGE*SNOWVG 4330. AJ(J,53) = AJ(J,53) + XBARE*SNOWBG + XVEGE*SNOWVG 4331. AIJ(I,J,2) = AIJ(I,J,2) + PTYPE 4332. AIJ(I,J,3) = AIJ(I,J,3) + XBARE*SNOWBG + XVEGE*SNOWVG 4333. C**** 4334. C**** Calaculate surface layer parameters, drag coefficients, and 4335. C**** eddy diffusion coefficient after calling the PBL table. 4336. C**** 4337. C**** Define ground variables and branch on stability 4338. 300 TKG = TCG(ITYPE) + TKF 4339. BETAG = RKAP*TKG*SHCD/PG 4340. IF(TKG*SHCD.lt.HS) GO TO 350 4341. C**** 4342. C**** TKG*SHCD > HS: boundary layer has unstable stratification 4343. C**** 4344. RIGS = 0. 4345. CDM = CDN*2. 4346. C CDH = CDM 4347. CIA = TWOPI*.0625/(1.+WS*CIAX) 4348. GO TO 400 4349. C**** 4350. C**** TKG*SHCD < HS: atmosphere is stable with respect to the ground 4351. C**** 4352. 350 RIGS = (PG-PS)*(BETAS-BETAG) / (WSSQ+.1) 4353. CDM = CDN*2. / (1.+4.*RIGS*RIGS) 4354. C CDH = CDM 4355. CIA = TWOPI*(.09375-.03125/(1.+4.*RIGS*RIGS)) / (1.+WS*CIAX) 4356. C**** 4357. 400 COSC = COS(CIA) 4358. SINC = SIN(CIA)*HEMI 4359. US = (U1*COSC - V1*SINC)*WSBYW1 4360. VS = (V1*COSC + U1*SINC)*WSBYW1 4361. C**** 4362. C**** Calculate fluxes of sensible heat, latent heat, thermal 4363. C**** radiation, and ground conduction (W/m^2) 4364. C**** 4365. CDMWSB = CDM*(WS+WSADD)/BETAS 4366. IF(ITYPE.eq.4) GO TO 450 4367. QG = QSAT(TKG,PG) 4368. SNH = (HS-TKG*SHCD)*CDMWSB 4369. C EVH = (QS-QG)*(ELHE+TCG(ITYPE)*SHCV)*CDMWSB 4370. EVH = (QS-QG)*(ELHE )*CDMWSB 4371. TRH = TRHR(I,J,0)-STBO*(TKG*TKG)*(TKG*TKG) 4372. IF(ITYPE.eq.1) GO TO 420 4373. C**** Calculate fluxes using implicit time step for nonocean points 4374. F0 = SRH+TRH+SNH+EVH 4375. F1 = (TCG(ITYPE)-TCG2)*dF1dTG 4377. DSNDTG = -SHCD*CDMWSB 4378. DQGDTG = QG*ELHX/(RVAP*TKG*TKG) 4379. DEVDTG = -DQGDTG*ELHE*CDMWSB 4380. DTRDTG = -4.*STBO*TKG*TKG*TKG 4381. DF0DTG = DSNDTG+DEVDTG+DTRDTG 4383. C DSNDHS = CDMWSB 4384. DEVDQS = ELHE*CDMWSB 4385. HSDEN = (1.+2.*S1BYG1)*ALPHA*DTSURF*PGK*CDMWSB + 4385.1 * MA(I,J,1)*PKDN(I,J,1) 4386. QSDEN = (1.+2.*S1BYG1)*ALPHA*DTSURF*DEVDQS + MA(I,J,1)*ELHE 4387. HSCON = -(1.+2.*S1BYG1)*DTSURF*PGK*SNH / HSDEN 4388. QSCON = -(1.+2.*S1BYG1)*DTSURF*EVH / QSDEN 4389. HSMUL = -(1.+2.*S1BYG1)*ALPHA*DTSURF*PGK*DSNDTG / HSDEN 4390. QSMUL = -(1.+2.*S1BYG1)*ALPHA*DTSURF*DEVDTG / QSDEN 4390.5 T2DEN = HCG2 + ALPHA*DTSURF*dF1dTG 4390.6 T2CON = DTSURF*F1 / T2DEN 4390.7 T2MUL = ALPHA*DTSURF*dF1dTG / T2DEN 4391. TGDEN = HCG1 - ALPHA*DTSURF* 4392. * (DF0DTG-dF1dTG + HSMUL*CDMWSB + QSMUL*DEVDQS + T2MUL*dF1dTG) 4393. DTG = DTSURF*(F0-F1 + 4393.1 + ALPHA*(HSCON*CDMWSB+QSCON*DEVDQS+T2CON*DF1DTG)) / TGDEN 4394. IF(TCG(ITYPE)+DTG.gt.0.) DTG = - TCG(ITYPE) 4395. DHS = HSCON + HSMUL*DTG 4396. DQS = QSCON + QSMUL*DTG 4396.5 DT2 = T2CON + T2MUL*DTG 4397. SNHDT = DTSURF*(SNH + ALPHA*(DTG*DSNDTG+DHS*CDMWSB)) 4398. EVHDT = DTSURF*(EVH + ALPHA*(DTG*DEVDTG+DQS*DEVDQS)) 4399. TRHDT = DTSURF*(TRH + ALPHA* DTG*DTRDTG) 4400. F1DT = DTSURF*(F1 + ALPHA*(DTG*dF1dTG-DT2*dF1dTG)) 4401. TCG(ITYPE) = TCG(ITYPE) + DTG 4402. TGRN2(I,J,ITYPE) = TCG2 + DT2 4404. GO TO 430 4405. C**** Calculate fluxes using explicit time step for ocean points 4406. 420 DEVDQS= ELHE*CDMWSB 4407. DHS = -(1.+2.*S1BYG1)*DTSURF*PGK*SNH / 4408. * ((1.+2.*S1BYG1)*DTSURF*PGK*CDMWSB + MA(I,J,1)*PKDN(I,J,1)) 4409. DQS = -(1.+2.*S1BYG1)*DTSURF*EVH / 4410. * ((1.+2.*S1BYG1)*DTSURF*DEVDQS + MA(I,J,1)*ELHE) 4411. SNHDT = DTSURF*(SNH+DHS*CDMWSB) 4412. EVHDT = DTSURF*(EVH+DQS*DEVDQS) 4413. TRHDT = DTSURF*TRH 4414. C**** Limit amount of dew 4415. C 430 DEW(ITYPE) = EVHDT/(ELHE+TCG(ITYPE)*SHCV) 4416. 430 DEW(ITYPE) = EVHDT/ ELHE 4417. IF(DEW(ITYPE).le..5*Q1*MA(I,J,1)) GO TO 440 4417.5 DEW(ITYPE) = .5*Q1*MA(I,J,1) 4418. EVHDT0= EVHDT 4419. C EVHDT = DEW(ITYPE)*(ELHE+TCG(ITYPE)*SHCV) 4420. EVHDT = DEW(ITYPE)* ELHE 4421. IF(ITYPE.ne.1) TCG(ITYPE) = TCG(ITYPE) + (EVHDT-EVHDT0)/HCG1 4422. C**** Store surface fluxes into arrays for subsequent subroutines 4423. 440 E0(I,J,ITYPE) = E0(I,J,ITYPE) + (SRHDT+TRHDT+SNHDT+EVHDT) 4424. E1(I,J,ITYPE) = E1(I,J,ITYPE) + F1DT 4425. W0(I,J,ITYPE) = W0(I,J,ITYPE) + DEW(ITYPE) 4425.5 TGRND(I,J,ITYPE) = TCG(ITYPE) 4426. GO TO 470 4427. C**** Ground hydrology routines return surface fluxes and update 4428. C**** ground prognostic variables 4429. 450 CALL GHIJ (I,J,TCG(4),WEARTH) 4429.1 C EVAP = AEVAPB+AEVAPD+AEVAPW 4429.2 C ALHG = ALHG + SHCV*TCG(4)*EVAP 4429.3 C AERUNU = AERUNU - SHCV*TCG(4)*EVAP 4430. SBG(I,J) = SNOWBG 4431. SVG(I,J) = SNOWVG 4432. DO 460 L=1,LMG 4433. WBG(I,J,L) = WATRBG(L) 4434. WVG(I,J,L) = WATRVG(L) 4435. HBG(I,J,L) = HEATBG(L) 4436. 460 HVG(I,J,L) = HEATVG(L) 4437. WVG(I,J,0) = WATRVG(0) 4438. HVG(I,J,0) = HEATVG(0) 4439. TRHDT = DTSURF*TRHDN - ATRG 4440. SNHDT = - ASHG 4441. EVHDT = - ALHG 4442. DEW(ITYPE) = - (AEVAPB+AEVAPD+AEVAPW) 4443. C**** 4445. 470 HDT(ITYPE) = SNHDT 4446. DMU(ITYPE) = DTSURF*CDMWSB*US 4447. DMV(ITYPE) = DTSURF*CDMWSB*VS 4448. SSX = WS*MA(I,J,1) / 4449. * (WS*MA(I,J,1) + SQRT(DMU(ITYPE)**2+DMV(ITYPE)**2)) 4450. DMU(ITYPE) = DMU(ITYPE)*SSX 4451. DMV(ITYPE) = DMV(ITYPE)*SSX 4452. C**** Sum quantities over surface types 4453. C * TRHDT,SNHDT,EVHDT 4454. DEWS = DEWS + PTYPE*DEW(ITYPE) 4455. USS = USS + PTYPE*US 4456. VSS = VSS + PTYPE*VS 4457. TCS = HS/SHCD - TKF 4458. TCSS = TCSS + PTYPE*TCS 4459. QSS = QSS + PTYPE*QS 4460. DMUS = DMUS + PTYPE*DMU(ITYPE) 4461. DMVS = DMVS + PTYPE*DMV(ITYPE) 4462. TCGS = TCGS + PTYPE*TCG(ITYPE) 4463. QGS = QGS + PTYPE*QG 4463.5 SRHDTS= SRHDTS+ PTYPE*SRHDT 4464. TRHDTS= TRHDTS+ PTYPE*TRHDT 4465. SNHDTS= SNHDTS+ PTYPE*SNHDT 4466. EVHDTS= EVHDTS+ PTYPE*EVHDT 4467. CIAS = CIAS + PTYPE*CIA 4468. RIGSS = RIGSS + PTYPE*RIGS 4469. CDMS = CDMS + PTYPE*CDM 4470. C CDHS = CDHS + PTYPE*CDH 4470.1 TKES = TKES + PTYPE*(US*DMU(ITYPE) + VS*DMV(ITYPE)) 4471. C**** 4472. GO TO (500,520,540,560),ITYPE 4473. C**** 4474. C**** OCEAN 4475. C**** 4476. 500 DMUA(I,J,1) = DMUA(I,J,1) + DMU(ITYPE) 4477. DMVA(I,J,1) = DMVA(I,J,1) + DMV(ITYPE) 4477.5 IF(FLAKE(I,J).gt.0.) then 4477.6 CJ(J, 6) = CJ(J, 6) + FLAKE(I,J)*(1.-RSI(I,J))*SRHDT 4478. CJ(J, 9) = CJ(J, 9) + FLAKE(I,J)*(1.-RSI(I,J))*TRHDT 4480. CJ(J,13) = CJ(J,13) + FLAKE(I,J)*(1.-RSI(I,J))*SNHDT 4481. CJ(J,14) = CJ(J,14) + FLAKE(I,J)*(1.-RSI(I,J))*EVHDT 4482. C CJ(J,19) = CJ(J,19) + FLAKE(I,J)*(1.-RSI(I,J))*DEW(1) 4482.1 C CJ(J,22) = CJ(J,22) + FLAKE(I,J)*(1.-RSI(I,J))*DEW(1)*ZATMO( 4483. CJ(J,23) = CJ(J,23) + FLAKE(I,J)*(1.-RSI(I,J))*TCS 4483.1 CJ(J,26) = CJ(J,26) + FLAKE(I,J)*(1.-RSI(I,J))*QA1UP 4483.2 CJ(J,27) = CJ(J,27) + FLAKE(I,J)*(1.-RSI(I,J))*QA1DN 4483.3 CJ(J,28) = CJ(J,28) + FLAKE(I,J)*(1.-RSI(I,J))*QS 4483.4 CJ(J,29) = CJ(J,29) + FLAKE(I,J)*(1.-RSI(I,J))*QG 4483.5 CJ(J,35) = CJ(J,35) + FLAKE(I,J)*(1.-RSI(I,J))*CIA 4483.6 CJ(J,43) = CJ(J,43) + FLAKE(I,J)*(1.-RSI(I,J))*TA1UP 4483.7 CJ(J,44) = CJ(J,44) + FLAKE(I,J)*(1.-RSI(I,J))*TA1DN 4483.8 CJ(J,48) = CJ(J,48) + FLAKE(I,J)*(1.-RSI(I,J))*WS 4483.9 CJ(J,64) = CJ(J,64) + FLAKE(I,J)*(1.-RSI(I,J))*BETAS 4483.91 CJ(J,65) = CJ(J,65) + FLAKE(I,J)*(1.-RSI(I,J))*RIGS 4483.92 CJ(J,66) = CJ(J,66) + FLAKE(I,J)*(1.-RSI(I,J))*CDM 4483.95 endif 4483.96 IF(FOCEAN(I,J).gt.0.) then 4483.97 EJ(J, 6) = EJ(J, 6) + FOCEAN(I,J)*(1.-RSI(I,J))*SRHDT 4484. EJ(J, 9) = EJ(J, 9) + FOCEAN(I,J)*(1.-RSI(I,J))*TRHDT 4486. EJ(J,13) = EJ(J,13) + FOCEAN(I,J)*(1.-RSI(I,J))*SNHDT 4487. EJ(J,14) = EJ(J,14) + FOCEAN(I,J)*(1.-RSI(I,J))*EVHDT 4488. C EJ(J,19) = EJ(J,19) + FOCEAN(I,J)*(1.-RSI(I,J))*DEW(1) 4488.1 C EJ(J,22) = EJ(J,22) + FOCEAN(I,J)*(1.-RSI(I,J))*DEW(1)*ZATMO( 4489. EJ(J,23) = EJ(J,23) + FOCEAN(I,J)*(1.-RSI(I,J))*TCS 4489.1 EJ(J,26) = EJ(J,26) + FOCEAN(I,J)*(1.-RSI(I,J))*QA1UP 4489.2 EJ(J,27) = EJ(J,27) + FOCEAN(I,J)*(1.-RSI(I,J))*QA1DN 4489.3 EJ(J,28) = EJ(J,28) + FOCEAN(I,J)*(1.-RSI(I,J))*QS 4489.4 EJ(J,29) = EJ(J,29) + FOCEAN(I,J)*(1.-RSI(I,J))*QG 4489.5 EJ(J,35) = EJ(J,35) + FOCEAN(I,J)*(1.-RSI(I,J))*CIA 4489.6 EJ(J,43) = EJ(J,43) + FOCEAN(I,J)*(1.-RSI(I,J))*TA1UP 4489.7 EJ(J,44) = EJ(J,44) + FOCEAN(I,J)*(1.-RSI(I,J))*TA1DN 4489.8 EJ(J,48) = EJ(J,48) + FOCEAN(I,J)*(1.-RSI(I,J))*WS 4489.9 EJ(J,64) = EJ(J,64) + FOCEAN(I,J)*(1.-RSI(I,J))*BETAS 4489.91 EJ(J,65) = EJ(J,65) + FOCEAN(I,J)*(1.-RSI(I,J))*RIGS 4489.92 EJ(J,66) = EJ(J,66) + FOCEAN(I,J)*(1.-RSI(I,J))*CDM 4489.95 endif 4490. GO TO 220 4491. C**** 4492. C**** OCEAN ICE 4493. C**** 4494. 520 DMUA(I,J,2) = DMUA(I,J,2) + DMU(ITYPE) 4495. DMVA(I,J,2) = DMVA(I,J,2) + DMV(ITYPE) 4495.5 IF(FLAKE(I,J).gt.0.) then 4495.6 DJ(J, 6) = DJ(J, 6) + FLAKE(I,J)*RSI(I,J)*SRHDT 4496. DJ(J, 9) = DJ(J, 9) + FLAKE(I,J)*RSI(I,J)*TRHDT 4498. DJ(J,13) = DJ(J,13) + FLAKE(I,J)*RSI(I,J)*SNHDT 4499. DJ(J,14) = DJ(J,14) + FLAKE(I,J)*RSI(I,J)*EVHDT 4500. C DJ(J,19) = DJ(J,19) + FLAKE(I,J)*RSI(I,J)*DEW(2) 4500.1 C DJ(J,22) = DJ(J,22) + FLAKE(I,J)*RSI(I,J)*DEW(2)*ZATMO(I,J) 4501. DJ(J,23) = DJ(J,23) + FLAKE(I,J)*RSI(I,J)*TCS 4501.1 DJ(J,26) = DJ(J,26) + FLAKE(I,J)*RSI(I,J)*QA1UP 4501.2 DJ(J,27) = DJ(J,27) + FLAKE(I,J)*RSI(I,J)*QA1DN 4501.3 DJ(J,28) = DJ(J,28) + FLAKE(I,J)*RSI(I,J)*QS 4501.4 DJ(J,29) = DJ(J,29) + FLAKE(I,J)*RSI(I,J)*QG 4501.5 DJ(J,35) = DJ(J,35) + FLAKE(I,J)*RSI(I,J)*CIA 4501.6 DJ(J,43) = DJ(J,43) + FLAKE(I,J)*RSI(I,J)*TA1UP 4501.7 DJ(J,44) = DJ(J,44) + FLAKE(I,J)*RSI(I,J)*TA1DN 4501.8 DJ(J,48) = DJ(J,48) + FLAKE(I,J)*RSI(I,J)*WS 4501.9 DJ(J,64) = DJ(J,64) + FLAKE(I,J)*RSI(I,J)*BETAS 4501.91 DJ(J,65) = DJ(J,65) + FLAKE(I,J)*RSI(I,J)*RIGS 4501.92 DJ(J,66) = DJ(J,66) + FLAKE(I,J)*RSI(I,J)*CDM 4501.95 endif 4501.96 IF(FOCEAN(I,J).gt.0.) then 4501.97 FJ(J, 6) = FJ(J, 6) + FOCEAN(I,J)*RSI(I,J)*SRHDT 4502. FJ(J, 9) = FJ(J, 9) + FOCEAN(I,J)*RSI(I,J)*TRHDT 4504. FJ(J,13) = FJ(J,13) + FOCEAN(I,J)*RSI(I,J)*SNHDT 4505. FJ(J,14) = FJ(J,14) + FOCEAN(I,J)*RSI(I,J)*EVHDT 4506. C FJ(J,19) = FJ(J,19) + FOCEAN(I,J)*RSI(I,J)*DEW(2) 4506.1 C FJ(J,22) = FJ(J,22) + FOCEAN(I,J)*RSI(I,J)*DEW(2)*ZATMO(I,J) 4507. FJ(J,23) = FJ(J,23) + FOCEAN(I,J)*RSI(I,J)*TCS 4507.1 FJ(J,26) = FJ(J,26) + FOCEAN(I,J)*RSI(I,J)*QA1UP 4507.2 FJ(J,27) = FJ(J,27) + FOCEAN(I,J)*RSI(I,J)*QA1DN 4507.3 FJ(J,28) = FJ(J,28) + FOCEAN(I,J)*RSI(I,J)*QS 4507.4 FJ(J,29) = FJ(J,29) + FOCEAN(I,J)*RSI(I,J)*QG 4507.5 FJ(J,35) = FJ(J,35) + FOCEAN(I,J)*RSI(I,J)*CIA 4507.6 FJ(J,43) = FJ(J,43) + FOCEAN(I,J)*RSI(I,J)*TA1UP 4507.7 FJ(J,44) = FJ(J,44) + FOCEAN(I,J)*RSI(I,J)*TA1DN 4507.8 FJ(J,48) = FJ(J,48) + FOCEAN(I,J)*RSI(I,J)*WS 4507.9 FJ(J,64) = FJ(J,64) + FOCEAN(I,J)*RSI(I,J)*BETAS 4507.91 FJ(J,65) = FJ(J,65) + FOCEAN(I,J)*RSI(I,J)*RIGS 4507.92 FJ(J,66) = FJ(J,66) + FOCEAN(I,J)*RSI(I,J)*CDM 4507.95 endif 4508. GO TO 240 4509. C**** 4510. C**** Glacial Ice 4511. C**** 4511.5 540 BJ(J, 6) = BJ(J, 6) + PTYPE*SRHDT 4512. BJ(J, 9) = BJ(J, 9) + PTYPE*TRHDT 4514. BJ(J,13) = BJ(J,13) + PTYPE*SNHDT 4515. BJ(J,14) = BJ(J,14) + PTYPE*EVHDT 4516. C BJ(J,19) = BJ(J,19) + PTYPE*DEW(3) 4516.1 C BJ(J,22) = BJ(J,22) + PTYPE*DEW(3)*ZATMO(I,J) 4517. BJ(J,23) = BJ(J,23) + PTYPE*TCS 4517.1 BJ(J,26) = BJ(J,26) + PTYPE*QA1UP 4517.2 BJ(J,27) = BJ(J,27) + PTYPE*QA1DN 4517.3 BJ(J,28) = BJ(J,28) + PTYPE*QS 4517.4 BJ(J,29) = BJ(J,29) + PTYPE*QG 4517.5 BJ(J,35) = BJ(J,35) + PTYPE*CIA 4517.6 BJ(J,43) = BJ(J,43) + PTYPE*TA1UP 4517.7 BJ(J,44) = BJ(J,44) + PTYPE*TA1DN 4517.8 BJ(J,48) = BJ(J,48) + PTYPE*WS 4517.9 BJ(J,64) = BJ(J,64) + PTYPE*BETAS 4517.91 BJ(J,65) = BJ(J,65) + PTYPE*RIGS 4517.92 BJ(J,66) = BJ(J,66) + PTYPE*CDM 4518. GO TO 260 4519. C**** 4520. C**** Ground 4521. C**** 4522. 560 MO (I,J,1) = MO (I,J,1) + PTYPE*(ARUNS+ARUNU)*DXYP(J) 4523. G0M(I,J,1) = G0M(I,J,1) + PTYPE* AERUNU *DXYP(J) 4523.5 AJ(J, 6) = AJ(J, 6) + PTYPE*SRHDT 4524. AJ(J, 9) = AJ(J, 9) + PTYPE*TRHDT 4526. AJ(J,13) = AJ(J,13) + PTYPE*SNHDT 4527. AJ(J,14) = AJ(J,14) + PTYPE*EVHDT 4528. AJ(J,19) = AJ(J,19) + PTYPE*DEW(4) 4529. AJ(J,23) = AJ(J,23) + PTYPE*TCS 4529.1 AJ(J,26) = AJ(J,26) + PTYPE*QA1UP 4529.2 AJ(J,27) = AJ(J,27) + PTYPE*QA1DN 4529.3 AJ(J,28) = AJ(J,28) + PTYPE*QS 4529.4 AJ(J,29) = AJ(J,29) + PTYPE*QG 4529.5 AJ(J,35) = AJ(J,35) + PTYPE*CIA 4529.6 AJ(J,43) = AJ(J,43) + PTYPE*TA1UP 4529.7 AJ(J,44) = AJ(J,44) + PTYPE*TA1DN 4529.8 AJ(J,48) = AJ(J,48) + PTYPE*WS 4529.9 AJ(J,64) = AJ(J,64) + PTYPE*BETAS 4529.91 AJ(J,65) = AJ(J,65) + PTYPE*RIGS 4529.92 AJ(J,66) = AJ(J,66) + PTYPE*CDM 4530. AJ(J,34) = AJ(J,34) - PTYPE*ARUNS 4531. AJ(J,38) = AJ(J,38) - PTYPE*ARUNU 4532. AJ(J,42) = AJ(J,42) - PTYPE*AERUNS 4533. AJ(J,47) = AJ(J,47) - PTYPE*(AERUNU-AERUNS) 4533.5 AJ(J,22) = AJ(J,22) + PTYPE*(DEW(4)-ARUNS-ARUNU)*ZATMO(I,J) 4533.6 CJ(J,22) = CJ(J,22) + PTYPE*(ARUNS+ARUNU)*ZATMO(I,J) 4533.7 AIJ(I,J,32) = AIJ(I,J,32) + (AEPB+AEPC) 4533.8 AIJ(I,J,33) = AIJ(I,J,33) + AEPP 4534. AIJ(I,J,44) = AIJ(I,J,44) + ARUNS 4535. AIJ(I,J,45) = AIJ(I,J,45) + ARUNU 4536. AIJ(I,J,56) = AIJ(I,J,56) - DEW(ITYPE) 4537. AIJ(I,J,60) = AIJ(I,J,60) + (SRHDT+TRHDT+SNHDT+EVHDT) 4538. AIJ(I,J, 7) = AIJ(I,J, 7) + WEARTH 4539. BLDATA(I,J,1) = TCG(4) 4540. BLDATA(I,J,2) = WEARTH 4541. BLDATA(I,J,3) = WS 4542. C**** 4543. C**** Test for dry convection and calculate vertically mixed variables 4544. C**** 4545. 600 MABY2 = .5*MA(I,J,1) 4545.1 QMUPO = .5*(Q0M(I,J,1)+.5*QZM(I,J,1)) 4546. QMDNSO = .5*(Q0M(I,J,1)-.5*QZM(I,J,1)) 4547. ENUPO = .5*(H0M(I,J,1)+.5*HZM(I,J,1))*PKUP(I,J,1) 4548. ENDNSO = .5*(H0M(I,J,1)-.5*HZM(I,J,1))*PKDN(I,J,1) 4548.1 C MDNSI = MABY2 - DEWS 4548.2 C MAN = MA(I,J,1) - DEWS 4548.3 C MSUMN = MSUM(I,J) - DEWS 4548.4 MDNSI = MABY2 4548.5 MAN = MA(I,J,1) 4548.6 MSUMN = MSUM(I,J) 4548.7 PKDNN = (GRAV*(MSUMN-.25*MAN))**RKAP 4548.8 PKUPN = (GRAV*(MSUMN-.75*MAN))**RKAP 4549. C**** Copy H0M,Q0M to H0MT,Q0MT which will be fixed while H0M changes 4550. DO 610 L=1,LBLM 4551. H0MT(L) = H0M(I,J,L) 4551.1 HXMT(L) = HXM(I,J,L) 4551.2 HYMT(L) = HYM(I,J,L) 4552. HZMT(L) = HZM(I,J,L) 4553. Q0MT(L) = Q0M(I,J,L) 4553.1 QXMT(L) = QXM(I,J,L) 4553.2 QYMT(L) = QYM(I,J,L) 4554. 610 QZMT(L) = QZM(I,J,L) 4555. C**** Loop over surface types 4556. DO 700 ITYPE=1,4 4557. IF(FTYPE(ITYPE).le.0.) GO TO 700 4558. PTYPE = FTYPE(ITYPE) 4559. C MDNO = MABY2 - DEWS + DEW(ITYPE) 4559.1 MDNO = MABY2 4560. QMDNI = QMDNSO*MDNO/MABY2 - DEW(ITYPE)*DXYP(J) 4561. ENDNI = ENDNSO*MDNO/MABY2 - 4562. C * (HDT(ITYPE) + SHCV*TKG*DEW(ITYPE))*DXYP(J) 4562.1 * (HDT(ITYPE) )*DXYP(J) 4564. C**** Test whether first layer is unstable with respect to second layer 4564.4 MK22 = MA(I,J,2)*(PKDN(I,J,2)+PKUP(I,J,2)) 4564.5 EN22 = (H0MT(2)-.5*HZMT(2))*PKDN(I,J,2) + 4564.6 * (H0MT(2)+.5*HZMT(2))*PKUP(I,J,2) 4565. IF(2.*(ENDNI+ENUPO)*MK22 .gt. EN22*MAN*(PKDNN+PKUPN)) GO TO 620 4566. C**** First layer is stable with respect to second layer 4567. ENDNN = ENDNI - (ENDNI-ENUPO)*.5*DEWS/MA(I,J,1) 4567.1 ENUPN = ENUPO + (ENDNI-ENUPO)*.5*DEWS/MA(I,J,1) 4567.2 H0MN = ENDNN/PKDNN + ENUPN/PKUPN 4567.3 HZMN = (ENUPN/PKUPN - ENDNN/PKDNN)*2. 4567.5 Q0MN = QMDNI + QMUPO 4567.6 QZMN = (QMUPO - QMDNI)*(1.-DEWS/MA(I,J,L)) 4568. H0M(I,J,1) = H0M(I,J,1) + PTYPE*(H0MN-H0MT(1)) 4569. HZM(I,J,1) = HZM(I,J,1) + PTYPE*(HZMN-HZMT(1)) 4570. Q0M(I,J,1) = Q0M(I,J,1) + PTYPE*(Q0MN-Q0MT(1)) 4571. QZM(I,J,1) = QZM(I,J,1) + PTYPE*(QZMN-QZMT(1)) 4572. MAS = MAN 4572.5 LMAX=1 4573. GO TO 660 4574. C**** Mix air through the first two unstable layers 4575. 620 MAS = MAN + MA(I,J,2) 4576. MKS2 = MAN*(PKDNN+PKUPN) + MK22 4577. ENS2 = 2.*(ENDNI+ENUPO) + EN22 4578. EXS2 = HXMT(1)*(PKDNN+PKUPN) + HXMT(2)*(PKDN(I,J,2)+PKUP(I,J,2)) 4579. EYS2 = HYMT(1)*(PKDNN+PKUPN) + HYMT(2)*(PKDN(I,J,2)+PKUP(I,J,2)) 4580. QMS = (QMDNI+QMUPO) + Q0MT(2) 4581. QXS = QXMT(1) + QXMT(2) 4582. QYS = QYMT(1) + QYMT(2) 4583. C**** Mix air through subsequent unstable layers 4584. L=3 4584.5 630 MKL2 = MA(I,J,L)*(PKDN(I,J,L)+PKUP(I,J,L)) 4585. ENL2 = (H0MT(L)-.5*HZMT(L))*PKDN(I,J,L) + 4585.1 * (H0MT(L)+.5*HZMT(L))*PKUP(I,J,L) 4586. IF(ENS2*MKL2 .le. ENL2*MKS2) GO TO 640 4587. MAS = MAS + MA(I,J,L) 4588. MKS2 = MKS2 + MKL2 4589. ENS2 = ENS2 + ENL2 4590. EXS2 = EXS2 + HXMT(L)*(PKDN(I,J,L)+PKUP(I,J,L)) 4591. EYS2 = EYS2 + HYMT(L)*(PKDN(I,J,L)+PKUP(I,J,L)) 4592. QMS = QMS + Q0MT(L) 4593. QXS = QXS + QXMT(L) 4594. QYS = QYS + QYMT(L) 4595. L=L+1 4596. IF(L.le.LBLM) GO TO 630 4597. C**** 4598. C**** Store thoroughly mixed air into prognostic variables 4599. C**** 4600. 640 LMAX=L-1 4601. H0M(I,J,1) = H0M(I,J,1) + PTYPE*(MAN*ENS2/MKS2-H0MT(1)) 4602. HXM(I,J,1) = HXM(I,J,1) + PTYPE*(MAN*EXS2/MKS2-HXMT(1)) 4603. HYM(I,J,1) = HYM(I,J,1) + PTYPE*(MAN*EYS2/MKS2-HYMT(1)) 4604. HZM(I,J,1) = HZM(I,J,1) - PTYPE*HZMT(1) 4605. Q0M(I,J,1) = Q0M(I,J,1) + PTYPE*(MAN*QMS/MAS-Q0MT(1)) 4606. QXM(I,J,1) = QXM(I,J,1) + PTYPE*(MAN*QXS/MAS-QXMT(1)) 4607. QYM(I,J,1) = QYM(I,J,1) + PTYPE*(MAN*QYS/MAS-QYMT(1)) 4608. QZM(I,J,1) = QZM(I,J,1) - PTYPE*QZMT(1) 4610. DO 650 L=2,LMAX 4611. H0M(I,J,L) = H0M(I,J,L) + PTYPE*(MA(I,J,L)*ENS2/MKS2-H0MT(L)) 4612. HXM(I,J,L) = HXM(I,J,L) + PTYPE*(MA(I,J,L)*EXS2/MKS2-HXMT(L)) 4613. HYM(I,J,L) = HYM(I,J,L) + PTYPE*(MA(I,J,L)*EYS2/MKS2-HYMT(L)) 4614. HZM(I,J,L) = HZM(I,J,L) - PTYPE*HZMT(L) 4615. Q0M(I,J,L) = Q0M(I,J,L) + PTYPE*(MA(I,J,L)*QMS/MAS-Q0MT(L)) 4616. QXM(I,J,L) = QXM(I,J,L) + PTYPE*(MA(I,J,L)*QXS/MAS-QXMT(L)) 4617. QYM(I,J,L) = QYM(I,J,L) + PTYPE*(MA(I,J,L)*QYS/MAS-QYMT(L)) 4618. 650 QZM(I,J,L) = QZM(I,J,L) - PTYPE*QZMT(L) 4620. C**** 4621. C**** Apply surface friction and mix momentum through boundary layer 4622. C**** 4623. 660 IF(QPOLE) GO TO 670 4624. UVMS(1) = UAT(ID(1),1,1)*MA(I,J,1) - DMU(ITYPE) 4625. UVMS(2) = UAT(ID(2),1,1)*MA(I,J,1) - DMU(ITYPE) 4626. UVMS(3) = UAT(ID(3),1,1)*MA(I,J,1) - DMV(ITYPE) 4627. UVMS(4) = UAT(ID(4),1,1)*MA(I,J,1) - DMV(ITYPE) 4628. GO TO 680 4629. 670 IF(J.eq.1) THEN 4630. DO 671 II=1,IM 4631. 671 UVMS(II) = VAT(II,JVPO,1)*MA(I,J,1) - 4632. * (DMV(ITYPE)*COSI(II)+DMU(ITYPE)*SINI(II)) 4633. ELSE 4634. DO 672 II=1,IM 4635. 672 UVMS(II) = VAT(II,JVPO,1)*MA(I,J,1) - 4636. * (DMV(ITYPE)*COSI(II)-DMU(ITYPE)*SINI(II)) 4637. ENDIF 4638. DVDOTV = DMU(ITYPE)*U1 + DMV(ITYPE)*V1 4639. IF(DVDOTV.le.0) THEN 4640. KMAX=IM 4641. ELSE 4642. KMAX=IM+1 4643. DMPOLE = DVDOTV/(U1*U1+V1*V1) 4644. IF(DMPOLE.gt..5*MA(I,J,1)) DMPOLE= .5*MA(I,J,1) 4645. UVMS(KMAX) = UAT(1,J,1)*(MA(I,J,1) - DMPOLE) 4646. DMUA(2,J,ITYPE) = DMUA(2,J,ITYPE) + DMPOLE 4647. ENDIF 4648. 680 DO 682 K=1,KMAX 4649. DO 681 L=2,LMAX 4650. 681 UVMS(K) = UVMS(K) + UAT(ID(K),1,L)*MA(I,J,L) 4651. UVMS(K) = UVMS(K)/MAS 4652. DO 682 L=1,LMAX 4653. 682 UA(ID(K),1,L) = UA(ID(K),1,L)+PTYPE*(UVMS(K)-UAT(ID(K),1,L))*RA(K) 4654. C**** Accumulate boundary layer diagnostics due to dry convection 4655. IF(QPOLE) THEN 4656. IF(KMAX.le.IM) GO TO 693 4657. DO 691 L=1,LMAX 4658. 691 AJL(J,L,30) = AJL(J,L,30) + 4659. + PTYPE*(UVMS(IM+1)-UAT(1,J,L))*MA(I,J,L)*2. 4660. ELSE 4661. DO 692 L=1,LMAX 4662. 692 AJL(J,L,30) = AJL(J,L,30) + 4663. + PTYPE*(UVMS(1)+UVMS(2)-UAT(IM1,J,L)-UAT(I,J,L))*MA(I,J,L) 4664. ENDIF 4665. 693 IF(.not.QDIAGD .or. LMAX.le.1) GO TO 696 4666. DO 694 KR=1,4 4667. IF(I.eq.IJDD(1,KR) .and. J.eq.IJDD(2,KR)) GO TO 695 4668. 694 CONTINUE 4669. GO TO 696 4670. 695 ADAILY(JHOUR,47,KR) = ADAILY(JHOUR,47,KR) + PTYPE 4671. ADAILY(JHOUR,48,KR) = ADAILY(JHOUR,48,KR) + PTYPE*LMAX 4672. 696 CONTINUE 4673. 700 CONTINUE 4674. C MA(I,J,1) = MAN 4675. C PKDN(I,J,1) = PKDNN 4676. C PKUP(I,J,1) = PKUPN 4676.5 C MSUM(I,J) = MSUMN 4676.6 UAS(I,J) = USS 4676.7 VAS(I,J) = VSS 4677. C**** 4678. C**** Accumulate diagnostics 4679. C**** 4680. C**** Quantities accumulated for latitude-longitude maps in DIAGIJ 4681. AIJ(I,J, 4) = AIJ(I,J, 4) + SNHDTS 4682. AIJ(I,J, 6) = AIJ(I,J, 6) - DEWS 4683. AIJ(I,J,22) = AIJ(I,J,22) + TRHDTS 4684. AIJ(I,J,23) = AIJ(I,J,23) + (SRHDTS+TRHDTS+SNHDTS+EVHDTS) 4684.5 AIJ(I,J,35) = AIJ(I,J,35) + TCSS 4685. IF(QRAD) AIJ(I,J,21) = AIJ(I,J,21) + TRHDTS/DTS 4686. IF(TCSS.lt.TMNMX(I,J,1)) TMNMX(I,J,1) = TCSS 4687. IF(TCSS.gt.TMNMX(I,J,2)) TMNMX(I,J,2) = TCSS 4688. IF(TCGS.lt.TMNMX(I,J,3)) TMNMX(I,J,3) = TCGS 4689. IF(TCGS.gt.TMNMX(I,J,4)) TMNMX(I,J,4) = TCGS 4690. TMNMX(I,J,6) = TMNMX(I,J,6) + TCSS 4691. IF(.not.QDASRF) GO TO 820 4692. AIJ(I,J,28) = AIJ(I,J,28) + TCGS 4694. AIJ(I,J,36) = AIJ(I,J,36) + USS 4695. AIJ(I,J,37) = AIJ(I,J,37) + VSS 4696. AIJ(I,J,47) = AIJ(I,J,47) + CIAS 4697. AIJ(I,J,48) = AIJ(I,J,48) + DMUS 4698. AIJ(I,J,49) = AIJ(I,J,49) + DMVS 4699. AIJ(I,J,50) = AIJ(I,J,50) + QSS 4700. AIJ(I,J,51) = AIJ(I,J,51) + WS 4700.1 AIJ(I,J,76) = AIJ(I,J,76) + TKES 4701. C**** Quantities accumulated hourly for DIAGD 4702. 820 IF(.NOT.QDIAGD) GO TO 900 4703. DO 830 KR=1,4 4704. IF(I.eq.IJDD(1,KR) .and. J.eq.IJDD(2,KR)) GO TO 840 4705. 830 CONTINUE 4706. GO TO 900 4707. 840 CONTINUE 4708. ADAILY(JHOUR, 6,KR) = ADAILY(JHOUR, 6,KR) + PS 4709. ADAILY(JHOUR, 7,KR) = ADAILY(JHOUR, 7,KR) + H0M(I,J,5)*PGK 4710. ADAILY(JHOUR, 8,KR) = ADAILY(JHOUR, 8,KR) + H0M(I,J,4)*PGK 4711. ADAILY(JHOUR, 9,KR) = ADAILY(JHOUR, 9,KR) + H0M(I,J,3)*PGK 4712. ADAILY(JHOUR,10,KR) = ADAILY(JHOUR,10,KR) + H0M(I,J,2)*PGK 4713. ADAILY(JHOUR,11,KR) = ADAILY(JHOUR,11,KR) + H0M(I,J,1)*PGK 4714. ADAILY(JHOUR,12,KR) = ADAILY(JHOUR,12,KR) + TCSS 4715. ADAILY(JHOUR,13,KR) = ADAILY(JHOUR,13,KR) + TCGS 4716. ADAILY(JHOUR,14,KR) = ADAILY(JHOUR,14,KR) + Q0M(I,J,5) 4717. ADAILY(JHOUR,15,KR) = ADAILY(JHOUR,15,KR) + Q0M(I,J,4) 4718. ADAILY(JHOUR,16,KR) = ADAILY(JHOUR,16,KR) + Q0M(I,J,3) 4719. ADAILY(JHOUR,17,KR) = ADAILY(JHOUR,17,KR) + Q0M(I,J,2) 4720. ADAILY(JHOUR,18,KR) = ADAILY(JHOUR,18,KR) + Q0M(I,J,1) 4721. ADAILY(JHOUR,19,KR) = ADAILY(JHOUR,19,KR) + QSS 4722. ADAILY(JHOUR,20,KR) = ADAILY(JHOUR,20,KR) + QGS 4723. ADAILY(JHOUR,28,KR) = ADAILY(JHOUR,28,KR) + SRHDTS 4724. ADAILY(JHOUR,29,KR) = ADAILY(JHOUR,29,KR) + TRHDTS 4725. ADAILY(JHOUR,30,KR) = ADAILY(JHOUR,30,KR) + SNHDTS 4726. ADAILY(JHOUR,31,KR) = ADAILY(JHOUR,31,KR) + EVHDTS 4727. ADAILY(JHOUR,32,KR) = ADAILY(JHOUR,32,KR) + 4728. * (SRHDTS+TRHDTS+SNHDTS+EVHDTS) 4729. ADAILY(JHOUR,33,KR) = ADAILY(JHOUR,33,KR) + U1 4730. ADAILY(JHOUR,34,KR) = ADAILY(JHOUR,34,KR) + V1 4731. ADAILY(JHOUR,36,KR) = ADAILY(JHOUR,36,KR) + USS 4732. ADAILY(JHOUR,37,KR) = ADAILY(JHOUR,37,KR) + VSS 4733. ADAILY(JHOUR,38,KR) = ADAILY(JHOUR,38,KR) + WS 4734. ADAILY(JHOUR,39,KR) = ADAILY(JHOUR,39,KR) + CIAS 4735. ADAILY(JHOUR,41,KR) = ADAILY(JHOUR,41,KR) + RIGSS 4736. ADAILY(JHOUR,42,KR) = ADAILY(JHOUR,42,KR) + CDMS 4737. C ADAILY(JHOUR,43,KR) = ADAILY(JHOUR,43,KR) + CDHS 4738. ADAILY(JHOUR,50,KR) = ADAILY(JHOUR,50,KR) - DEWS 4739. C**** End of outside I loop 4740. 900 IM1=I 4741. C**** End of outside J loop 4742. 910 CONTINUE 4743. C**** End of outside loop over surface time steps 4744. 920 CONTINUE 4745. RETURN 4746. END 7000. 7001. SUBROUTINE DRYCNV 7002. C**** 7003. C**** DRYCNV causes a thorough overtunring of the air throughout 7004. C**** the layers which were originally thermally unstable. Dry 7005. C**** convection originating from the first layer is calculated 7006. C**** in subroutine SURFCE. Dry convection originating above the 7007. C**** first layer is calculated in DRYCNV. 7008. C**** 7009. INCLUDE 'C070.COM' 7010. LOGICAL*4 QPOLE 7011. COMMON /WORK00/ PKDN(IM,JM,LMA),PKUP(IM,JM,LMA) 7012. COMMON /WORK01/ UT(IM,JM,LMA),VT(IM,JM,LMA), 7013. * UMS(IM+1),RA(IM+1),ID(IM+1) 7014. C**** Determine highest layer from which convection may start 7015. LMINM=LMA-1 7017. C**** Load U,V into UT,VT. UT,VT will be fixed during dry 7018. C**** convection while U,V will be updated. 7019. DO 10 L=2,LMA 7020. DO 10 J=1,JM 7021. DO 10 I=1,IM 7022. UT(I,J,L) = UA(I,J,L) 7023. 10 VT(I,J,L) = VA(I,J,L) 7024. C**** 7025. C**** Main J loop 7026. C**** 7027. DO 400 J=1,JM 7028. QPOLE = J.eq.1 .OR. J.eq.JM 7029. IF(QPOLE) GO TO 120 7030. C**** Conditions at non-polar points 7031. IMAX=IM 7032. KMAX=4 7033. RA(1) = .5 7034. RA(2) = .5 7035. RA(3) = RAMVS(J) 7036. RA(4) = RAMVN(J) 7037. GO TO 200 7038. C**** Conditions at the poles 7039. 120 IMAX=1 7040. KMAX=IM+1 7041. JVPO=1 7042. IF(J.eq.JM) JVPO=JM-1 7043. DO 140 I=1,IM 7044. RA(I) = RAMVN(1) 7045. 140 ID(I) = I + (JVPO-1)*IM + IM*JM*LMA 7046. RA(IM+1) = 1. 7047. ID(IM+1) = 1 + (J-1)*IM 7048. C**** 7049. C**** Main I loop 7050. C**** 7051. 200 IM1=IM 7052. DO 400 I=1,IMAX 7053. C**** 7054. C**** Loop over dry convection starting layers, test for convection, 7055. C**** and calculate vertically mixed variables 7056. C**** 7057. LMIN=2 7057.2 210 MKS2 = MA(I,J,LMIN)* (PKDN(I,J,LMIN) +PKUP(I,J,LMIN)) 7057.3 ENS2 = (H0M(I,J,LMIN)-.5*HZM(I,J,LMIN))*PKDN(I,J,LMIN) + 7057.4 + (H0M(I,J,LMIN)+.5*HZM(I,J,LMIN))*PKUP(I,J,LMIN) 7057.5 MKL2 = MA(I,J,LMIN+1)* (PKDN(I,J,LMIN+1) +PKUP(I,J,LMIN+1)) 7057.6 ENL2 = (H0M(I,J,LMIN+1)-.5*HZM(I,J,LMIN+1))*PKDN(I,J,LMIN+1) + 7057.7 + (H0M(I,J,LMIN+1)+.5*HZM(I,J,LMIN+1))*PKUP(I,J,LMIN+1) 7058. IF(ENS2*MKL2 .gt. ENL2*MKS2) GO TO 220 7060. IF(LMIN.ge.LMINM) GO TO 400 7061. LMIN=LMIN+1 7062. GO TO 210 7063. C**** Mix air through the first two unstable layers 7063.5 220 MAS = MA(I,J,LMIN) + MA(I,J,LMIN+1) 7064. MKS2 = MKS2 + MKL2 7065. ENS2 = ENS2 + ENL2 7066. EXS2 = HXM(I,J,LMIN )*(PKDN(I,J,LMIN )+PKUP(I,J,LMIN )) + 7066.1 * HXM(I,J,LMIN+1)*(PKDN(I,J,LMIN+1)+PKUP(I,J,LMIN+1)) 7067. EYS2 = HYM(I,J,LMIN )*(PKDN(I,J,LMIN )+PKUP(I,J,LMIN )) + 7067.1 * HYM(I,J,LMIN+1)*(PKDN(I,J,LMIN+1)+PKUP(I,J,LMIN+1)) 7068. QMS = Q0M(I,J,LMIN) + Q0M(I,J,LMIN+1) 7069. QXS = QXM(I,J,LMIN) + QXM(I,J,LMIN+1) 7070. QYS = QYM(I,J,LMIN) + QYM(I,J,LMIN+1) 7071. C**** Mix air through subsequent unstable layers 7072. DO 230 L=LMIN+2,LMA 7072.5 MKL2 = MA(I,J,L)* (PKDN(I,J,L)+ PKUP(I,J,L)) 7073. ENL2 = (H0M(I,J,L)-.5*HZM(I,J,L))*PKDN(I,J,L) + 7073.1 + (H0M(I,J,L)+.5*HZM(I,J,L))*PKUP(I,J,L) 7074. IF(ENS2*MKL2 .le. ENL2*MKS2) GO TO 300 7074.5 MAS = MAS + MA(I,J,L) 7075. MKS2 = MKS2 + MKL2 7076. ENS2 = ENS2 + ENL2 7077. EXS2 = EXS2 + HXM(I,J,L)*(PKDN(I,J,L)+PKUP(I,J,L)) 7078. EYS2 = EYS2 + HYM(I,J,L)*(PKDN(I,J,L)+PKUP(I,J,L)) 7079. QMS = QMS + Q0M(I,J,L) 7080. QXS = QXS + QXM(I,J,L) 7081. 230 QYS = QYS + QYM(I,J,L) 7083. C L=LMA+1 7083.1 IF(L.ne.LMA+1) STOP 'DRYCNV 7083' 7084. C**** 7085. C**** Store thoroughly mixed air into prognostic variables 7086. C**** 7087. 300 LMAX=L-1 7089. DO 310 L=LMIN,LMAX 7091. H0M(I,J,L) = MA(I,J,L)*ENS2/MKS2 7092. HXM(I,J,L) = MA(I,J,L)*EXS2/MKS2 7093. HYM(I,J,L) = MA(I,J,L)*EYS2/MKS2 7094. HZM(I,J,L) = 0. 7095. Q0M(I,J,L) = MA(I,J,L)*QMS/MAS 7096. QXM(I,J,L) = MA(I,J,L)*QXS/MAS 7097. QYM(I,J,L) = MA(I,J,L)*QYS/MAS 7098. 310 QZM(I,J,L) = 0. 7099. C**** Mix momentum 7100. IF(QPOLE) GO TO 320 7101. ID(1) = IM1+(J-1)*IM 7102. ID(2) = I +(J-1)*IM 7103. ID(4) = ID(2)+IM*JM*LMA 7104. ID(3) = ID(4)-IM 7105. 320 DO 340 K=1,KMAX 7106. UMS(K) = 0. 7107. DO 330 L=LMIN,LMAX 7108. 330 UMS(K) = UMS(K) + UT(ID(K),1,L)*MA(I,J,L) 7109. UMS(K) = UMS(K)/MAS 7110. DO 340 L=LMIN,LMAX 7111. 340 UA(ID(K),1,L) = UA(ID(K),1,L) + (UMS(K)-UT(ID(K),1,L))*RA(K) 7112. IF(QPOLE) GO TO 360 7113. DO 350 L=LMIN,LMAX 7114. 350 AJL(J,L,30) = AJL(J,L,30) + 7115. * (UMS(1)+UMS(2)-UT(IM1,J,L)-UT(I,J,L))*MA(I,J,L) 7116. GO TO 390 7117. 360 DO 370 L=LMIN,LMAX 7118. 370 AJL(J,L,30) = AJL(J,L,30) + 2.*(UMS(IM+1)-UT(1,J,L))*MA(1,J,L) 7119. C**** End of loop over dry convection starting layer 7120. 390 LMIN=LMAX+1 7121. IF(LMIN.le.LMINM) GO TO 210 7122. 400 IM1=I 7123. C**** 7124. C**** Dry convection within a layer, force HZ to be non-negative 7125. C**** 7126. DO 510 L=1,LMA 7127. DO 510 J=1,JM 7128. IMAX=IM 7129. IF(J.eq.1 .or. J.eq.JM) IMAX=1 7130. DO 510 I=1,IMAX 7131. IF(HZM(I,J,L).ge.0.) GO TO 510 7131.5 H0M(I,J,L) = ((H0M(I,J,L)-.5*HZM(I,J,L))*PKDN(I,J,L) + 7131.6 * (H0M(I,J,L)+.5*HZM(I,J,L))*PKUP(I,J,L)) / 7131.7 * (PKDN(I,J,L) + PKUP(I,J,L)) 7132. HZM(I,J,L) = 0. 7133. QZM(I,J,L) = 0. 7134. 510 CONTINUE 7135. RETURN 7136. END 8000. 8001. SUBROUTINE ASDRAG 8002. C**** 8003. C**** SDRAG puts a drag on the model's top layer 8004. C**** 8005. INCLUDE 'C070.COM' 8006. PARAMETER (IQ1=IM/4+1, IQ2=IM/2+1, IQ3=3*IM/4+1) 8007. COMMON /WORK00/ PKDN(IM,JM,LMA),PKUP(IM,JM,LMA) 8008. COMMON /WORK01/ WLM(IM,JM),X(IM,JM),UX(IM,JM),VX(IM,JM) 8012. C**** Calculate wind magnitude for layer LMA 8013. WLM(1,1) = SQRT(.5*(VA( 1,1,LMA)**2+VA(IQ1,1,LMA)**2 8014. * +VA(IQ2,1,LMA)**2+VA(IQ3,1,LMA)**2)) 8015. WLM(1,JM) = SQRT(.5*(VA( 1,JM-1,LMA)**2+VA(IQ1,JM-1,LMA)**2 8016. * +VA(IQ2,JM-1,LMA)**2+VA(IQ3,JM-1,LMA)**2)) 8017. DO 10 J=2,JM-1 8018. IM1=IM 8019. DO 10 I=1,IM 8020. WLM(I,J) = SQRT(.5*(UA(IM1,J,LMA)**2+UA(I,J,LMA)**2 8021. * +VA(I,J-1,LMA)**2+VA(I,J,LMA)**2)) 8022. 10 IM1=I 8023. C**** Calculate the stratospheric drag factor X 8024. DO 20 J=1,JM 8025. IMAX=IM 8026. IF(J.eq.1 .OR. J.eq.JM) IMAX=1 8027. DO 20 I=1,IMAX 8028. CDN = .0005d0 + .00005d0*WLM(I,J) 8029. 20 X(I,J) = DTS*CDN*WLM(I,J)*GRAV*MSTRAT*DXYP(J)/ 8030. * (RKAP*H0M(I,J,LMA)*.5*(PKDN(I,J,LMA)+PKUP(I,J,LMA))) 8031. C**** Save current velocity components for layer LMA 8032. DO 30 J=1,JM 8033. DO 30 I=1,IM 8034. UX(I,J) = UA(I,J,LMA) 8035. 30 VX(I,J) = VA(I,J,LMA) 8036. C**** Reduce momentum components by stratospheric drag 8037. DO 40 I=1,IM 8038. VA(I,1,LMA) = VA(I,1,LMA) - X(1,1)*VX(I,1)*RAMVN(1) 8039. 40 VA(I,JM-1,LMA) = VA(I,JM-1,LMA) - X(1,JM)*VX(I,JM-1)*RAMVS(JM) 8040. DO 50 J=2,JM-1 8041. I=IM 8042. DO 50 IP1=1,IM 8043. UA(I,J,LMA) = UA(I,J,LMA) - (X(I,J)+X(IP1,J))*UX(I,J)*.5 8044. VA(I,J-1,LMA) = VA(I,J-1,LMA) - X(I,J)*VX(I,J-1)*RAMVS(J) 8045. VA(I,J,LMA) = VA(I,J,LMA) - X(I,J)*VX(I,J)*RAMVN(J) 8046. 50 I=IP1 8047. DO 60 I=1,IM 8048. UA(I, 1,LMA) = UA(I, 1,LMA) - X(1, 1)*UX(I, 1) 8049. 60 UA(I,JM,LMA) = UA(I,JM,LMA) - X(1,JM)*UX(I,JM) 8051. C**** Reduce the U-wind at the poles with 8-day time constant 8052. FAC = 1. - 1./(NDAY*8) 8053. DO 80 L=1,LMA 8054. DO 80 I=1,IM 8055. UA(I,1,L) = FAC*UA(I,1,L) 8056. 80 UA(I,JM,L) = FAC*UA(I,JM,L) 8057. RETURN 8058. END