1. C**** 2. C**** C070L.S Fortran source code for Lake routines 2003/12/02 3. C**** 4. C**** NASA/GISS Climate Model III 5. C**** 6. C**** NASA/Goddard Space Flight Center Programmed by: 7. C**** Institute for Space Studies Gary L. Russell 8. C**** 2880 Broadway, New York, NY 10025 Gavin A. Schmidt 9. C**** U.S.A. 10. C**** 11. SUBROUTINE PRECOL 12. C**** 13. C**** PRECLO adds the atmospheric precipitation to the open lake 14. C**** 15. INCLUDE 'C070.COM' 16. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 17. C**** 18. C**** OCENCB MO Liquid lake mass above sill depth (kg) 19. C**** G0M Liquid lake enthalpy (J) 20. C**** 21. C**** SEAICB RSI Horizontal ratio of lake ice to lake (1) 22. C**** 23. C**** FIXDCB FLAKE Lake fraction (1) 24. C**** 25. C**** PRECCB PREC Precipitation from atmosphere (kg/m**2) 26. C**** EPRE Energy of precipitation (J/m**2) 27. C**** 28. DO 10 J=2,JM-1 29. DO 10 I=1,IM 30. FOLAKE = FLAKE(I,J)*(1.-RSI(I,J)) 31. IF(FOLAKE.le.0.) GO TO 10 32. MO(I,J,1) = MO(I,J,1) + FOLAKE*DXYP(J)*PREC(I,J,1) 33. G0M(I,J,1) = G0M(I,J,1) + FOLAKE*DXYP(J)*EPRE(I,J,1) 34. 10 continue 37. RETURN 38. END 1000. 1001. SUBROUTINE OLAKE 1002. C**** 1003. C**** OLAKE receives the atmospheric fluxes of heat and water and 1004. C**** applies them to the open lake, updating the lake heat and mass 1005. C**** 1006. INCLUDE 'C070.COM' 1007. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 1008. * W0(IM,JM,4), E0(IM,JM,4) 1009. C**** 1010. C**** OCENCB MO Liquid lake mass above sill depth (kg) 1011. C**** G0M Liquid lake enthalpy (J) 1012. C**** 1013. C**** SEAICB RSI Horizontal ratio of lake ice to lake (1) 1014. C**** 1015. C**** FIXDCB FLAKE Lake fraction (1) 1016. C**** 1017. C**** FLUXCB W0(1) Dew received from atmosphere (kg/m**2) 1018. C**** E0(1) Energy received from atmosphere (J/m**2) 1019. C**** 1020. C**** Outside loop over J and I 1021. DO 90 J=2,JM 1022. IMAX=IM 1023. IF(J.eq.JM) IMAX=1 1024. DO 90 I=1,IMAX 1025. FOLAKE = FLAKE(I,J)*(1.-RSI(I,J)) 1026. IF(FOLAKE.le.0.) GO TO 90 1027. C**** Accumulate diagnostics 1028. CJ(J,18) = CJ(J,18) + FOLAKE*GZM(I,J,1) 1029. CJ(J,19) = CJ(J,19) + FOLAKE*W0(I,J,1) 1030. CJ(J,22) = CJ(J,22) + FOLAKE*W0(I,J,1)*ZATMO(I,J) 1030.1 CJ(J,50) = CJ(J,50) + MO(I,J,1)*BYDXYP(J) 1031. AIJ(I,J,53) = AIJ(I,J,53) - FOLAKE*W0(I,J,1) 1032. AIJ(I,J,57) = AIJ(I,J,57) + FOLAKE*E0(I,J,1) 1033. MO(I,J,1) = MO(I,J,1) + FOLAKE*W0(I,J,1)*DXYP(J) 1034. G0M(I,J,1) = G0M(I,J,1) + FOLAKE*E0(I,J,1)*DXYP(J) 1035. 90 continue 1036. RETURN 1037. END 2000. 2001. SUBROUTINE PRECLI 2002. C**** 2003. C**** PRECLI adds atmospheric precipitation to lake ice 2004. C**** 2005. INCLUDE 'C070.COM' 2006. PARAMETER (SNOMAX=91.66d0, dSNdRN=2., R2=1/XSI2,R3=1/XSI3) 2007. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 2008. C**** 2009. C**** OCENCB MO Liquid lake mass above sill depth (kg) 2010. C**** 2011. C**** SEAICE RSI Horizontal ratio of lake ice to lake (1) 2012. C**** MSI Mass of lake ice (kg/m^2) 2013. C**** HSI Enthalpy minus latent heat of lake ice (J/m^2) 2014. C**** 2015. C**** FIXDCB FLAKE Lake fraction (1) 2016. C**** 2017. C**** PRECCB PREC Precipitation from atmosphere (kg/m^2) 2018. C**** EPRE Energy of precipitation (J/m^2) 2019. C**** 2020. C**** Outside loop over J and I 2021. C**** 2022. DO 500 J=2,JM 2023. IMAX=IM 2024. IF(J.eq.JM) IMAX=1 2025. DO 500 I=1,IMAX 2026. FICE = FLAKE(I,J)*RSI(I,J) 2027. IF(PREC(I,J,1).le.0. .or. FICE.le.0.) GO TO 500 2028. C**** 2029. C**** Apply precipitation heat flux to HSI1 2030. C**** 2031. HSI(I,J,1) = HSI(I,J,1) + EPRE(I,J,1) 2032. IF(EPRE(I,J,1) .le. -PREC(I,J,1)*ELHM) GO TO 300 2033. IF(EPRE(I,J,1) .le. 0.) GO TO 200 2034. C**** 2035. C**** All precipitation is rain above 0øC 2036. C**** Rain compresses snow amount into ice 2037. C**** Calculate snow and ice melted or rain frozen in layer 1 2038. C**** 2039. RAIN = PREC(I,J,1) 2040. IF(HSI(I,J,1)/ELHM+XSI1*MSI(I,J,1) .le. 0.) GO TO 140 2041. C**** Warm rain cools to 0øC and melts some snow or ice 2042. MELT1 = HSI(I,J,1)/ELHM + XSI1*MSI(I,J,1) 2043. MO(I,J,1) = MO(I,J,1) + FICE*(RAIN+MELT1)*DXYP(J) 2044. CJ(J,21) = CJ(J,21) + FICE*(RAIN+MELT1)*ZATMO(I,J) 2045. DJ(J,21) = DJ(J,21) - FICE*(RAIN+MELT1)*ZATMO(I,J) 2046. C CJ(J,33) = CJ(J,33) + FICE*(RAIN+MELT1) ! in DIAGJ 2047. DJ(J,33) = DJ(J,33) - FICE*(RAIN+MELT1) 2048. IF(MSI(I,J,1)-ACE1I .le. MELT1) GO TO 120 2049. C**** Rain and melting compresses snow into ice 2050. FMSI2 = MIN (dSNdRN*(RAIN+MELT1), MSI(I,J,1)-ACE1I-MELT1) 2051. C FMSI1 = XSI1*CMPRS - XSI2*MELT1 2052. FHSI1 = - ELHM*(XSI1*FMSI2 - XSI2*MELT1) 2053. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 2054. MSI(I,J,1) = MSI(I,J,1) - MELT1 - FMSI2 2055. GO TO 400 2056. C**** All snow and some ice has melted 2057. C FMSI1 = XSI1*(MSI1-ACE1I) - MELT1 ! < 0 2058. 120 FMSI2 = MSI(I,J,1) - ACE1I - MELT1 ! < 0 2059. FHSI1 = HSI(I,J,2)*((XSI1/XSI2)*FMSI2-MELT1) / MSI(I,J,1) 2060. FHSI2 = HSI(I,J,3)*FMSI2*R3 / MSI(I,J,2) 2061. FHSI3 = HSI(I,J,4)*FMSI2 / MSI(I,J,2) 2062. MSI(I,J,1) = ACE1I 2063. GO TO 420 2064. C**** Rain compresses snow into ice, some rain will freeze 2065. 140 CMPRS = MIN (dSNdRN*RAIN, MSI(I,J,1)-ACE1I) 2066. IF(-HSI(I,J,1)/ELHM-XSI1*MSI(I,J,1) .lt. RAIN) GO TO 150 2067. C**** All rain freezes in layer 1 2068. C FREZ1 = RAIN 2069. C FMSI1 = XSI1*CMPRS + FREZ1 2070. FMSI2 = CMPRS + RAIN 2071. FHSI1 = HSI(I,J,1)*(XSI1*CMPRS+RAIN) / (XSI1*MSI(I,J,1)+RAIN) 2072. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 2073. MSI(I,J,1) = MSI(I,J,1) - CMPRS 2074. GO TO 400 2075. C**** Not all rain freezes in layer 1 2076. 150 FREZ1 = - HSI(I,J,1)/ELHM - XSI1*MSI(I,J,1) 2077. MO(I,J,1) = MO(I,J,1) + FICE*(RAIN-FREZ1)*DXYP(J) 2078. CJ(J,21) = CJ(J,21) + FICE*(RAIN-FREZ1)*ZATMO(I,J) 2079. DJ(J,21) = DJ(J,21) - FICE*(RAIN-FREZ1)*ZATMO(I,J) 2080. C CJ(J,33) = CJ(J,33) + FICE*(RAIN-FREZ1) ! in DIAGJ 2081. DJ(J,33) = DJ(J,33) - FICE*(RAIN-FREZ1) 2082. C FMSI1 = XSI1*CMPRS + FREZ1 2083. FMSI2 = CMPRS + FREZ1 2084. FHSI1 = - ELHM*(XSI1*CMPRS + FREZ1) 2085. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 2086. MSI(I,J,1) = MSI(I,J,1) - CMPRS 2087. GO TO 400 2088. C**** 2089. C**** Precipitation is a mixture of rain and snow at 0øC 2090. C**** Rain compresses snow into ice, some rain will freeze 2091. C**** Note that snow amount may exceed SNOMAX 2092. C**** 2093. 200 SNOW = -EPRE(I,J,1)/ELHM 2094. RAIN = PREC(I,J,1) - SNOW 2095. CMPRS = MIN (dSNdRN*RAIN, MSI(I,J,1)+SNOW-ACE1I) 2096. IF(-HSI(I,J,1)/ELHM-XSI1*MSI(I,J,1)-SNOW .lt. RAIN) GO TO 210 2097. C**** All rain freezes in layer 1 2098. C FREZ1 = RAIN 2099. FMSI1 = XSI1*CMPRS + XSI2*SNOW + RAIN 2100. FMSI2 = CMPRS + RAIN 2101. FHSI1 = HSI(I,J,1)*FMSI1 / (XSI1*MSI(I,J,1)+PREC(I,J,1)) 2102. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 2103. MSI(I,J,1) = MSI(I,J,1) + (SNOW-CMPRS) 2104. GO TO 400 2105. C**** Not all rain freezes in layer 1 2106. 210 FREZ1 = - HSI(I,J,1)/ELHM - XSI1*MSI(I,J,1) - SNOW 2107. MO(I,J,1) = MO(I,J,1) + FICE*(RAIN-FREZ1)*DXYP(J) 2108. CJ(J,21) = CJ(J,21) + FICE*(RAIN-FREZ1)*ZATMO(I,J) 2109. DJ(J,21) = DJ(J,21) - FICE*(RAIN-FREZ1)*ZATMO(I,J) 2110. C CJ(J,33) = CJ(J,33) + FICE*(RAIN-FREZ1) ! in DIAGJ 2111. DJ(J,33) = DJ(J,33) - FICE*(RAIN-FREZ1) 2112. C FMSI1 = XSI1*CMPRS + XSI2*SNOW + FREZ1 2113. FMSI2 = CMPRS + FREZ1 2114. FHSI1 = - ELHM*(XSI1*CMPRS + XSI2*SNOW + FREZ1) 2115. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 2116. MSI(I,J,1) = MSI(I,J,1) + (SNOW-CMPRS) 2117. GO TO 400 2118. C**** 2119. C**** All precipitation is snow which increases snow amount 2120. C**** 2121. 300 IF(MSI(I,J,1)+PREC(I,J,1) .gt. SNOMAX+ACE1I) GO TO 310 2122. C FMSI1 = XSI2*PREC ! > 0 2123. FHSI1 = HSI(I,J,1)*XSI2*PREC(I,J,1) /(XSI1*MSI(I,J,1)+PREC(I,J,1)) 2124. MSI(I,J,1) = MSI(I,J,1) + PREC(I,J,1) 2125. HSI(I,J,1) = HSI(I,J,1) - FHSI1 2126. HSI(I,J,2) = HSI(I,J,2) + FHSI1 2127. GO TO 500 2128. C**** Too much snow has accumulated, compress some into ice 2129. 310 FMSI2 = MSI(I,J,1) + PREC(I,J,1) - (.9d0*SNOMAX+ACE1I) ! > 0 2130. C FMSI1 = XSI2*PREC + XSI1*FMSI2 ! > 0 2131. FHSI1 = HSI(I,J,1)*(XSI2*PREC(I,J,1)+XSI1*FMSI2) / 2132. / (XSI1*MSI(I,J,1)+PREC(I,J,1)) 2133. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 2134. MSI(I,J,1) = .9d0*SNOMAX+ACE1I 2135. C**** 2136. C**** Advect ice (usually downward) 2137. C**** 2138. 400 FHSI3 = HSI(I,J,3)*FMSI2*(XSI4/XSI3) / MSI(I,J,2) 2139. 420 MSI(I,J,2) = MSI(I,J,2) + FMSI2 2140. HSI(I,J,1) = HSI(I,J,1) - FHSI1 2141. HSI(I,J,2) = HSI(I,J,2) + (FHSI1 - FHSI2) 2142. HSI(I,J,3) = HSI(I,J,3) + (FHSI2 - FHSI3) 2143. HSI(I,J,4) = HSI(I,J,4) + FHSI3 2144. 500 continue 2145. RETURN 2146. END 3000. 3001. SUBROUTINE LAKICE 3002. C**** 3003. C**** LAKICE receives the atmospheric fluxes of heat and water and 3004. C**** applies them to the lake ice, updating the heat, mass and snow 3005. C**** 3006. INCLUDE 'C070.COM' 3007. PARAMETER (RHOI=910., ALAMI=2.1762d0, ALPHA=1., TCFO=0., 3008. * dSNdML=2., R1=1/XSI1,R2=1/XSI2,R3=1/XSI3,R4=1/XSI4) 3009. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 3010. * W0(IM,JM,4), E0(IM,JM,4), E1(IM,JM,4) 3011. C**** 3012. C**** OCENCB MO Liquid lake mass above sill depth (kg) 3013. C**** 3014. C**** SEAICB RSI Horizontal ratio of lake ice to lake (1) 3015. C**** MSI Mass of sea ice (kg/m^2) 3016. C**** HSI Enthalpy minus latent heat of lake ice (J/m^2) 3017. C**** 3018. C**** FIXDCB FLAKE Lake fraction (1) 3019. C**** 3020. C**** FLUXCB W0(2) Dew received from atmosphere (kg/m^2) 3021. C**** E0(2) Energy received from atmosphere (J/m^2) 3022. C**** E1(2) Energy from layer 1 to layer 2 (J/m^2) 3023. C**** 3024. C**** Outside loop over J and I 3025. C**** 3026. DO 400 J=2,JM 3027. IMAX=IM 3028. IF(J.eq.JM) IMAX=1 3029. DO 400 I=1,IMAX 3030. FICE = FLAKE(I,J)*RSI(I,J) 3031. IF(FICE.le.0.) GO TO 400 3032. TSI1 = (HSI(I,J,1)*R1/MSI(I,J,1) + ELHM)/SHCI 3033. TSI2 = (HSI(I,J,2)*R2/MSI(I,J,1) + ELHM)/SHCI 3034. TSI3 = (HSI(I,J,3)*R3/MSI(I,J,2) + ELHM)/SHCI 3035. TSI4 = (HSI(I,J,4)*R4/MSI(I,J,2) + ELHM)/SHCI 3036. C**** Accumulate diagnostics 3037. DJ(J,30) = DJ(J,30) + FICE 3038. DJ(J,51) = DJ(J,51) + FICE*MSI(I,J,2) 3039. DJ(J,18) = DJ(J,18) + FICE*TSI1 3040. DJ(J,17) = DJ(J,17) + FICE*TSI2 3041. DJ(J,16) = DJ(J,16) + FICE*TSI3 3042. DJ(J,52) = DJ(J,52) + FICE*TSI4 3043. DJ(J,19) = DJ(J,19) + FICE*W0(I,J,2) 3044. DJ(J,22) = DJ(J,22) + FICE*W0(I,J,2)*ZATMO(I,J) 3045. AIJ(I,J, 1) = AIJ(I,J, 1) + FICE 3046. AIJ(I,J,54) = AIJ(I,J,54) - FICE*W0(I,J,2) 3047. AIJ(I,J,58) = AIJ(I,J,58) + FICE*E0(I,J,2) 3048. IF(MSI(I,J,1).le.ACE1I) GO TO 10 3049. DJ(J,31) = DJ(J,31) + FICE 3050. DJ(J,53) = DJ(J,53) + FICE*(MSI(I,J,1)-ACE1I) 3051. AIJ(I,J,2) = AIJ(I,J,2) + FICE 3052. AIJ(I,J,3) = AIJ(I,J,3) + FICE*(MSI(I,J,1)-ACE1I) 3053. C**** 3054. C**** Calculate and apply diffusive and surface energy fluxes 3055. C**** 3056. 10 HCI2 = SHCI*XSI2*MSI(I,J,1) 3057. HCI3 = SHCI*XSI3*MSI(I,J,2) 3058. HCI4 = SHCI*XSI4*MSI(I,J,2) 3059. dF2dTI = ALAMI*RHOI*DTS / (.5*XSI2*MSI(I,J,1)+.5*XSI3*MSI(I,J,2)) 3060. dF3dTI = ALAMI*RHOI*DTS*2. / MSI(I,J,2) 3061. dF4dTI = ALAMI*RHOI*DTS*2.*R4 / MSI(I,J,2) 3062. CEXP F2 = dF2dTI*(TSI2-TSI3) 3063. CEXP F3 = dF3dTI*(TSI3-TSI4) 3064. CEXP F4 = dF4dTI*(TSI4-TCFO) 3065. F2 = dF2dTI*(HCI2*(TSI2-TSI3) + ALPHA*E1(I,J,2)) / 3066. / (HCI2 + ALPHA*dF2dTI) 3067. F3 = dF3dTI*(HCI3*(TSI3-TSI4) + ALPHA*F2) / (HCI3 + ALPHA*dF3dTI) 3068. F4 = dF4dTI*(HCI4*(TSI4-TCFO) + ALPHA*F3) / (HCI4 + ALPHA*dF4dTI) 3069. HSI(I,J,1) = HSI(I,J,1) + (E0(I,J,2) - E1(I,J,2)) 3070. HSI(I,J,2) = HSI(I,J,2) + (E1(I,J,2) - F2) 3071. HSI(I,J,3) = HSI(I,J,3) + (F2 - F3) 3072. HSI(I,J,4) = HSI(I,J,4) + (F3 - F4) 3073. G0M(I,J,1) = G0M(I,J,1) + FICE*F4*DXYP(J) 3074. DJ(J,15) = DJ(J,15) - FICE*F4 3075. IF(HSI(I,J,1)/ELHM+XSI1*MSI(I,J,1)+W0(I,J,2) .le. 0.) GO TO 100 3076. C**** 3077. C**** Fluxes heat layer 1 to freezing point and melt some snow or ice 3078. C**** 3079. MELT1 = HSI(I,J,1)/ELHM + XSI1*MSI(I,J,1) + W0(I,J,2) 3080. MO(I,J,1) = MO(I,J,1) + FICE*MELT1*DXYP(J) 3081. CJ(J,22) = CJ(J,22) + FICE*MELT1*ZATMO(I,J) 3082. DJ(J,22) = DJ(J,22) - FICE*MELT1*ZATMO(I,J) 3083. C CJ(J,34) = CJ(J,34) + FICE*MELT1 ! in DIAGJ 3084. DJ(J,34) = DJ(J,34) - FICE*MELT1 3085. IF(W0(I,J,2) .gt. 0.) GO TO 30 3086. C**** Evaporation reduces snow or ice amount, MELT1 > 0 3087. IF(MSI(I,J,1)-ACE1I+W0(I,J,2) .gt. MELT1) GO TO 20 3088. C**** All snow and some ice evaporates and melts 3089. C**** Ice advection is upward into layer 2 from layer 3 3090. FMSI1 = XSI1*(MSI(I,J,1)-ACE1I) + W0(I,J,2) - MELT1 ! < 0 3091. FMSI2 = MSI(I,J,1)-ACE1I + W0(I,J,2) - MELT1 ! < 0 3092. GO TO 110 3093. C**** Evaporation and melting reduce snow amount 3094. C**** Ice advection is downward from layer 2 into layer 3 3095. 20 FMSI2 = MIN (dSNdML*MELT1, MSI(I,J,1)-ACE1I+W0(I,J,2)-MELT1) ! >0 3096. FMSI1 = XSI1*FMSI2 + XSI2*(W0(I,J,2) - MELT1) 3097. IF(FMSI1.lt.0.) FHSI1 = HSI(I,J,2)*FMSI1*R2 / MSI(I,J,1) 3098. IF(FMSI1.ge.0.) FHSI1 = - ELHM*FMSI1 3099. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 3100. MSI(I,J,1) = MSI(I,J,1) + W0(I,J,2) - MELT1 - FMSI2 3101. GO TO 310 3102. C**** Dew increases ice amount, MELT1 > 0 3103. 30 IF(MSI(I,J,1)-ACE1I .gt. MELT1) GO TO 40 3104. C**** All snow and some ice melts 3105. FMSI1 = XSI1*(MSI(I,J,1)-ACE1I) + W0(I,J,2) - MELT1 3106. FMSI2 = MSI(I,J,1) - ACE1I + W0(I,J,2) - MELT1 3107. IF(FMSI2.le.0.) GO TO 110 3108. C**** Ice advection is downward from layer 2 into layer 3 3109. IF(FMSI1.lt.0.) FHSI1 = HSI(I,J,2)*FMSI1*R2 / MSI(I,J,1) 3110. IF(FMSI1.ge.0.) FHSI1 = - ELHM*FMSI1 3111. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 3112. MSI(I,J,1) = ACE1I 3113. GO TO 310 3114. C**** Melting reduces snow amount 3115. C**** Ice advection is downward from layer 2 into layer 3 3116. 40 CMPRS = MIN (dSNdML*MELT1, MSI(I,J,1)-ACE1I-MELT1) 3117. FMSI1 = W0(I,J,2) + XSI1*CMPRS - XSI2*MELT1 3118. FMSI2 = W0(I,J,2) + CMPRS ! > 0 3119. IF(FMSI1.lt.0.) FHSI1 = HSI(I,J,2)*FMSI1*R2 / MSI(I,J,1) 3120. IF(FMSI1.ge.0.) FHSI1 = - ELHM*FMSI1 3121. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 3122. MSI(I,J,1) = MSI(I,J,1) - MELT1 - CMPRS 3123. GO TO 310 3124. C**** 3125. C**** No snow or ice melts in layer 1 3126. C**** 3127. 100 IF(W0(I,J,2) .gt. 0.) GO TO 300 3128. IF(MSI(I,J,1)-ACE1I+W0(I,J,2) .ge. 0.) GO TO 200 3129. C**** 3130. C**** All snow and some ice evaporates 3131. C**** Ice advection is upward into layer 2 from layer 3 3132. C**** 3133. FMSI1 = XSI1*(MSI(I,J,1)-ACE1I) + W0(I,J,2) ! < 0 3134. FMSI2 = MSI(I,J,1) - ACE1I + W0(I,J,2) ! < 0 3135. 110 FHSI1 = HSI(I,J,2)*FMSI1*R2 / MSI(I,J,1) 3136. FHSI2 = HSI(I,J,3)*FMSI2*R3 / MSI(I,J,2) 3137. MSI(I,J,1) = ACE1I 3138. IF(HSI(I,J,4)/ELHM+XSI4*MSI(I,J,2) .le. 0.) GO TO 120 3139. C**** Fluxes heat layer 4 to freezing point and melt some ice 3140. MELT4 = HSI(I,J,4)/ELHM + XSI4*MSI(I,J,2) 3141. MO(I,J,1) = MO(I,J,1) + FICE*MELT4*DXYP(J) 3142. CJ(J,22) = CJ(J,22) + FICE*MELT4*ZATMO(I,J) 3143. DJ(J,22) = DJ(J,22) - FICE*MELT4*ZATMO(I,J) 3144. C CJ(J,38) = CJ(J,38) + FICE*MELT4 ! in DIAGJ 3145. DJ(J,38) = DJ(J,38) - FICE*MELT4 3146. FMSI3 = XSI4*FMSI2 + XSI3*MELT4 3147. IF(FMSI3.le.0.) FHSI3 = - ELHM*FMSI3 3148. IF(FMSI3.gt.0.) FHSI3 = HSI(I,J,3)*FMSI3*R3 / MSI(I,J,2) 3149. MSI(I,J,2) = MSI(I,J,2) + FMSI2 - MELT4 3150. GO TO 330 3151. C**** No ice melts in layer 4 3152. C FMSI3 = XSI4*FMSI2 3153. 120 FHSI3 = HSI(I,J,4)*FMSI2 / MSI(I,J,2) 3154. MSI(I,J,2) = MSI(I,J,2) + FMSI2 3155. GO TO 330 3156. C**** 3157. C**** Some snow evaporates 3158. C**** No ice advection between layers 2 and 3 3159. C**** 3160. C FMSI1 = XSI2*W0(I,J,2) ! < 0 3161. C FMSI2 = 0 3162. 200 FHSI1 = HSI(I,J,2)*W0(I,J,2) / MSI(I,J,1) 3163. MSI(I,J,1) = MSI(I,J,1) + W0(I,J,2) 3164. HSI(I,J,1) = HSI(I,J,1) - FHSI1 3165. HSI(I,J,2) = HSI(I,J,2) + FHSI1 3166. IF(HSI(I,J,4)/ELHM+XSI4*MSI(I,J,2) .le. 0.) GO TO 400 3167. C**** Fluxes heat layer 4 to freezing point and melt some ice 3168. MELT4 = HSI(I,J,4)/ELHM + XSI4*MSI(I,J,2) 3169. MO(I,J,1) = MO(I,J,1) + FICE*MELT4*DXYP(J) 3170. CJ(J,22) = CJ(J,22) + FICE*MELT4*ZATMO(I,J) 3171. DJ(J,22) = DJ(J,22) - FICE*MELT4*ZATMO(I,J) 3172. C CJ(J,38) = CJ(J,38) + FICE*MELT4 ! in DIAGJ 3173. DJ(J,38) = DJ(J,38) - FICE*MELT4 3174. C FMSI3 = XSI3*MELT4 ! > 0 3175. FHSI3 = HSI(I,J,3)*MELT4 / MSI(I,J,2) 3176. MSI(I,J,2) = MSI(I,J,2) - MELT4 3177. HSI(I,J,3) = HSI(I,J,3) - FHSI3 3178. HSI(I,J,4) = HSI(I,J,4) + FHSI3 3179. GO TO 400 3180. C**** 3181. C**** Dew increases ice amount, advect ice downward 3182. C**** Ice advection is downward from layer 2 into layer 3 3183. C**** 3184. C FMSI1 = W0(I,J,2) ! > 0 3185. 300 FMSI2 = W0(I,J,2) ! > 0 3186. FHSI1 = HSI(I,J,1)*FMSI2 / (XSI1*MSI(I,J,1)+W0(I,J,2)) 3187. FHSI2 = HSI(I,J,2)*FMSI2*R2 / MSI(I,J,1) 3188. 310 IF(HSI(I,J,4)/ELHM+XSI4*MSI(I,J,2) .le. 0.) GO TO 320 3189. C**** Fluxes heat layer 4 to freezing point and melt some ice 3190. MELT4 = HSI(I,J,4)/ELHM + XSI4*MSI(I,J,2) 3191. MO(I,J,1) = MO(I,J,1) + FICE*MELT4*DXYP(J) 3192. CJ(J,22) = CJ(J,22) + FICE*MELT4*ZATMO(I,J) 3193. DJ(J,22) = DJ(J,22) - FICE*MELT4*ZATMO(I,J) 3194. C CJ(J,38) = CJ(J,38) + FICE*MELT4 ! in DIAGJ 3195. DJ(J,38) = DJ(J,38) - FICE*MELT4 3196. C FMSI3 = XSI4*FMSI2 + XSI3*MELT4 ! > 0 3197. FHSI3 = HSI(I,J,3)*((XSI4/XSI3)*FMSI2+MELT4) / MSI(I,J,2) 3198. MSI(I,J,2) = MSI(I,J,2) + FMSI2 - MELT4 3199. GO TO 330 3200. C**** No ice melts in layer 4 3201. C FMSI3 = XSI4*FMSI2 ! > 0 3202. 320 FHSI3 = HSI(I,J,3)*(XSI4/XSI3)*FMSI2 / MSI(I,J,2) 3203. MSI(I,J,2) = MSI(I,J,2) + FMSI2 3204. 330 HSI(I,J,1) = HSI(I,J,1) - FHSI1 3205. HSI(I,J,2) = HSI(I,J,2) + (FHSI1 - FHSI2) 3206. HSI(I,J,3) = HSI(I,J,3) + (FHSI2 - FHSI3) 3207. HSI(I,J,4) = HSI(I,J,4) + FHSI3 3208. 400 continue 3209. RETURN 3210. END 4000. 4001. SUBROUTINE OCLIM 4002. C**** 4003. C**** OCLIM reads in monthly climatological lake or ocean values of 4004. C**** temperature, ice fraction, and ice thickness, fits them by a 4005. C**** quadratic or piecewise linear distribution, and evaluates the 4006. C**** fitted quantities each day. 4007. C**** 4008. C**** Output: RSI = horizontal ratio of lake ice to lake 4009. C**** GZM (C) = lake surface temperature 4010. C**** MSI (kg/m^2)= sea ice mass 4011. C**** HSI (J/m^2) = lake ice enthalpy minus latent energy 4012. C**** MO (kg) = liquid lake mass above sill depth 4013. C**** G0M (J) = liquid lake enthalpy 4014. C**** 4015. INCLUDE 'C070.COM' 4016. PARAMETER (ACE2SI=273., ACESI=ACE1I+ACE2SI, RHOI=910., FLEAD=.1d0) 4017. INTEGER*4 JDPERM(12) 4018. CHARACTER TITLE*80, CMON(12)*4 4019. COMMON /OCLICB/ 4020. * AOST(IM,JM),EOST1(IM,JM),EOST0(IM,JM),BOST(IM,JM),COST(IM,JM), 4021. * ARSI(IM,JM),ERSI1(IM,JM),ERSI0(IM,JM),BRSI(IM,JM),CRSI(IM,JM), 4022. * AMSI(IM,JM),EMSI1(IM,JM),EMSI0(IM,JM),BMSI(IM,JM),CMSI(IM,JM), 4023. * KRSI(IM,JM), KMSI(IM,JM) 4024. C**** 4025. C**** OCLICB AOST monthly mean ocean surface temperature (C) 4026. C**** BOST first order moment of ocean surface temperature 4027. C**** COST second order moment of ocean surface temperature 4028. C**** EOST0 ocean surface temperature at beginning of month 4029. C**** EOST1 ocean surface temperature at end of month 4030. C**** 4031. C**** ARSI monthly mean ratio of sea (lake) ice to water 4032. C**** BRSI first order moment of ratio of sea ice to water 4033. C**** or .5*[(ERSI0-Limit)^2 + (ERSI1-Limit)^2] / (ARSI-Limit) 4034. C**** CRSI second order moment of ratio of sea ice to water 4035. C**** ERSI0 ratio of sea ice to water at beginning of month 4036. C**** ERSI1 ratio of sea ice to water at end of month 4037. C**** KRSI -1: continuous piecewise linear fit at lower limit 4038. C**** 0: quadratic fit 4039. C**** 1: continuous piecewise linear fit at upper limit 4040. C**** 4041. C**** AMSI monthly mean total sea (lake) ice mass (kg/m^2) 4042. C**** BMSI first order moment of total sea ice mass 4043. C**** or .5*[(EMSI0-Limit)^2 + (EMSI1-Limit)^2] / (AMSI-Limit) 4044. C**** CMSI second order moment of total sea ice mass 4045. C**** EMSI0 total sea ice mass at beginning of month 4046. C**** EMSI1 total sea ice mass at end of month 4047. C**** KMSI -1: continuous piecewise linear fit at lower limit 4048. C**** 0: quadratic fit 4049. C**** 4050. DATA JDPERM /31,28,31, 30,31,30, 31,31,30, 31,30,31/ 4051. DATA CMON /'Jan','Feb','Mar','Apr','May','Jun', 4052. * 'Jul','Aug','Sep','Oct','Nov','Dec'/ 4053. DATA JMLAST /0/ 4054. C**** Kill coding for determination of JMONX, replace JMONX with JMON 4055. C**** when JMON becomes integer 4056. DO 10 KMONTH=1,12 4057. 10 IF(CMON(KMONTH).eq.JMON) GO TO 20 4058. STOP 'OCLIM 10:' 4059. 20 JMONX = KMONTH 4060. C**** 4061. IF(JMLAST.eq.JMONX) GO TO 400 4062. IF(JMLAST.gt.0) GO TO 200 4063. C**** 4064. C**** At restart, read beginning of month quantities EOST0, ERSI0, EMSI0 4065. C**** 4066. JMLAST = JMONX - 1 4067. IF(JMLAST.le.0) JMLAST = 12 4068. DO 110 KMONTH=1,JMLAST-1 4069. READ (21) 4070. READ (31) 4071. 110 READ (32) 4072. CALL READR4 (21,IM*JM*2,EOST1,EOST1) ! EOST0 follows EOST1 4073. CALL READR4 (31,IM*JM*2,ERSI1,ERSI1) ! ERSI0 follows ERSI1 4074. CALL READR4 (32,IM*JM*2,EMSI1,EMSI1) ! EMSI0 follows EMSI1 4075. C**** Convert ERSI0 from percentage to fraction 4076. DO 120 I=1,IM*JM 4077. 120 ERSI0(I,1) = ERSI0(I,1)/100. 4078. GO TO 220 4079. C**** 4080. C**** New month encountered 4081. C**** 4082. C**** Copy EOST1 to EOST0, ERSI1 to ERSI0, EMSI1 to EMSI0 4083. 200 DO 210 I=1,IM*JM 4084. EOST0(I,1) = EOST1(I,1) 4085. ERSI0(I,1) = ERSI1(I,1) 4086. 210 EMSI0(I,1) = EMSI1(I,1) 4087. C**** Read data for new month: A0ST, EOST1, ARSI, ERSI1, AMSI, EMSI1 4088. 220 IF(JMLAST.eq.12) REWIND (21) 4089. IF(JMLAST.eq.12) REWIND (31) 4090. IF(JMLAST.eq.12) REWIND (32) 4091. CALL READR4 (21,IM*JM*2,AOST,AOST) ! EOST1 follows AOST 4092. CALL READR4 (31,IM*JM*2,ARSI,ARSI) ! ERSI1 follows ARSI 4093. CALL READR4 (32,IM*JM*2,AMSI,AMSI) ! EMSI1 follows AMSI 4094. JMLAST = JMONX 4095. C**** Convert ARSI and ERSI1 from percentage to fraction 4096. DO 230 I=1,IM*JM 4097. ARSI(I,1) = ARSI(I,1)/100. 4098. 230 ERSI1(I,1) = ERSI1(I,1)/100. 4099. C**** 4100. C**** Determine whether new month will use quadratic fit or 4101. C**** continuous piecewise linear fit 4102. C**** 4103. DO 320 J=1,JM 4104. IMAX=IM 4105. IF(J.eq.1 .or. J.eq.JM) IMAX=1 4106. DO 320 I=1,IMAX 4106.5 IF(FLAKE(I,J).le.0.) GO TO 320 4107. C**** OST always uses quadratic fit 4108. C**** OST(t) = AOST + BOST*t + COST*(t*t - 1/12) , -.5 < t < .5 4109. BOST(I,J) = EOST1(I,J) - EOST0(I,J) 4110. COST(I,J) = 3.*EOST0(I,J) + 3.*EOST1(I,J) - 6.*AOST(I,J) 4111. C**** Assume RSI will use quadratic fit 4112. BRSI(I,J) = ERSI1(I,J) - ERSI0(I,J) 4113. CRSI(I,J) = 3.*ERSI0(I,J) + 3.*ERSI1(I,J) - 6.*ARSI(I,J) 4114. KRSI(I,J) = 0 4115. IF(ABS(CRSI(I,J)) .le. ABS(BRSI(I,J))) GO TO 310 4116. RSICSQ = CRSI(I,J)* ! RSI*CRSI^2 at apex of parabola 4117. * (ARSI(I,J)*CRSI(I,J) - .25*BRSI(I,J)**2 - CRSI(I,J)**2/12.) 4118. IF(RSICSQ.lt.0.) then 4119. C**** RSI uses piecewise linear fit because quadratic fit at apex < 0 4120. KRSI(I,J) = -1 4120.5 BRSI(I,J) = 1d20 4120.6 IF(ARSI(I,J).gt.0) 4121. * BRSI(I,J) = .5*(ERSI0(I,J)**2 + ERSI1(I,J)**2) / ARSI(I,J) 4122. ELSEIF(RSICSQ.gt.CRSI(I,J)**2) then 4123. C**** RSI uses piecewise linear fit because quadratic fit at apex > 1 4124. KRSI(I,J) = 1 4124.5 BRSI(I,J) = -1d20 4124.6 IF(ARSI(I,J).lt.1) 4125. * BRSI(I,J) = .5*((ERSI0(I,J)-1.)**2 + (ERSI1(I,J)-1.)**2) / 4126. / (ARSI(I,J)-1.) 4127. endif 4128. C**** Assume MSI will use quadratic fit 4129. 310 BMSI(I,J) = EMSI1(I,J) - EMSI0(I,J) 4130. CMSI(I,J) = 3.*EMSI0(I,J) + 3.*EMSI1(I,J) - 6.*AMSI(I,J) 4131. KMSI(I,J) = 0 4132. IF(ABS(CMSI(I,J)) .le. ABS(BMSI(I,J))) GO TO 320 4133. MSICSQ = CMSI(I,J)* ! MSI*CMSI^2 at apex of parabola 4134. * (AMSI(I,J)*CMSI(I,J) - .25*BMSI(I,J)**2 - CMSI(I,J)**2/12.) 4135. IF(MSICSQ.lt.ACESI*CMSI(I,J)**2) then 4136. C**** MSI uses piecewise linear fit because quadratic fit < ACESI 4137. KMSI(I,J) = -1 4137.5 BMSI(I,J) = 1d20 4137.6 IF(AMSI(I,J).gt.ACESI) 4138. * BMSI(I,J) = .5*((EMSI0(I,J)-ACESI)**2 + (EMSI1(I,J)-ACESI)**2) 4139. / / (AMSI(I,J)-ACESI) 4140. endif 4141. 320 continue 4142. C**** 4143. C**** Calculate OST, RSI and MSI for current day 4144. C**** 4145. 400 DO 540 J=1,JM 4146. IMAX=IM 4147. IF(J.eq.1 .or. J.eq.JM) IMAX=1 4148. DO 540 I=1,IMAX 4149. IF(FLAKE(I,J).le.0.) GO TO 540 4150. TIME = (JDATE-.5d0)/JDPERM(JMONX) - .5 ! -.5 < TIME < .5 4151. C**** OST always uses quadratic fit 4152. GZM(I,J,1) = AOST(I,J) + BOST(I,J)*TIME + 4153. + COST(I,J)*(TIME*TIME - 1./12.d0) 4154. IF(KRSI(I,J)) 410,430,420 4155. C**** RSI uses piecewise linear fit because quadratic fit at apex < 0 4156. 410 IF(ERSI0(I,J)-BRSI(I,J)*(TIME+.5) .gt. 0.) then 4157. RSINEW = ERSI0(I,J) - BRSI(I,J)*(TIME+.5) ! TIME < T0 4158. ELSEIF(ERSI1(I,J)-BRSI(I,J)*(.5-TIME) .gt. 0.) then 4159. RSINEW = ERSI1(I,J) - BRSI(I,J)*(.5-TIME) ! T1 < TIME 4160. ELSE 4161. RSINEW = 0. ! T0 < TIME < T1 4162. endif 4163. GO TO 440 4164. C**** RSI uses piecewise linear fit because quadratic fit at apex > 1 4165. 420 IF(ERSI0(I,J)-BRSI(I,J)*(TIME+.5) .lt. 1.) then 4166. RSINEW = ERSI0(I,J) - BRSI(I,J)*(TIME+.5) ! TIME < T0 4167. ELSEIF(ERSI1(I,J)-BRSI(I,J)*(.5-TIME) .lt. 1.) then 4168. RSINEW = ERSI1(I,J) - BRSI(I,J)*(.5-TIME) ! T1 < TIME 4169. ELSE 4170. RSINEW = 1. ! T0 < TIME < T1 4171. endif 4172. GO TO 440 4173. C**** RSI uses quadratic fit 4174. 430 RSINEW = ARSI(I,J) + BRSI(I,J)*TIME + 4175. + CRSI(I,J)*(TIME*TIME - 1./12.d0) 4176. 440 IF(KMSI(I,J).ge.0) GO TO 450 4177. C**** MSI uses piecewise linear fit because quadratic fit < ACESI 4178. IF(EMSI0(I,J)-BMSI(I,J)*(TIME+.5) .gt. ACESI) then 4179. MSINEW = EMSI0(I,J) - BMSI(I,J)*(TIME+.5) ! TIME < T0 4180. ELSEIF(EMSI1(I,J)-BMSI(I,J)*(.5-TIME) .gt. ACESI) then 4181. MSINEW = EMSI1(I,J) - BMSI(I,J)*(.5-TIME) ! T1 < TIME 4182. ELSE 4183. MSINEW = ACESI ! T0 < TIME < T1 4184. endif 4185. GO TO 500 4186. C**** MSI uses quadratic fit 4187. 450 MSINEW = AMSI(I,J) + BMSI(I,J)*TIME + 4188. + CMSI(I,J)*(TIME*TIME - 1./12.d0) 4189. C**** 4190. C**** Changes in ice values add or subtract liquid water to lakes 4191. C**** 4192. 500 IF(RSI(I,J).le.0.) GO TO 530 4193. IF(RSINEW.gt.0.) GO TO 510 4194. C**** All lake ice has melted: RSINEW = 0, RSIold > 0 4195. MO(I,J,1) = MO(I,J,1) + FWATER(I,J)*RSI(I,J)* 4196. * (MSI(I,J,1)+MSI(I,J,2))*DXYP(J) 4197. G0M(I,J,1) = G0M(I,J,1) + FWATER(I,J)*RSI(I,J)* 4198. * (HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4))*DXYP(J) 4199. DJ(J,36) = DJ(J,36) - FWATER(I,J)*RSI(I,J)* 4200. * (MSI(I,J,1)+MSI(I,J,2)) 4201. DJ(J,45) = DJ(J,45) - FWATER(I,J)*RSI(I,J)* 4202. * (HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4)) 4203. RSI(I,J) = 0. 4204. MSI(I,J,1) = ACE1I 4205. MSI(I,J,2) = ACESI-ACE1I 4206. HSI(I,J,1) = - ELHM*XSI1*ACE1I 4207. HSI(I,J,2) = - ELHM*XSI2*ACE1I 4208. HSI(I,J,3) = - ELHM*XSI3*(ACESI-ACE1I) 4209. HSI(I,J,4) = - ELHM*XSI4*(ACESI-ACE1I) 4210. GO TO 540 4211. C**** Both current day and prior day have sea ice: RSINEW > 0, RSIold >0 4212. 510 RSINEW = MIN(RSINEW,1.-FLEAD*RHOI/MSINEW) 4213. MSI2N = MSINEW - ACE1I 4214. C**** Calcualte mass and heat fluxes of sea ice at edges 3 and 4 4215. FMSI4 = MSI(I,J,2) - MSI2N 4216. C FMSI3 = XSI3*FMSI4 4217. FHSI3 = HSI(I,J,3)*FMSI4 / MSI(I,J,2) ! FMSI4 > 0 4218. FHSI4 = HSI(I,J,4)*FMSI4 / MSI(I,J,2)/XSI4 4219. IF(FMSI4.lt.0.) then 4220. FHSI3 = HSI(I,J,4)*FMSI4*(XSI3/XSI4) / MSI(I,J,2) ! FMSI4 < 0 4221. FHSI4 = - ELHM*FMSI4 4222. endif 4223. IF(FMSI4.lt.-XSI4*MSI2N) ! FMSI4 < - XSI4*MSI2N < 0 4224. * FHSI3 = - HSI(I,J,4) - ELHM*(FMSI4+XSI4*MSI2N) 4225. IF(FMSI4.gt.XSI4*MSI(I,J,2)) ! FMSI4 > XSI4*MSI(I,J,2) > 0 4226. * FHSI4 = HSI(I,J,4) + 4227. + HSI(I,J,3)*(FMSI4-XSI4*MSI(I,J,2)) / MSI(I,J,2)/XSI3 4228. IF(RSINEW.gt.RSI(I,J)) GO TO 520 4229. C**** Some lake ice has melted: RSINEW < RSIold, usually FMSI4 > 0 4230. MO(I,J,1) = MO(I,J,1) + FWATER(I,J)*(RSINEW*FMSI4 + 4231. + (RSI(I,J)-RSINEW)*(MSI(I,J,1)+MSI(I,J,2)))*DXYP(J) 4232. G0M(I,J,1) = G0M(I,J,1) + FWATER(I,J)*(RSINEW*FHSI4 + 4233. + (RSI(I,J)-RSINEW)*(HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4))) 4234. * *DXYP(J) 4235. DJ(J,36) = DJ(J,36) - FWATER(I,J)*(RSINEW*FMSI4 + 4236. + (RSI(I,J)-RSINEW)*(MSI(I,J,1)+MSI(I,J,2))) 4237. DJ(J,45) = DJ(J,45) - FWATER(I,J)*(RSINEW*FHSI4 + 4238. + (RSI(I,J)-RSINEW)*(HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4))) 4239. RSI(I,J) = RSINEW 4240. MSI(I,J,2) = MSI2N 4241. HSI(I,J,3) = HSI(I,J,3) - FHSI3 4242. HSI(I,J,4) = HSI(I,J,4) + (FHSI3 - FHSI4) 4243. GO TO 540 4244. C**** Additional lake ice has been frozen: RSINEW > RSIold 4245. 520 FREZ12 = (RSINEW-RSI(I,J))*ACE1I 4246. FREZ34 = (RSINEW-RSI(I,J))*MSI2N 4247. MO(I,J,1) = MO(I,J,1) + FWATER(I,J)*(RSI(I,J)*FMSI4 + 4248. + (RSI(I,J)-RSINEW)*MSINEW)*DXYP(J) 4249. G0M(I,J,1) = G0M(I,J,1) + FWATER(I,J)*(RSI(I,J)*FHSI4 - 4250. - (RSI(I,J)-RSINEW)*MSINEW*ELHM)*DXYP(J) 4251. DJ(J,37) = DJ(J,37) - FWATER(I,J)*(RSI(I,J)*FMSI4 + 4252. + (RSI(I,J)-RSINEW)*MSINEW) 4253. DJ(J,46) = DJ(J,46) - FWATER(I,J)*(RSI(I,J)*FHSI4 - 4254. - (RSI(I,J)-RSINEW)*MSINEW*ELHM) 4255. MSI(I,J,1) = (RSI(I,J)*MSI(I,J,1)+FREZ12) / RSINEW 4256. MSI(I,J,2) = MSI2N 4257. HSI(I,J,1) = (RSI(I,J)* HSI(I,J,1) - ELHM*XSI1*FREZ12) / RSINEW 4258. HSI(I,J,2) = (RSI(I,J)* HSI(I,J,2) - ELHM*XSI2*FREZ12) / RSINEW 4259. HSI(I,J,3) = (RSI(I,J)*(HSI(I,J,3)-FHSI3) - ELHM*XSI3*FREZ34) / 4260. / RSINEW 4261. HSI(I,J,4) = (RSI(I,J)*(HSI(I,J,4)+FHSI3-FHSI4) - 4262. - ELHM*XSI4*FREZ34) / RSINEW 4263. RSI(I,J) = RSINEW 4264. GO TO 540 4265. 530 IF(RSINEW.le.0.) GO TO 540 4266. C**** First time lake ice has been frozen: RSINEW > 0, RSIold = 0 4267. MO(I,J,1) = MO(I,J,1) - FWATER(I,J)*RSINEW*MSINEW*DXYP(J) 4268. G0M(I,J,1) = G0M(I,J,1) + FWATER(I,J)*RSINEW*MSINEW*DXYP(J)*ELHM 4269. DJ(J,37) = DJ(J,37) + FWATER(I,J)*RSINEW*MSINEW 4270. DJ(J,46) = DJ(J,46) - FWATER(I,J)*RSINEW*MSINEW*ELHM 4271. RSI(I,J) = RSINEW 4272. MSI(I,J,1) = ACE1I 4273. MSI(I,J,2) = MSINEW-ACE1I 4274. HSI(I,J,1) = - ELHM*XSI1*ACE1I 4275. HSI(I,J,2) = - ELHM*XSI2*ACE1I 4276. HSI(I,J,3) = - ELHM*XSI3*(MSINEW-ACE1I) 4277. HSI(I,J,4) = - ELHM*XSI4*(MSINEW-ACE1I) 4278. 540 continue 4279. C**** 4280. C**** Replicate polar values to all longitudes 4281. C**** 4291. RETURN 4292. END 5000. 5001. SUBROUTINE RIVERF 5002. C**** 5003. C**** RIVERF transports lake water from each GCM grid box to its 5004. C**** downstream neighbor according to the river direction file. 5005. C**** 5006. INCLUDE 'C070.COM' 5008. REAL*8 XK(8),YK(8) 5009. COMMON /RIVRCB/ RATE(IM,JM), IFLOW(IM,JM),JFLOW(IM,JM), 5010. * IDPOLE,JDPOLE, KDIREC(IM,JM) 5011. COMMON /WORK01/ FLOW(IM,JM),EFLOW(IM,JM) 5012. C**** Location of river entry in downstream ocean box 5013. DATA XK /-.5, .0, .5, .5, .5, .0,-.5,-.5/, 5014. * YK /-.5,-.5,-.5, .0, .5, .5, .5, .0/ 5015. C**** 5016. C**** OCENCB MO Liquid lake mass above sill depth (kg) 5017. C**** G0M Liquid lake enthalpy above sill depth (J) 5018. C**** GZM Lake surface temperature (C) 5019. C**** 5020. C**** RIVRCB RATE river flow rate 5021. C**** IFLOW,JFLOW downstream river flow grid box 5022. C**** KDIREC river direction (0 to 8) 5118. C**** 5119. C**** Calculate net mass and energy changes due to river flow 5120. C**** 5121. DO 410 I=1,IM*(JM-1) 5122. FLOW(I,1) = 0. 5123. 410 EFLOW(I,1) = 0. 5124. DO 490 JU=2,JM-1 5125. DO 490 IU=1,IM 5126. IF(KDIREC(IU,JU).le.0) GO TO 490 5127. JD=JFLOW(IU,JU) 5128. ID=IFLOW(IU,JU) 5129. IF(MO(IU,JU,1).le.0.) GO TO 490 5130. DMM = MO(IU,JU,1)*RATE(IU,JU) 5131. C DGM = G0M(IU,JU,1)*DMM/MO(IU,JU,1) 5132. DGM = GZM(IU,JU,1)*DMM*SHCW 5133. FLOW(IU,JU) = FLOW(IU,JU) - DMM 5134. EFLOW(IU,JU) = EFLOW(IU,JU) - DGM 5134.5 CJ(JU,54) = CJ(JU,54) - BYDXYP(JU)* DMM 5134.6 CJ(JU,25) = CJ(JU,25) - BYDXYP(JU)*(DGM+GRAV*DMM*ZATMO(IU,JU)) 5135. AIJ(IU,JU,34) = AIJ(IU,JU,34) + DMM 5136. AIJ(IU,JU,20) = AIJ(IU,JU,20) + (DGM+GRAV*DMM*ZATMO(IU,JU)) 5138. IF(FOCEAN(ID,JD).gt.0.) GO TO 450 5139. FLOW(ID,JD) = FLOW(ID,JD) + DMM 5140. EFLOW(ID,JD) = EFLOW(ID,JD) + 5140.1 + DGM+GRAV*DMM*(ZATMO(IU,JU)-ZATMO(ID,JD)) 5140.5 CJ(JD,54) = CJ(JD,54) + BYDXYP(JD)* DMM 5140.6 CJ(JD,25) = CJ(JD,25) + BYDXYP(JD)*(DGM+GRAV*DMM*ZATMO(IU,JU)) 5141. GO TO 490 5142. C**** 5143. C**** Apply river mouth flow to ocean arrays 5144. C**** 5145. 450 K=KDIREC(IU,JU) 5146. GXM(ID,JD,1) = GXM(ID,JD,1) + 3.*XK(K)*DGM + 5147. * (GXM(ID,JD,1)*(.5-1.5*XK(K)*XK(K))-3.*XK(K)*G0M(ID,JD,1))* 5147.1 * (DMM*BYDXYP(JD)/MO(ID,JD,1)) 5148. GYM(ID,JD,1) = GYM(ID,JD,1) + 3.*YK(K)*DGM + 5149. * (GYM(ID,JD,1)*(.5-1.5*YK(K)*YK(K))-3.*YK(K)*G0M(ID,JD,1))* 5149.1 * (DMM*BYDXYP(JD)/MO(ID,JD,1)) 5150. SXM(ID,JD,1) = SXM(ID,JD,1) + (DMM*BYDXYP(JD)/MO(ID,JD,1))* 5151. * (SXM(ID,JD,1)*(.5-1.5*XK(K)*XK(K))-3.*XK(K)*S0M(ID,JD,1)) 5152. SYM(ID,JD,1) = SYM(ID,JD,1) + (DMM*BYDXYP(JD)/MO(ID,JD,1))* 5153. * (SYM(ID,JD,1)*(.5-1.5*YK(K)*YK(K))-3.*YK(K)*S0M(ID,JD,1)) 5154. MO(ID,JD,1) = MO(ID,JD,1) + BYDXYP(JD)* DMM 5155. G0M(ID,JD,1) = G0M(ID,JD,1) + 5155.1 + DGM+GRAV*DMM*(ZATMO(IU,JU)-ZATMO(ID,JD)) 5155.5 EJ(JD,54) = EJ(JD,54) + BYDXYP(JD)* DMM 5155.6 EJ(JD,25) = EJ(JD,25) + BYDXYP(JD)*(DGM+GRAV*DMM*ZATMO(IU,JU)) 5156. 490 continue 5157. C**** 5158. C**** Calculate river flow at the South Pole 5159. C**** 5160. FLOW(1,1) = FLOW(1,1)/IM 5161. EFLOW(1,1) = EFLOW(1,1)/IM 5162. IF(MO(1,1,1).le.0.) GO TO 600 5163. DMM = MO(1,1,1)*RATE(1,1) 5164. C DGM = G0M(1,1,1)*DMM/MO(1,1,1) 5165. DGM = GZM(1,1,1)*DMM*SHCW 5166. FLOW(1,1) = FLOW(1,1) - DMM 5167. EFLOW(1,1) = EFLOW(1,1) - DGM 5167.5 CJ(1,54) = CJ(1,54) - BYDXYP(1)* DMM 5167.6 CJ(1,25) = CJ(1,25) - BYDXYP(1)*(DGM+GRAV*DMM*ZATMO(1,1)) 5168. AIJ(1,1,34) = AIJ(1,1,34) + DMM 5169. AIJ(1,1,20) = AIJ(1,1,20) + (DGM+GRAV*DMM*ZATMO(1,1)) 5170. FLOW(IDPOLE,JDPOLE) = FLOW(IDPOLE,JDPOLE) + IM*DMM 5171. EFLOW(IDPOLE,JDPOLE) = EFLOW(IDPOLE,JDPOLE) + 5171.1 + IM*(DGM+GRAV*DMM*(ZATMO(1,1)-ZATMO(IDPOLE,JDPOLE))) 5171.5 CJ(JDPOLE,54) = CJ(JDPOLE,54) + IM*DMM*BYDXYP(JDPOLE) 5171.6 CJ(JDPOLE,25) = CJ(JDPOLE,25) + 5171.7 + IM*(DGM+GRAV*DMM*ZATMO(1,1))*BYDXYP(JDPOLE) 5172. C**** 5173. C**** Apply net river flow to continental reservoirs 5174. C**** 5175. 600 DO 610 J=2,JM-1 5176. DO 610 I=1,IM 5177. IF(FOCEAN(I,J).gt.0.) GO TO 610 5178. MO(I,J,1) = MO(I,J,1) + FLOW(I,J) 5179. G0M(I,J,1) = G0M(I,J,1) + EFLOW(I,J) 5180. 610 continue 5181. MO(1,1,1) = MO(1,1,1) + FLOW(1,1) 5182. G0M(1,1,1) = G0M(1,1,1) + EFLOW(1,1) 5183. RETURN 5184. C**** 5185. 910 FORMAT (A72) 5186. 911 FORMAT (72A1) 5187. END 5500. 5501. SUBROUTINE RIVER0 5502. C**** 5503. C**** RIVER0 is called from INPUT to do the following: 5504. C**** 1. read in the river direction file from unit 54 5505. C**** 2. calculate the downstream grid box for each continental box 5506. C**** 3. calculate the river flow RATE = DTS * SPEED * dHORZ 5507. C**** 4. determine whether ocean box receives flow from ice shelf 5508. C**** 5509. INCLUDE 'C070.COM' 5510. LOGICAL*4 QSHELF 5511. CHARACTER TITLEI*80, CDIREC*1 5512. COMMON /RIVRCB/ RATE(IM,JM), IFLOW(IM,JM),JFLOW(IM,JM), 5513. * IDPOLE,JDPOLE, KDIREC(IM,JM), QSHELF(IM,JM) 5514. COMMON /WORK01/ DHORZ(IM,JM), CDIREC(IM,JM) 5515. C**** 5516. C**** RIVRCB RATE river flow rate 5517. C**** IFLOW,JFLOW downstream river flow grid box 5518. C**** KDIREC river direction (0 to 8) 5519. C**** QSHELF whether ocean box receives flow from ice shelf 5520. C**** 5521. C**** Read in CDIREC: Number = river direction (0 to 8) 5522. C**** s = ocean receives flow from ice shelf 5523. C**** 5524. READ (54,910) TITLEI 5525. WRITE (6,*) 'River Direction file read on unit 54: ',TITLEI 5526. READ (54,910) 5527. DO 110 I72=72,IM,72 5528. DO 110 J=JM,1,-1 5529. 110 READ (54,911) (CDIREC(I,J),I=I72-71,I72) 5530. CLOSE (54) 5531. C**** Create integral direction array KDIREC from CDIREC 5532. DO 120 J=1,JM 5533. DO 120 I=1,IM 5534. QSHELF(I,J) = CDIREC(I,J) .eq. 's' 5535. KDIREC(I,J) = 0. 5536. IF(FOCEAN(I,J).le.0.) KDIREC(I,J) = ICHAR(CDIREC(I,J)) - 48 5537. IF(KDIREC(I,J).ge.0 .and. KDIREC(I,J).le.8) GO TO 120 5538. WRITE (6,*) 'Land box has invalid river direction: I,J,KDIREC=', 5539. * I,J,KDIREC(I,J) 5540. STOP 'RIVER0 120:' 5541. 120 continue 5542. C**** 5543. C**** From each box calculate the downstream river box 5544. C**** 5545. DO 290 J=2,JM-1 5546. DO 290 I=1,IM 5547. GO TO (210,220,230,240,250,260,270,280),KDIREC(I,J) 5548. IFLOW(I,J) = I 5549. JFLOW(I,J) = J 5550. DHORZ(I,J) = 0. 5551. GO TO 290 5552. 210 IFLOW(I,J) = I+1 5553. JFLOW(I,J) = J+1 5554. DHORZ(I,J) = SQRT(DXV(J)*DXV(J)+DYV(J)*DYV(J)) 5555. IF(I.eq.IM) IFLOW(I,J) = 1 5556. GO TO 290 5557. 220 IFLOW(I,J) = I 5558. JFLOW(I,J) = J+1 5559. DHORZ(I,J) = DYV(J) 5560. GO TO 290 5561. 230 IFLOW(I,J) = I-1 5562. JFLOW(I,J) = J+1 5563. DHORZ(I,J) = SQRT(DXV(J)*DXV(J)+DYV(J)*DYV(J)) 5564. IF(I.eq.1) IFLOW(I,J) = IM 5565. GO TO 290 5566. 240 IFLOW(I,J) = I-1 5567. JFLOW(I,J) = J 5568. DHORZ(I,J) = DXP(J) 5569. IF(I.eq.1) IFLOW(I,J) = IM 5570. GO TO 290 5571. 250 IFLOW(I,J) = I-1 5572. JFLOW(I,J) = J-1 5573. DHORZ(I,J) = SQRT(DXV(J-1)*DXV(J-1)+DYV(J-1)*DYV(J-1)) 5574. IF(I.eq.1) IFLOW(I,J) = IM 5575. GO TO 290 5576. 260 IFLOW(I,J) = I 5577. C IF(J.eq.2) IFLOW(I,J) = 1 ! corrected under "South Pole .." 5578. JFLOW(I,J) = J-1 5579. DHORZ(I,J) = DYV(J-1) 5580. GO TO 290 5581. 270 IFLOW(I,J) = I+1 5582. JFLOW(I,J) = J-1 5583. DHORZ(I,J) = SQRT(DXV(J-1)*DXV(J-1)+DYV(J-1)*DYV(J-1)) 5584. IF(I.eq.IM) IFLOW(I,J) = 1 5585. GO TO 290 5586. 280 IFLOW(I,J) = I+1 5587. JFLOW(I,J) = J 5588. DHORZ(I,J) = DXP(J) 5589. IF(I.eq.IM) IFLOW(I,J) = 1 5590. 290 continue 5591. C**** South Pole is a special case 5592. DO 295 I=1,IM 5593. IF(KDIREC(I,1).eq.2) then 5594. IDPOLE = I 5595. JDPOLE = 2 5596. IFLOW(1,1) = I 5597. JFLOW(1,1) = 2 5598. DHORZ(1,1) = DYV(1) 5599. endif 5600. IF(KDIREC(I,2).eq.6) then 5601. IFLOW(I,2) = 1 5602. JFLOW(I,2) = 1 5603. endif 5604. 295 continue 5605. C**** 5606. C**** Calculate river flow RATE (per source time step) 5607. C**** 5608. SPEED0= .35d0 5609. SPMIN = .15d0 5610. SPMAX = 5. 5611. DZDH1 = .00005d0 5612. DO 310 JU=1,JM-1 5613. IMAX=IM 5614. IF(JU.eq.1) IMAX=1 5615. DO 310 IU=1,IMAX 5616. IF(KDIREC(IU,JU).le.0) GO TO 310 5617. JD=JFLOW(IU,JU) 5618. ID=IFLOW(IU,JU) 5619. DZDH = (ZATMO(IU,JU)-ZATMO(ID,JD)) / DHORZ(IU,JU) 5620. SPEED = SPEED0*DZDH/DZDH1 5621. IF(SPEED.lt.SPMIN) SPEED = SPMIN 5622. IF(SPEED.gt.SPMAX) SPEED = SPMAX 5623. RATE(IU,JU) = DTS*SPEED/DHORZ(IU,JU) 5624. 310 continue 5625. RETURN 5626. C**** 5627. 910 FORMAT (A72) 5628. 911 FORMAT (72A1) 5629. END 6000. 6001. SUBROUTINE PRECGI 6002. C**** 6003. C**** PRECGI adds the atmospheric precipitation to the glacial ice 6004. C**** 6005. INCLUDE 'C070.COM' 6006. PARAMETER (SNOMAX=91.66d0, dSNdRN=2., R2=1/XGI2) 6007. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 6008. C**** 6009. C**** GICECB MGI Mass of glacial ice and snow (kg/m^2) 6010. C**** HGI Enthalpy minus latent heat of glacial ice (J/m^2) 6011. C**** 6012. C**** OCENCB MO River and lake mass above sill depth (kg) 6013. C**** G0M River and lake energy above sill depth (J) 6014. C**** 6015. C**** FIXDCB FGICE Glacial ice fraction 6016. C**** 6017. C**** PRECCB PREC Precipitation from atmosphere (kg/m^2) 6018. C**** EPRE Energy of precipitation (J/m^2) 6019. C**** 6020. C**** Outside loop over J and I 6021. C**** 6022. DO 500 J=1,JM-1 6023. IMAX=IM 6024. IF(J.eq.1) IMAX=1 6025. DO 500 I=1,IMAX 6026. IF(PREC(I,J,1).le.0. .or. FGICE(I,J).le.0.) GO TO 500 6027. FICE = FGICE(I,J) 6028. C**** 6029. C**** Apply precipitation heat flux to HGI1 6030. C**** 6031. HGI(I,J,1) = HGI(I,J,1) + EPRE(I,J,1) 6032. IF(EPRE(I,J,1).le.-PREC(I,J,1)*ELHM) GO TO 300 6033. IF(EPRE(I,J,1).le.0.) GO TO 200 6034. C**** 6035. C**** All precipitation is rain above 0øC 6036. C**** Rain compresses snow amount into ice 6037. C**** Calculate snow and ice melted or rain frozen in layer 1 6038. C**** 6039. RAIN = PREC(I,J,1) 6040. IF(HGI(I,J,1)/ELHM+XGI1*MGI(I,J) .le. 0.) GO TO 140 6041. C**** Warm rain cools to 0øC and melts some snow or ice. 6042. MELT1 = HGI(I,J,1)/ELHM + XGI1*MGI(I,J) 6043. IF(MGI(I,J)-ACE1I .le. MELT1) GO TO 120 6044. C**** Rain and melting compresses snow into ice 6045. CMPRS = MIN (dSNdRN*(RAIN+MELT1), MGI(I,J)-ACE1I-MELT1) 6046. C FMGI1 = XGI1*CMPRS - XGI2*MELT1 6047. FHGI1 = - ELHM*(XGI1*CMPRS - XGI2*MELT1) 6048. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. RAIN+MELT1) GO TO 110 6049. C**** All rain and melt water is frozen in layer 2 6050. C FREZ2 = RAIN + MELT1 6051. FMGI4 = RAIN + MELT1 + CMPRS ! > 0 6052. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+RAIN+MELT1) 6053. MGI(I,J) = MGI(I,J) - MELT1 - CMPRS 6054. GO TO 400 6055. C**** Only some rain and melt water is frozen in layer 2 6056. 110 FREZ2 = -HGI(I,J,2)/ELHM-XGI2*MGI(I,J) 6057. RAIN2 = RAIN + MELT1 - FREZ2 6058. FMGI2 = CMPRS + FREZ2 ! > 0 6059. FHGI2 = - ELHM*FMGI2 6060. MGI(I,J) = MGI(I,J) - MELT1 - CMPRS 6061. GO TO 230 6062. C**** All snow and some ice has melted, what can be frozen in layer 2? 6063. C FMGI1 = XGI1*(MGI-ACE1I) - MELT1 ! < 0 6064. 120 FHGI1 = - ELHM*(XGI1*(MGI(I,J)-ACE1I)-MELT1) 6065. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. RAIN+MELT1) GO TO 130 6066. C**** All rain and melt water is frozen in layer 2 6067. C FREZ2 = RAIN + MELT1 6068. FMGI4 = MGI(I,J) - ACE1I + RAIN ! > 0 6069. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+RAIN+MELT1) 6070. MGI(I,J) = ACE1I 6071. GO TO 400 6072. C**** Only some rain and melt water is frozen in layer 2 6073. 130 FREZ2 = -HGI(I,J,2)/ELHM - XGI2*MGI(I,J) 6074. RAIN2 = RAIN + MELT1 - FREZ2 6075. FMGI2 = MGI(I,J) - ACE1I - MELT1 + FREZ2 6076. FHGI2 = - ELHM*FMGI2 6077. MGI(I,J) = ACE1I 6078. GO TO 230 6079. C**** Rain compresses snow into ice, some rain will freeze 6080. 140 CMPRS = MIN (dSNdRN*RAIN, MGI(I,J)-ACE1I) 6081. IF(-HGI(I,J,1)/ELHM-XGI1*MGI(I,J) .lt. RAIN) GO TO 150 6082. C**** All rain freezes in layer 1 6083. C FREZ1 = RAIN 6084. C FMGI1 = XGI1*CMPRS + FREZ1 ! > 0 6085. FMGI4 = CMPRS + RAIN ! > 0 6086. FHGI1 = HGI(I,J,1)*(XGI1*CMPRS+FREZ1) / (XGI1*MGI(I,J)+RAIN) 6087. FHGI2 = HGI(I,J,2)*FMGI4*R2 / MGI(I,J) 6088. MGI(I,J) = MGI(I,J) - CMPRS 6089. GO TO 400 6090. C**** Not all rain freezes in layer 1 6091. 150 FREZ1 = - HGI(I,J,1)/ELHM - XGI1*MGI(I,J) 6092. RAIN1 = RAIN - FREZ1 6093. C FMGI1 = XGI1*CMPRS + FREZ1 ! > 0 6094. FHGI1 = - ELHM*(XGI1*CMPRS + FREZ1) 6095. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. RAIN1) GO TO 160 6096. C**** All remaining rain freezes in layer 2 6097. C FREZ2 = RAIN1 6098. FMGI4 = CMPRS + RAIN ! > 0 6099. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+RAIN1) 6100. MGI(I,J) = MGI(I,J) - CMPRS 6101. GO TO 400 6102. C**** Not all remaining rain freezes in layer 2 6103. 160 FREZ2 = - HGI(I,J,2)/ELHM - XGI2*MGI(I,J) 6104. RAIN2 = RAIN1 - FREZ2 6105. FMGI2 = CMPRS + FREZ1 + FREZ2 ! > 0 6106. FHGI2 = - ELHM*FMGI2 6107. MGI(I,J) = MGI(I,J) - CMPRS 6108. GO TO 230 6109. C**** 6110. C**** Precipitation is a mixture of rain and snow at 0øC 6111. C**** Rain compresses snow into ice, most rain will freeze 6112. C**** Note that snow amount may exceed SNOMAX 6113. C**** 6114. 200 SNOW = -EPRE(I,J,1)/ELHM 6115. RAIN = PREC(I,J,1) - SNOW 6116. CMPRS = MIN (dSNdRN*RAIN, MGI(I,J)+SNOW-ACE1I) 6117. IF(-HGI(I,J,1)/ELHM-XGI1*MGI(I,J)-SNOW .lt. RAIN) GO TO 210 6118. C**** All rain freezes in layer 1 6119. C FREZ1 = RAIN 6120. FMGI1 = XGI1*CMPRS + XGI2*SNOW + RAIN ! > 0 6121. FMGI4 = CMPRS + RAIN ! > 0 6122. FHGI1 = HGI(I,J,1)*FMGI1 / (XGI1*MGI(I,J)+PREC(I,J,1)) 6123. FHGI2 = HGI(I,J,2)*FMGI4*R2 / MGI(I,J) 6124. MGI(I,J) = MGI(I,J) + (SNOW-CMPRS) 6125. GO TO 400 6126. C**** Not all rain freezes in layer 1 6127. 210 FREZ1 = - HGI(I,J,1)/ELHM - XGI1*MGI(I,J) - SNOW 6128. RAIN1 = RAIN - FREZ1 6129. C FMGI1 = XGI1*CMPRS + XGI2*SNOW + FREZ1 ! > 0 6130. FHGI1 = - ELHM*(XGI1*CMPRS + XGI2*SNOW + FREZ1) 6131. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. RAIN1) GO TO 220 6132. C**** All remaining rain freezes in layer 2 6133. C FREZ2 = RAIN1 6134. FMGI4 = CMPRS + RAIN ! > 0 6135. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+RAIN1) 6136. MGI(I,J) = MGI(I,J) + (SNOW-CMPRS) 6137. GO TO 400 6138. C**** Not all remaining rain freezes in layer 2 6139. 220 FREZ2 = - HGI(I,J,2)/ELHM - XGI2*MGI(I,J) 6140. RAIN2 = RAIN1 - FREZ2 6141. FMGI2 = CMPRS + FREZ1 + FREZ2 ! > 0 6142. FHGI2 = - ELHM*FMGI2 6143. MGI(I,J) = MGI(I,J) + (SNOW-CMPRS) 6144. 230 IF(-HGI(I,J,3)/ELHM-XGI3*ACE2GI .lt. RAIN2) GO TO 240 6145. C**** All remaining rain freezes in layer 3 6146. C FMGI3 = FMGI2 + RAIN2 ! > 0 6147. FMGI4 = FMGI2 + RAIN2 ! > 0 6148. FHGI3 = HGI(I,J,3)*FMGI4 / (XGI3*ACE2GI+RAIN2) 6149. GO TO 410 6150. C**** Not all remaining rain freezes in layer 3 6151. 240 FREZ3 = - HGI(I,J,3)/ELHM - XGI3*ACE2GI 6152. RAIN3 = RAIN2 - FREZ3 6153. FMGI3 = FMGI2 + FREZ3 6154. FHGI3 = - ELHM*FMGI3 6155. IF(-HGI(I,J,4)/ELHM-XGI4*ACE2GI .lt. RAIN3) GO TO 250 6156. C**** All remaining rain freezes in layer 4 6157. FMGI4 = FMGI3 + RAIN3 ! > 0 6158. FHGI4 = HGI(I,J,4)*FMGI4 / (XGI4*ACE2GI+RAIN3) 6159. GO TO 420 6160. C**** Not all remaining rain freezes in layer 4 6161. C**** Some rain could not be frozen, pour it into lake arrays 6162. 250 FREZ4 = - HGI(I,J,4)/ELHM - XGI4*ACE2GI 6163. RAIN4 = RAIN3 - FREZ4 6164. FMGI4 = FMGI3 + FREZ4 6165. FHGI4 = - ELHM*FMGI4 6166. MO(I,J,1) = MO(I,J,1) + FICE*RAIN4*DXYP(J) 6167. BJ(J,21) = BJ(J,21) - FICE*RAIN4*ZATMO(I,J) 6168. CJ(J,21) = CJ(J,21) + FICE*RAIN4*ZATMO(I,J) 6169. BJ(J,33) = BJ(J,33) - FICE*RAIN4 6170. C CJ(J,33) = CJ(J,33) + FICE*RAIN4 ! in DIAGJ 6171. AIJ(I,J,46) = AIJ(I,J,46) + RAIN4 6172. GO TO 420 6173. C**** 6174. C**** All precipitation is snow which increases snow amount 6175. C**** 6176. 300 FHGI1 = HGI(I,J,1)*XGI2*PREC(I,J,1) / (XGI1*MGI(I,J)+PREC(I,J,1)) 6177. MGI(I,J) = MGI(I,J) + PREC(I,J,1) 6178. HGI(I,J,1) = HGI(I,J,1) - FHGI1 6179. HGI(I,J,2) = HGI(I,J,2) + FHGI1 6180. IF(MGI(I,J) .le. SNOMAX+ACE1I) GO TO 500 6181. C**** Too much snow has accumulated, compress some into ice 6182. FMGI4 = MGI(I,J) - (.9d0*SNOMAX+ACE1I) ! > 0 6183. FHGI1 = HGI(I,J,1)*FMGI4 / MGI(I,J) 6184. FHGI2 = HGI(I,J,2)*FMGI4*R2 / MGI(I,J) 6185. MGI(I,J) = .9d0*SNOMAX+ACE1I 6186. C**** 6187. C**** Advect ice (usually downward) 6188. C**** 6189. 400 FHGI3 = HGI(I,J,3)*FMGI4 / (XGI3*ACE2GI) 6190. 410 FHGI4 = HGI(I,J,4)*FMGI4 / (XGI4*ACE2GI) 6191. 420 HGI(I,J,1) = HGI(I,J,1) - FHGI1 6192. HGI(I,J,2) = HGI(I,J,2) + (FHGI1 - FHGI2) 6193. HGI(I,J,3) = HGI(I,J,3) + (FHGI2 - FHGI3) 6194. HGI(I,J,4) = HGI(I,J,4) + (FHGI3 - FHGI4) 6195. BJ(J,36) = BJ(J,36) - FICE*FMGI4 6196. BJ(J,45) = BJ(J,45) - FICE*(FHGI4+GRAV*ZATMO(I,J)*FMGI4) 6197. AIJ(I,J,8) = AIJ(I,J,8) + FMGI4 6198. AIJ(I,J,9) = AIJ(I,J,9) + FHGI4 6199. 500 continue 6200. RETURN 6201. END 7000. 7001. SUBROUTINE GLAICE 7002. C**** 7003. C**** GLAICE receives the atmospheric fluxes of heat and water and 7004. C**** applies them to Glacial Ice, updating the heat, mass and snow 7005. C**** 7006. INCLUDE 'C070.COM' 7007. PARAMETER (RHOI=910., ALAMI=2.1762d0, ALPHA=1., dSNdML=2., 7008. * R1=1/XGI1,R2=1/XGI2,R3=1/XGI3,R4=1/XGI4) 7009. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 7010. * W0(IM,JM,4), E0(IM,JM,4), E1(IM,JM,4) 7011. C**** 7012. C**** GICECB MGI Mass of glacial ice and snow (kg/m^2) 7013. C**** HGI Enthalpy minus latent heat of glacial ice (J/m^2) 7014. C**** 7015. C**** OCENCB MO River and lake mass above sill depth (kg) 7016. C**** 7017. C**** FIXDCB FGICE Glacial ice fraction (1) 7018. C**** 7019. C**** FLUXCB W0(3) Dew received from atmosphere (kg/m^2) 7020. C**** E0(3) Energy received from atmosphere (J/m^2) 7021. C**** E1(3) Energy from layer 1 to layer 2 (J/m**2) 7022. C**** 7023. C**** Outside loop over J and I 7024. C**** 7025. DO 600 J=1,JM 7026. IMAX=IM 7027. IF(J.eq.1 .or. J.eq.JM) IMAX=1 7028. DO 600 I=1,IMAX 7029. IF(FGICE(I,J).le.0.) GO TO 600 7030. FICE = FGICE(I,J) 7031. TGI1 = (HGI(I,J,1)*R1/MGI(I,J) + ELHM) / SHCI 7032. TGI2 = (HGI(I,J,2)*R2/MGI(I,J) + ELHM) / SHCI 7033. TGI3 = (HGI(I,J,3)/(XGI3*ACE2GI) + ELHM) / SHCI 7034. TGI4 = (HGI(I,J,4)/(XGI4*ACE2GI) + ELHM) / SHCI 7035. C**** Accumulate diagnostics 7036. BJ(J,18) = BJ(J,18) + FICE*TGI1 7037. BJ(J,17) = BJ(J,17) + FICE*TGI2 7038. BJ(J,16) = BJ(J,16) + FICE*TGI3 7039. BJ(J,52) = BJ(J,52) + FICE*TGI4 7040. BJ(J,19) = BJ(J,19) + FICE*W0(I,J,3) 7041. BJ(J,22) = BJ(J,22) + FICE*W0(I,J,3)*ZATMO(I,J) 7042. AIJ(I,J,55) = AIJ(I,J,55) - W0(I,J,3) 7043. AIJ(I,J,59) = AIJ(I,J,59) + E0(I,J,3) 7044. IF(MGI(I,J).le.ACE1I) GO TO 10 7045. BJ(J,31) = BJ(J,31) + FICE 7046. BJ(J,53) = BJ(J,53) + FICE*(MGI(I,J)-ACE1I) 7047. AIJ(I,J,2) = AIJ(I,J,2) + FICE 7048. AIJ(I,J,3) = AIJ(I,J,3) + FICE*(MGI(I,J)-ACE1I) 7049. C**** 7050. C**** Calculate and apply diffusive and surface energy fluxes 7051. C**** 7052. 10 HCI2 = SHCI*XGI2*MGI(I,J) 7053. HCI3 = SHCI*XGI3*ACE2GI 7054. dF2dTI = ALAMI*RHOI*DTS / (.5*XGI2*MGI(I,J)+.5*XGI3*ACE2GI) 7055. dF3dTI = ALAMI*RHOI*DTS / (.5*ACE2GI) 7056. CEXP F2 = dF2dTI*(TGI2-TGI3) 7057. CEXP F3 = dF3dTI*(TGI3-TGI4) 7058. F2 = dF2dTI*(HCI2*(TGI2-TGI3) + ALPHA*E1(I,J,3)) / 7059. / (HCI2 + ALPHA*dF2dTI) 7060. F3 = dF3dTI*(HCI3*(TGI3-TGI4) + ALPHA*F2) / (HCI3 + ALPHA*dF3dTI) 7061. HGI(I,J,1) = HGI(I,J,1) + (E0(I,J,3) - E1(I,J,3)) 7062. HGI(I,J,2) = HGI(I,J,2) + (E1(I,J,3) - F2) 7063. HGI(I,J,3) = HGI(I,J,3) + (F2 - F3) 7064. HGI(I,J,4) = HGI(I,J,4) + F3 7065. IF(W0(I,J,3) .gt. 0.) GO TO 400 7066. IF(MGI(I,J)-ACE1I+W0(I,J,3) .ge. 0.) GO TO 300 7067. C**** 7068. C**** All snow and some ice evaporates 7069. C**** 7070. IF(HGI(I,J,1)/ELHM+XGI1*MGI(I,J)+W0(I,J,3) .le. 0.) GO TO 200 7071. C**** Fluxes heat HGI1 to freezing point and melt some ice 7072. MELT1 = HGI(I,J,1)/ELHM + XGI1*MGI(I,J) + W0(I,J,3) 7073. C FMGI1 = XGI1*(MGI-ACE1I) + W0 - MELT1 ! < 0 7074. FHGI1 = - ELHM*(XGI1*(MGI(I,J)-ACE1I)+W0(I,J,3)-MELT1) 7075. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. MELT1) GO TO 130 7076. C**** All melt water is frozen in layer 2 7077. C FREZ2 = MELT1 7078. FMGI4 = MGI(I,J) - ACE1I + W0(I,J,3) ! < 0 7079. GO TO 210 7080. C**** Only some melt water is frozen in layer 2 7081. 130 FREZ2 = - HGI(I,J,2)/ELHM - XGI2*MGI(I,J) 7082. WATR2 = MELT1 - FREZ2 7083. FMGI2 = MGI(I,J) - ACE1I + W0(I,J,3) - WATR2 ! < 0 7084. FHGI2 = - ELHM*FMGI2 7085. MGI(I,J) = ACE1I 7086. IF(-HGI(I,J,3)/ELHM-XGI3*ACE2GI .lt. WATR2) GO TO 150 7087. C**** All remaining melt water freezes in layer 3 7088. C FREZ3 = WATR2 7089. FMGI4 = FMGI2 + WATR2 ! < 0 7090. GO TO 220 7091. C**** Only some remaining melt water freezes in layer 3 7092. 150 FREZ3 = - HGI(I,J,3)/ELHM - XGI3*ACE2GI 7093. WATR3 = WATR2 - FREZ3 7094. FMGI3 = FMGI2 + FREZ3 ! < 0 7095. FHGI3 = - ELHM*FMGI3 7096. IF(-HGI(I,J,4)/ELHM-XGI4*ACE2GI .lt. WATR3) GO TO 460 7097. C**** All remaining melt water freezes in layer 4 7098. FMGI4 = FMGI3 + WATR3 ! < 0 7099. FHGI4 = HGI(I,J,4)*FMGI4 / (XGI4*ACE2GI+WATR3) 7100. GO TO 530 7101. C**** 7102. C**** No ice melts, advect ice upward 7103. C**** 7104. C FMGI1 = XGI1*(MGI-ACE1I) + W0 ! < 0 7105. 200 FMGI4 = MGI(I,J) - ACE1I + W0(I,J,3) ! < 0 7106. FHGI1 = HGI(I,J,2)*(XGI1*(MGI(I,J)-ACE1I)+W0(I,J,3))*R2 / MGI(I,J) 7107. 210 MGI(I,J) = ACE1I 7108. FHGI2 = HGI(I,J,3)*FMGI4 / (XGI3*ACE2GI) 7109. 220 FHGI4 = HGI(I,J,4)*FMGI4 / (XGI4*ACE2GI) 7110. HGI(I,J,1) = HGI(I,J,1) - FHGI1 7111. HGI(I,J,2) = HGI(I,J,2) + (FHGI1 - FHGI2) 7112. HGI(I,J,3) = HGI(I,J,3) + (FHGI2 - FHGI4) 7113. GO TO 550 7114. C**** 7115. C**** Evaporation reduces snow amount 7116. C**** 7117. 300 IF(HGI(I,J,1)/ELHM+XGI1*MGI(I,J)+W0(I,J,3) .le. 0.) GO TO 320 7118. C**** Fluxes heat HGI1 to freezing point and melt some snow or ice 7119. MELT1 = HGI(I,J,1)/ELHM + XGI1*MGI(I,J) + W0(I,J,3) 7120. IF(MGI(I,J)-ACE1I+W0(I,J,3) .le. MELT1) GO TO 420 7121. C**** Melting compresses snow into ice 7122. CMPRS = MIN (dSNdML*MELT1, MGI(I,J)-ACE1I+W0(I,J,3)-MELT1) 7123. C FMGI1 = XGI1*CMPRS + XGI2*(W0(I,J,3) - MELT1) 7124. FHGI1 = - ELHM*(XGI1*CMPRS + XGI2*(W0(I,J,3) - MELT1)) 7125. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. MELT1) GO TO 310 7126. C**** All melt water is frozen in layer 2 7127. C FREZ2 = MELT1 7128. FMGI4 = CMPRS + MELT1 ! > 0 7129. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+MELT1) 7130. MGI(I,J) = MGI(I,J) + W0(I,J,3) - FMGI4 7131. GO TO 510 7132. C**** Only some melt water is frozen in layer 2 7133. 310 FREZ2 = -HGI(I,J,2)/ELHM-XGI2*MGI(I,J) 7134. WATR2 = MELT1 - FREZ2 7135. FMGI2 = CMPRS + FREZ2 ! > 0 7136. FHGI2 = - ELHM*FMGI2 7137. MGI(I,J) = MGI(I,J) + W0(I,J,3) - MELT1 - CMPRS 7138. GO TO 440 7139. C**** No ice melts, advect ice upward into layer 1 from layer 2 7140. C FMGI1 = XGI2*W0(I,J,3) ! < 0 7141. 320 FHGI1 = HGI(I,J,2)*W0(I,J,3) / MGI(I,J) 7142. MGI(I,J) = MGI(I,J) + W0(I,J,3) 7143. HGI(I,J,1) = HGI(I,J,1) - FHGI1 7144. HGI(I,J,2) = HGI(I,J,2) + FHGI1 7145. GO TO 600 7146. C**** 7147. C**** Dew increases ice amount 7148. C**** 7149. 400 IF(HGI(I,J,1)/ELHM+XGI1*MGI(I,J)+W0(I,J,3) .le. 0.) GO TO 500 7150. C**** Fluxes heat HGI1 to freezing point and melt some snow or ice 7151. MELT1 = HGI(I,J,1)/ELHM + XGI1*MGI(I,J) + W0(I,J,3) 7152. IF(MGI(I,J)-ACE1I .le. MELT1) GO TO 420 7153. C**** Melting compresses snow into ice 7154. CMPRS = MIN (dSNdML*MELT1, MGI(I,J)-ACE1I-MELT1) 7155. C FMGI1 = XGI1*CMPRS - XGI2*MELT1 + W0(I,J,3) 7156. FHGI1 = - ELHM*(XGI1*CMPRS - XGI2*MELT1 + W0(I,J,3)) 7157. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. MELT1) GO TO 410 7158. C**** All melt water is frozen in layer 2 7159. C FREZ2 = MELT1 7160. FMGI4 = MELT1 + CMPRS + W0(I,J,3) ! > 0 7161. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+MELT1) 7162. MGI(I,J) = MGI(I,J) - MELT1 - CMPRS 7163. GO TO 510 7164. C**** Only some melt water is frozen in layer 2 7165. 410 FREZ2 = -HGI(I,J,2)/ELHM-XGI2*MGI(I,J) 7166. WATR2 = MELT1 - FREZ2 7167. FMGI2 = CMPRS + FREZ2 + W0(I,J,3) ! > 0 7168. FHGI2 = - ELHM*FMGI2 7169. MGI(I,J) = MGI(I,J) - MELT1 - CMPRS 7170. GO TO 440 7171. C**** All snow and some ice has melted, what can be frozen in layer 2? 7172. C FMGI1 = XGI1*(MGI-ACE1I) + W0(I,J,3) - MELT1 7173. 420 FHGI1 = - ELHM*(XGI1*(MGI(I,J)-ACE1I)+W0(I,J,3)-MELT1) 7174. IF(-HGI(I,J,2)/ELHM-XGI2*MGI(I,J) .lt. MELT1) GO TO 430 7175. C**** All melt water is frozen in layer 2 7176. C FREZ2 = MELT1 7177. FMGI4 = MGI(I,J) - ACE1I + W0(I,J,3) ! > 0 7178. FHGI2 = HGI(I,J,2)*FMGI4 / (XGI2*MGI(I,J)+MELT1) 7179. MGI(I,J) = ACE1I 7180. GO TO 510 7181. C**** Only some melt water is frozen in layer 2 7182. 430 FREZ2 = - HGI(I,J,2)/ELHM - XGI2*MGI(I,J) 7183. WATR2 = MELT1 - FREZ2 7184. FMGI2 = MGI(I,J) - ACE1I + W0(I,J,3) - WATR2 7185. FHGI2 = - ELHM*FMGI2 7186. MGI(I,J) = ACE1I 7187. 440 IF(-HGI(I,J,3)/ELHM-XGI3*ACE2GI .lt. WATR2) GO TO 450 7188. C**** All remaining melt water freezes in layer 3 7189. C FREZ3 = WATR2 7190. FMGI4 = FMGI2 + WATR2 ! > 0 7191. FHGI3 = HGI(I,J,3)*FMGI4 / (XGI3*ACE2GI+WATR2) 7192. GO TO 520 7193. C**** Only some remaining melt water freezes in layer 3 7194. 450 FREZ3 = - HGI(I,J,3)/ELHM - XGI3*ACE2GI 7195. WATR3 = WATR2 - FREZ3 7196. FMGI3 = FMGI2 + FREZ3 ! > 0 7197. FHGI3 = - ELHM*FMGI3 7198. IF(-HGI(I,J,4)/ELHM-XGI4*ACE2GI .lt. WATR3) GO TO 460 7199. C**** All remaining melt water freezes in layer 4 7200. FMGI4 = FMGI3 + WATR3 ! > 0 7201. FHGI4 = HGI(I,J,4)*FMGI4 / (XGI4*ACE2GI+WATR3) 7202. GO TO 530 7203. C**** Only some remaining melt water freezes in layer 4 7204. C**** Some melt water could not be frozen, pour it into lake arrays 7205. 460 FREZ4 = - HGI(I,J,4)/ELHM - XGI4*ACE2GI 7206. WATR4 = WATR3 - FREZ4 7207. FMGI4 = FMGI3 + FREZ4 ! > 0 7208. FHGI4 = - ELHM*FMGI4 7209. MO(I,J,1) = MO(I,J,1) + FICE*WATR4*DXYP(J) 7210. BJ(J,22) = BJ(J,22) - FICE*WATR4*ZATMO(I,J) 7211. CJ(J,22) = CJ(J,22) + FICE*WATR4*ZATMO(I,J) 7212. BJ(J,38) = BJ(J,38) - FICE*WATR4 7213. C CJ(J,38) = CJ(J,38) + FICE*WATR4 ! in DIAGJ 7214. AIJ(I,J,46) = AIJ(I,J,46) + WATR4 7215. GO TO 530 7216. C**** 7217. C**** Advect ice (usually downward) 7218. C**** 7219. 500 FMGI4 = W0(I,J,3) ! > 0 7220. FHGI1 = HGI(I,J,1)*FMGI4 / (XGI1*MGI(I,J)+W0(I,J,3)) 7221. FHGI2 = HGI(I,J,2)*FMGI4*R2 / MGI(I,J) 7222. 510 FHGI3 = HGI(I,J,3)*FMGI4 / (XGI3*ACE2GI) 7223. 520 FHGI4 = HGI(I,J,4)*FMGI4 / (XGI4*ACE2GI) 7224. 530 HGI(I,J,1) = HGI(I,J,1) - FHGI1 7225. HGI(I,J,2) = HGI(I,J,2) + (FHGI1 - FHGI2) 7226. HGI(I,J,3) = HGI(I,J,3) + (FHGI2 - FHGI3) 7227. HGI(I,J,4) = HGI(I,J,4) + (FHGI3 - FHGI4) 7228. 550 BJ(J,37) = BJ(J,37) - FICE* FMGI4 7229. BJ(J,46) = BJ(J,46) - FICE*(FHGI4+GRAV*ZATMO(I,J)*FMGI4) 7230. AIJ(I,J,8) = AIJ(I,J,8) + FMGI4 7231. AIJ(I,J,9) = AIJ(I,J,9) + FHGI4 7232. 600 continue 7233. RETURN 7234. END