1. C**** 2. C**** C077hI.S Fortran source code for Sea Ice routines 2003/11/19 3. C**** 4. C**** NASA/GISS Climate Model III Gary L. Russell 5. C**** Gavin A. Schmidt 6. C**** NASA/Goddard Space Flight Center 7. C**** Institute for Space Studies 8. C**** 2880 Broadway, New York, NY 10025 9. C**** U.S.A. 10. C**** 11. SUBROUTINE DYNSI 12. C**** 13. C**** DYNSI receives the surface momentum stress and current sea ice 14. C**** velocity, and calculates the time integral of the velocity and 15. C**** the sea ice velocity at the end of the time step. Horizontal 16. C**** sea ice stresses are the gradient of internal sea ice pressure. 17. C**** 18. C**** M (kg/m^2) = MSIL = sea ice mass per unit area 19. C**** u,v (m/s) = USI,VSI = sea ice velocity components 20. C**** x,y (m/s) = UO,VO = ocean velocity components of layer 1 21. C**** p,q (m/s^2) = PGFU,PGFV = pressure gradient force 22. C**** g,h (m/s^2) = internal sea ice stress components 23. C**** BM (1/s) = coast line blocking factor * RSI * MSI2 24. C**** r,s (N/m^2) = DMUA/DTS,DMVA/DTS = surface stress components 25. C**** C/a (kg/m^3)= 1 + .0015*M = ocean drag times density 26. C**** 27. C**** du/dt = r/M - (u-x)|(u,v)-(x,y)|C/aM - uBM + (f+uùtaní/R)v - p - g 28. C**** dv/dt = s/M - (v-y)|(u,v)-(x,y)|C/aM - vBM - (f+uùtaní/R)u - q - h 29. C**** 30. C**** The above equations are rewritten as: 31. C**** 32. C**** du/dt = au + bv + c where a = -|(u,v)-(x,y)|C/aM - BM 33. C**** dv/dt = av - bu + d b = f + uùtaní/R 34. C**** c = r/M + x|(u,v)-(x,y)|C/aM - p - g 35. C**** d = s/M + y|(u,v)-(x,y)|C/aM - q - h 36. C**** 37. C**** The general solution to these equations is: 38. C**** 39. C**** u = X exp(At) sin(Bt) + Y exp(At) cos(Bt) + Z 40. C**** v = U exp(At) sin(Bt) + V exp(At) cos(Bt) + W 41. C**** 42. C**** and the particular solution is: 43. C**** 44. C**** Z = (-ca+db)/(a^2+b^2) W = (-cb-da)/(a^2+b^2) 45. C**** Y = u(t=0) - Z V = v(t=0) - W 46. C**** A = a B = b X = V U = -Y 47. C**** 48. INCLUDE 'C070.COM' 49. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 50. * W0(IM,JM,4), E0(IM,JM,4),E1(IM,JM,4), UAS(IM,JM),VAS(IM,JM), 51. * DMUI(IM,JM),DMVI(IM,JM), USIDT(IM,JM),VSIDT(IM,JM) 52. C**** 53. C**** SURFCB BLDATA(6) pressure at bottom of sea ice (Pa-101325) 54. C**** BLDATA(8) geopotential above water surface (m^2/s^2) 55. C**** 56. C**** FLUXCB DMUA(2) U momentum downward into sea ice (kg/m*s) 57. C**** DMVA(2) V momentum downward into sea ice (kg/m*s) 58. C**** DMUI U momentum downward from sea ice (kg/m*s) 59. C**** DMVI V momentum downward from sea ice (kg/m*s) 60. C**** USIDT U compon of time integrated sea ice velocity (m) 61. C**** VSIDT V compon of time integrated sea ice velocity (m) 62. C**** 63. PARAMETER (NSILIM=20, BLOCKX=.1d0/3600., ERRLIM=1., ACE2SI=182., 64. * RHOI=910., OMEGA=TWOPI*366./(365.*SDAY)) 65. REAL*8 BYDXP(JM),DYBYDX(JM),SINP2O(JM),TANPBR(JM), 66. * BYDYV(JM),DXBYDY(JM),SINV2O(JM),TANVBR(JM) 66.1 REAL*4 BLOCK(IM,JM) 66.2 CHARACTER*80 TITLE 67. LOGICAL*4 QSIM,QSIU,QSIV 68. COMMON /WORK01/ QSIM(IM,JM),QSIU(IM,JM),QSIV(IM,JM), 69. * PSIMAX(IM,JM),MSIL(IM,JM),MMSI(IM,JM),BRMU(IM,JM),BRMV(IM,JM), 70. * USIAV(IM,JM),DPSIU(IM,JM),PGFU(IM,JM),RSIU(IM,JM),ASIEW(IM,JM), 71. * VSIAV(IM,JM),DPSIV(IM,JM),PGFV(IM,JM),RSIV(IM,JM),ASINS(IM,JM), 72. * BYMMIU(IM,JM),ABYABU(IM,JM),XSNBTU(IM,JM),ZU(IM,JM), 73. * BYMMIV(IM,JM),ABYABV(IM,JM),XSNBTV(IM,JM),ZV(IM,JM), 74. * RCDDWU(IM,JM),BBYABU(IM,JM),XCSBTU(IM,JM),WU(IM,JM), 75. * RCDDWV(IM,JM),BBYABV(IM,JM),XCSBTV(IM,JM),WV(IM,JM) 76. C**** 77. DATA IFIRST /1/ 78. IF(IFIRST.eq.0) GO TO 100 79. IFIRST = 0 80. C**** 81. C**** Eliminate division by constant quantities 82. C**** 83. BYDTS = 1./DTS 84. BYDYP = 1./DYP(3) 85. DLAT = NINT(180./(JM-1))*TWOPI/360. 86. DO 10 J=2,JM-1 87. BYDXP(J) = 1./DXP(J) 88. DYBYDX(J) = DYP(J)/DXP(J) 89. SINP2O(J) = SINP(J)*2.*OMEGA 90. 10 TANPBR(J) = SINP(J)/(COSP(J)*RADIUS) 91. DO 20 J=1,JM-1 92. BYDYV(J) = 1./DYV(J) 93. DXBYDY(J) = DXV(J)/DYV(J) 94. SINV = SIN(DLAT*(J-JM/2)) 95. SINV2O(J) = SINV*2.*OMEGA 96. 20 TANVBR(J) = SINV/(COSV(J)*RADIUS) 98. C**** Read in island and coast line blocking factor 100. READ (33) TITLE,BLOCK 101. WRITE (6,*) 'Read on unit 25: ',TITLE 102. CLOSE (33) 124. C**** 125. C**** Quantities defined on A grid and are fixed for this call of DYNSI: 126. C**** 127. C**** MSIL (kg/m^2)= vertically integrated sea ice mass per unit area 128. C**** MMSI (kg) = vertically integrated sea ice mass 129. C**** PSIMAX (Pa) = maximum internal sea ice pressure 130. C**** 131. 100 DO 110 J=2,JM-1 132. DO 110 I=1,IM 133. MSIL(I,J) = 0. 134. MMSI(I,J) = 0. 135. QSIM(I,J) = LMM(I,J).gt.0. .and. RSI(I,J).gt.0. 136. IF(.not.QSIM(I,J)) GO TO 110 137. MSIL(I,J) = MSI(I,J,1)+MSI(I,J,2) 138. MMSI(I,J) = (MSI(I,J,1)+MSI(I,J,2))*RSI(I,J)*DXYP(J) 139. PSIMAX(I,J) = (27500./RHOI)*MSIL(I,J)*EXP(-20.+20.*RSI(I,J)) 141. AIJ(I,J,16) = AIJ(I,J,16) + RSI(I,J)*MSIL(I,J) 142. AIJ(I,J,38) = AIJ(I,J,38) + RSI(I,J)* PSI(I,J) 143. 110 continue 144. C**** Define or replicate values at the North Pole 145. MSIL(1,JM) = 0. 146. MMSI(1,JM) = 0. 147. PGFU(1,JM) = 0. 148. DPSIU(1,JM) = 0. 149. QSIM(1,JM) = LMM(1,JM).gt.0. .and. RSI(1,JM).gt.0. 150. C IF(.not.QSIM(1,JM)) GO TO 200 151. MSIL(1,JM) = MSI(1,JM,1)+MSI(1,JM,2) 152. MMSI(1,JM) = (MSI(1,JM,1)+MSI(1,JM,2))*RSI(1,JM)*DXYP(JM) 153. PSIMAX(1,JM)= (27500./RHOI)*MSIL(1,JM)*EXP(-20.+20.*RSI(1,JM)) 155. AIJ(1,JM,16) = AIJ(1,JM,16) + RSI(1,JM)*MSIL(1,JM) 156. AIJ(1,JM,38) = AIJ(1,JM,38) + RSI(1,JM)* PSI(1,JM) 157. DO 120 I=2,IM 158. RSI(I,JM) = RSI(1,JM) 159. DMUA(I,JM,2) = DMUA(1,JM,2)*COSI(I) + DMVA(1,JM,2)*SINI(I) 160. DMVA(I,JM,2) = DMVA(1,JM,2)*COSI(I) - DMUA(1,JM,2)*SINI(I) 161. MSIL(I,JM) = MSIL(1,JM) 162. MMSI(I,JM) = MMSI(1,JM) 163. PGFU(I,JM) = 0. 164. 120 DPSIU(I,JM) = 0. 165. C**** 166. C**** Quantities defined on C grid and are fixed for this call of DYNSI: 167. C**** 168. C**** BYMMIU (1/kg) = reciprocal of sea ice mass on C-grid points 168.5 C**** BRMU (1/s) = island and coast line blocking factor 169. C**** PGFU (m/s^2) = deceleration due to pressure gradient force 170. C**** DPSIU (m/s^2) = deceleration due to internal sea ice pressure 171. C**** 172. C**** U component 173. 200 I=IM 174. DO 220 J=2,JM-1 175. DO 220 IP1=1,IM 176. QSIU(I,J) = LMU(I,J).gt.0. .and. RSI(I,J)+RSI(IP1,J).gt.0. 177. IF(QSIU(I,J)) GO TO 210 178. USI(I,J) = 0. 179. USIDT(I,J) = 0. 179.5 BRMU(I,J) = 0. 180. PGFU(I,J) = 0. 181. RSIU(I,J) = 0. 182. ASIEW(I,J) = 0. 183. DPSIU(I,J) = 0. 184. GO TO 220 185. 210 BYMMIU(I,J)= 2./(MMSI(I,J)+MMSI(IP1,J)) 185.5 BRMU(I,J) = (BLOCK(I,J)*RSI(I,J)*MSI(I,J,2)+ 185.6 * BLOCK(IP1,J)*RSI(IP1,J)*MSI(IP1,J,2))*(BLOCKX*.5/ACE2SI) 186. PGFU(I,J) = BYDXP(J)*((BLDATA(IP1,J,6)-BLDATA(I,J,6))/RHOI + 187. * BLDATA(IP1,J,8)-BLDATA(I,J,8)) 188. DPSIU(I,J) = BYDXP(J)*(PSI(IP1,J)-PSI(I,J))/RHOI 189. AIJ(I,J,67) = AIJ(I,J,67) + (RSI(I,J)+RSI(IP1,J))*USI(I,J) 190. 220 I=IP1 191. C**** V component 192. DO 240 J=1,JM-1 193. DO 240 I=1,IM 194. QSIV(I,J) = LMV(I,J).gt.0. .and. RSI(I,J)+RSI(I,J+1).gt.0. 195. IF(QSIV(I,J)) GO TO 230 196. VSI(I,J) = 0. 197. VSIDT(I,J) = 0. 197.5 BRMV(I,J) = 0. 198. PGFV(I,J) = 0. 199. RSIV(I,J) = 0. 200. ASINS(I,J) = 0. 201. DPSIV(I,J) = 0. 202. GO TO 240 203. 230 BYMMIV(I,J)= 2./(MMSI(I,J)+MMSI(I,J+1)) 203.5 BRMV(I,J) = (BLOCK(I,J)*RSI(I,J)*MSI(I,J,2)+ 203.6 * BLOCK(I,J+1)*RSI(I,J+1)*MSI(I,J+1,2))*(BLOCKX*.5/ACE2SI) 204. PGFV(I,J) = BYDYV(J)*((BLDATA(I,J+1,6)-BLDATA(I,J,6))/RHOI + 205. * BLDATA(I,J+1,8)-BLDATA(I,J,8)) 206. DPSIV(I,J) = BYDYV(J)*(PSI(I,J+1)-PSI(I,J))/RHOI 207. AIJ(I,J,68) = AIJ(I,J,68) + (RSI(I,J)+RSI(I,J+1))*VSI(I,J) 208. 240 continue 209. C**** 210. C**** Quantities defined on C grid, are fixed for this call of DYNSI, 211. C**** and require four point averaging of the other component: 212. C**** 213. C**** VSIAV (m/s) = four point average of sea ice velocity component 214. C**** RCDDWU (kg/s*m^2)= magnitude of sea ice minus ocean velocity * C/a 215. C**** ABYABU (s) = a/(a^2+b^2) 216. C**** BBYABU (s) = b/(a^2+b^2) 217. C**** ZU (m/s) = (-ca+db)/(a^2+b^2) excluding sea ice pressure grad 218. C**** WU (m/s) = (-cb-da)/(a^2+b^2) excluding sea ice pressure grad 219. C**** XSNBTU = exp(a*DTS) sin(b*DTS) 220. C**** XCSBTU = exp(a*DTS) cos(b*DTS) 221. C**** 222. C**** U component 223. I=IM 224. DO 260 J=2,JM-1 225. DO 260 IP1=1,IM 226. IF(.not.QSIU(I,J)) GO TO 260 227. VSIAV(I,J) = .25*((VSI(I,J-1)+VSI(I ,J))*MMSI(I ,J) + 228. * (VSI(IP1,J-1)+VSI(IP1,J))*MMSI(IP1,J))*BYMMIU(I,J) 229. MSILU = .5*(MSIL(I,J)+MSIL(IP1,J)) 230. VOAV = .25*(VO(I,J-1,1)+VO(I,J,1)+VO(IP1,J-1,1)+VO(IP1,J,1)) 231. PGFVAV= .25*((PGFV(I ,J-1)+PGFV(I ,J))*MMSI(I ,J) + 232. * (PGFV(IP1,J-1)+PGFV(IP1,J))*MMSI(IP1,J))*BYMMIU(I,J) 233. DW = SQRT((USI(I,J)-UO(I,J,1))**2 + (VSIAV(I,J)-VOAV)**2) 234. RCDDWU(I,J) = (1.+.0015d0*MSILU)*DW 235. A = -(RCDDWU(I,J)/MSILU) - BRMU(I,J) 236. B = SINP2O(J) + USI(I,J)*TANPBR(J) 237. C = .5*(DMUA(I ,J,2)*RSI(I ,J)+ 238. * DMUA(IP1,J,2)*RSI(IP1,J))*(BYDTS*BYMMIU(I,J)*DXYP(J)) 239. * + UO(I,J,1)*(RCDDWU(I,J)/MSILU) - PGFU(I,J) 240. D = .5*(DMVA(I ,J,2)*RSI(I ,J)+ 241. * DMVA(IP1,J,2)*RSI(IP1,J))*(BYDTS*BYMMIU(I,J)*DXYP(J)) 242. * + VOAV*(RCDDWU(I,J)/MSILU) - PGFVAV 243. ABYABU(I,J) = A / (A*A + B*B) 244. BBYABU(I,J) = B / (A*A + B*B) 245. ZU(I,J) = - C*ABYABU(I,J) + D*BBYABU(I,J) 246. WU(I,J) = - C*BBYABU(I,J) - D*ABYABU(I,J) 247. EXPAT = EXP(A*DTS) 248. XSNBTU(I,J) = EXPAT*SIN(B*DTS) 249. XCSBTU(I,J) = EXPAT*COS(B*DTS) 250. 260 I=IP1 251. C**** V component 252. IM1=IM 253. DO 280 J=1,JM-1 254. DO 280 I=1,IM 255. IF(.not.QSIV(I,J)) GO TO 280 256. USIAV(I,J) = .25*((USI(IM1,J)+USI(I,J ))*MMSI(I,J ) + 257. * (USI(IM1,J+1)+USI(I,J+1))*MMSI(I,J+1))*BYMMIV(I,J) 258. MSILV = .5*(MSIL(I,J)+MSIL(I,J+1)) 259. UOAV = .25*(UO(IM1,J,1)+UO(I,J,1)+UO(IM1,J+1,1)+UO(I,J+1,1)) 260. PGFUAV= .25*((PGFU(IM1,J )+PGFU(I,J ))*MMSI(I,J ) + 261. * (PGFU(IM1,J+1)+PGFU(I,J+1))*MMSI(I,J+1))*BYMMIV(I,J) 262. DW = SQRT((USIAV(I,J)-UOAV)**2 + (VSI(I,J)-VO(I,J,1))**2) 263. RCDDWV(I,J) = (1.+.0015d0*MSILV)*DW 264. A = -(RCDDWV(I,J)/MSILV) - BRMV(I,J) 265. B = SINV2O(J) + USIAV(I,J)*TANVBR(J) 266. C = .5*(DMUA(I,J ,2)*RSI(I,J )*DXYP(J )+ 267. * DMUA(I,J+1,2)*RSI(I,J+1)*DXYP(J+1))*(BYDTS*BYMMIV(I,J)) 268. * + UOAV*(RCDDWV(I,J)/MSILV) - PGFUAV 269. D = .5*(DMVA(I,J ,2)*RSI(I,J )*DXYP(J )+ 270. * DMVA(I,J+1,2)*RSI(I,J+1)*DXYP(J+1))*(BYDTS*BYMMIV(I,J)) 271. * + VO(I,J,1)*(RCDDWV(I,J)/MSILV) - PGFV(I,J) 272. ABYABV(I,J) = A / (A*A + B*B) 273. BBYABV(I,J) = B / (A*A + B*B) 274. ZV(I,J) = - C*ABYABV(I,J) + D*BBYABV(I,J) 275. WV(I,J) = - C*BBYABV(I,J) - D*ABYABV(I,J) 276. EXPAT = EXP(A*DTS) 277. XSNBTV(I,J) = EXPAT*SIN(B*DTS) 278. XCSBTV(I,J) = EXPAT*COS(B*DTS) 279. 280 IM1=I 280. C**** 281. C**** Iteratively solve for PSI so that sea ice area entering a grid 282. C**** box converges to zero or is divergent 283. C**** 284. NSI=0 285. 300 NSI=NSI+1 286. C**** 287. C**** Calculate time integral of sea ice velocity (m) 288. C**** 289. C**** U component 290. I=IM 291. DO 320 J=2,JM-1 292. DO 320 IP1=1,IM 293. IF(.not.QSIU(I,J)) GO TO 320 294. DPSIVA= .25*((DPSIV(I,J-1)+DPSIV(I ,J))*MMSI(I ,J) + 295. * (DPSIV(IP1,J-1)+DPSIV(IP1,J))*MMSI(IP1,J))*BYMMIU(I,J) 296. ZNEW = ZU(I,J) + DPSIU(I,J)*ABYABU(I,J) - DPSIVA*BBYABU(I,J) 297. WNEW = WU(I,J) + DPSIU(I,J)*BBYABU(I,J) + DPSIVA*ABYABU(I,J) 298. Y = USI(I,J) - ZNEW 299. X = VSIAV(I,J) - WNEW 300. USIDT(I,J) = XSNBTU(I,J) *(X*ABYABU(I,J)+Y*BBYABU(I,J)) 301. * + (XCSBTU(I,J)-1.)*(Y*ABYABU(I,J)-X*BBYABU(I,J)) 302. * + ZNEW*DTS 303. 320 I=IP1 304. C**** V component 305. IM1=IM 306. DO 340 J=1,JM-1 307. DO 340 I=1,IM 308. IF(.not.QSIV(I,J)) GO TO 340 309. DPSIUA= .25*((DPSIU(IM1,J)+DPSIU(I,J ))*MMSI(I,J ) + 310. * (DPSIU(IM1,J+1)+DPSIU(I,J+1))*MMSI(I,J+1))*BYMMIV(I,J) 311. ZNEW = ZV(I,J) + DPSIUA*ABYABV(I,J) - DPSIV(I,J)*BBYABV(I,J) 312. WNEW = WV(I,J) + DPSIUA*BBYABV(I,J) + DPSIV(I,J)*ABYABV(I,J) 313. Y = USIAV(I,J) - ZNEW 314. X = VSI(I,J) - WNEW 315. VSIDT(I,J) = XSNBTV(I,J) *(X*BBYABV(I,J)-Y*ABYABV(I,J)) 316. * + (XCSBTV(I,J)-1.)*(Y*BBYABV(I,J)+X*ABYABV(I,J)) 317. * + WNEW*DTS 318. 340 IM1=I 319. C**** 320. C**** Calculate sea ice area fluxes (m^2) 321. C**** 322. C**** U component 323. I=IM 324. DO 420 J=2,JM-1 325. DO 420 IP1=1,IM 326. IF(.not.QSIU(I,J)) GO TO 420 327. IF(USIDT(I,J).ge.0.) GO TO 410 328. RSIU(I,J) = RSI(IP1,J) - 329. * (1.+USIDT(I,J)*DYP(J)*BYDXYP(J))*RSIX(IP1,J) 330. ASIEW(I,J) = USIDT(I,J)*DYP(J)*RSIU(I,J) 331. GO TO 420 332. 410 RSIU(I,J) = RSI(I,J) + 333. * (1.-USIDT(I,J)*DYP(J)*BYDXYP(J))*RSIX(I,J) 334. ASIEW(I,J) = USIDT(I,J)*DYP(J)*RSIU(I,J) 335. 420 I=IP1 336. C**** V component 337. DO 440 J=2,JM-1 338. DO 440 I=1,IM 339. IF(.not.QSIV(I,J)) GO TO 440 340. IF(VSIDT(I,J).ge.0.) GO TO 430 341. RSIV(I,J) = RSI(I,J+1) - 342. * (1.+VSIDT(I,J)*DXV(J)*BYDXYP(J+1))*RSIY(I,J+1) 343. ASINS(I,J) = VSIDT(I,J)*DXV(J)*RSIV(I,J) 344. GO TO 440 345. 430 RSIV(I,J) = RSI(I,J) + 346. * (1.-VSIDT(I,J)*DXV(J)*BYDXYP(J))*RSIY(I,J) 347. ASINS(I,J) = VSIDT(I,J)*DXV(J)*RSIV(I,J) 348. 440 continue 349. C**** 350. C**** Update sea ice pressure: PSI (Pa) 351. C**** 352. ERRMAX = 0. 353. IM1=IM 354. DO 540 J=2,JM-1 355. DO 540 I=1,IM 356. IF(.not.QSIM(I,J)) GO TO 540 357. CONVRG = ASIEW(IM1,J) - ASIEW(I,J) + ASINS(I,J-1) - ASINS(I,J) 358. IF(CONVRG) 510,540,520 359. C**** Sea ice is diverging from grid box 360. 510 IF(PSI(I,J).le.0.) GO TO 540 361. PSINEW = PSI(I,J) + 2.*RHOI*BYDTS*BYDTS*CONVRG / 362. * ((RSIU(IM1,J)+RSIU(I,J))*DYBYDX(J) 363. * + RSIV(I,J-1)*DXBYDY(J-1)+RSIV(I,J)*DXBYDY(J)) 364. IF(PSINEW.lt.0) PSINEW = 0. 365. GO TO 530 366. C**** Sea ice is converging into grid box 367. 520 IF(PSI(I,J).ge.PSIMAX(I,J)) GO TO 540 368. PSINEW = PSI(I,J) + 2.*RHOI*BYDTS*BYDTS*CONVRG / 369. * ((RSIU(IM1,J)+RSIU(I,J))*DYBYDX(J) 370. * + RSIV(I,J-1)*DXBYDY(J-1)+RSIV(I,J)*DXBYDY(J)) 371. IF(PSINEW.gt.PSIMAX(I,J)) PSINEW = PSIMAX(I,J) 372. C 530 IF(ABS(PSINEW-PSI(I,J)).gt.ERRMAX) ERRMAX = ABS(PSINEW-PSI(I,J)) 373. 530 IF(ABS(PSINEW-PSI(I,J)).le.ERRMAX) GO TO 535 374. ERRMAX = ABS(PSINEW-PSI(I,J)) 375. IMAX=I 376. JMAX=J 377. 535 continue 378. PSI(I,J) = PSINEW 379. 540 IM1=I 380. C**** 381. C**** Update sea ice pressure at North Pole 382. C**** 383. IF(.not.QSIM(1,JM)) GO TO 600 384. SRSIV = 0. 385. SCONV = 0. 386. DO 550 I=1,IM 387. SRSIV = SRSIV + RSIV(I,JM-1) 388. 550 SCONV = SCONV + ASINS(I,JM-1) 389. IF(SCONV) 560,600,570 390. C**** Sea ice is diverging from North Pole box 391. 560 IF(PSI(1,JM).le.0.) GO TO 600 392. PSINEW = PSI(1,JM) + 2.*RHOI*BYDTS*BYDTS*SCONV / 393. * (SRSIV*DXBYDY(JM-1)) 394. IF(PSINEW.lt.0.) PSINEW = 0. 395. GO TO 580 396. C**** Sea ice is converging into North Pole box 397. 570 IF(PSI(1,JM).ge.PSIMAX(1,JM)) GO TO 600 398. PSINEW = PSI(1,JM) + 2.*RHOI*BYDTS*BYDTS*SCONV / 399. * (SRSIV*DXBYDY(JM-1)) 400. IF(PSINEW.gt.PSIMAX(1,JM)) PSINEW = PSIMAX(1,JM) 401. C**** Replicate PSI(1,JM) to all longitudes 402. 580 DO 590 I=1,IM 403. 590 PSI(I,JM) = PSINEW 404. C**** 405. C**** Calculate internal sea ice stress on C grid: DPSIU (m/s^2) 406. C**** 407. 600 I=IM 408. DO 610 J=2,JM-1 409. DO 610 IP1=1,IM 410. C IF(.not.QSIU(I,J) .and. .not.QSIV(I,J)) GO TO 610 411. DPSIU(I,J) = (PSI(IP1,J)-PSI(I,J))*BYDXP(J)/RHOI 412. DPSIV(I,J) = (PSI(I,J+1)-PSI(I,J))*BYDYV(J)/RHOI 413. 610 I=IP1 414. C**** 415. IF(NSI.lt.NSILIM .and. ERRMAX.gt.ERRLIM) GO TO 300 415.5 C WRITE (6,*) 'NSI,IMAX,JMAX,ERRMAX=',NSI,IMAX,JMAX,ERRMAX 416. C**** 417. C**** Final iteration for new sea ice velocities. 418. C**** Calculate sea ice velocity at end of time step, 419. C**** calculate sea ice - ocean stress, 420. C**** accumulate sea ice diagnostics 421. C**** 422. C**** U component 423. I=IM 424. DO 710 J=2,JM-1 425. DO 710 IP1=1,IM 426. IF(.not.QSIU(I,J)) GO TO 710 427. DPSIVA= .25*((DPSIV(I,J-1)+DPSIV(I ,J))*MMSI(I ,J) + 428. * (DPSIV(IP1,J-1)+DPSIV(IP1,J))*MMSI(IP1,J))*BYMMIU(I,J) 429. ZNEW = ZU(I,J) + DPSIU(I,J)*ABYABU(I,J) - DPSIVA*BBYABU(I,J) 430. WNEW = WU(I,J) + DPSIU(I,J)*BBYABU(I,J) + DPSIVA*ABYABU(I,J) 431. Y = USI(I,J) - ZNEW 432. X = VSIAV(I,J) - WNEW 433. USIDT(I,J) = XSNBTU(I,J) *(X*ABYABU(I,J)+Y*BBYABU(I,J)) 434. * + (XCSBTU(I,J)-1.)*(Y*ABYABU(I,J)-X*BBYABU(I,J)) 435. * + ZNEW*DTS 436. USI(I,J) = X*XSNBTU(I,J) + Y*XCSBTU(I,J) + ZNEW 437. DMUI(I,J) = (USIDT(I,J)-UO(I,J,1)*DTS)*RCDDWU(I,J) 438. AIJ(I,J,69) = AIJ(I,J,69) + (RSI(I,J)+RSI(IP1,J))*DMUI(I,J) 439. 710 I=IP1 440. C**** V component 441. IM1=IM 442. DO 720 J=2,JM-1 443. DO 720 I=1,IM 444. IF(.not.QSIV(I,J)) GO TO 720 445. DPSIUA= .25*((DPSIU(IM1,J)+DPSIU(I,J ))*MMSI(I,J ) + 446. * (DPSIU(IM1,J+1)+DPSIU(I,J+1))*MMSI(I,J+1))*BYMMIV(I,J) 447. ZNEW = ZV(I,J) + DPSIUA*ABYABV(I,J) - DPSIV(I,J)*BBYABV(I,J) 448. WNEW = WV(I,J) + DPSIUA*BBYABV(I,J) + DPSIV(I,J)*ABYABV(I,J) 449. Y = USIAV(I,J) - ZNEW 450. X = VSI(I,J) - WNEW 451. VSIDT(I,J) = XSNBTV(I,J) *(X*BBYABV(I,J)-Y*ABYABV(I,J)) 452. * + (XCSBTV(I,J)-1.)*(Y*BBYABV(I,J)+X*ABYABV(I,J)) 453. * + WNEW*DTS 454. VSI(I,J) = -Y*XSNBTV(I,J) + X*XCSBTV(I,J) + WNEW 455. DMVI(I,J) = (VSIDT(I,J)-VO(I,J,1)*DTS)*RCDDWV(I,J) 456. AIJ(I,J,70) = AIJ(I,J,70) + (RSI(I,J)+RSI(I,J+1))*DMVI(I,J) 457. 720 IM1=I 458. RETURN 459. END 1000. 1001. SUBROUTINE ADVSI 1002. C**** 1003. C**** ADVSI uses the time integrated velocities to advect sea ice 1004. C**** from one grid box to another and to change the north-south 1005. C**** center of sea ice mass within a grid box. 1006. C**** 1007. INCLUDE 'C070.COM' 1008. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 1009. * W0(IM,JM,4), E0(IM,JM,4),E1(IM,JM,4), UAS(IM,JM),VAS(IM,JM), 1010. * DMUI(IM,JM),DMVI(IM,JM), USIDT(IM,JM),VSIDT(IM,JM) 1011. COMMON /WORK01/ FAW(IM),FASI(IM),FXSI(IM),FYSI(IM), 1012. * FMSI(2+LMSI,IM),SFMSI(2+LMSI),AMSI(2+LMSI) 1013. C**** 1014. C**** FLUXCB USIDT U compon of time integrated sea ice velocity (m) 1015. C**** VSIDT V compon of time integrated sea ice velocity (m) 1016. C**** 1017. C**** WORK01 FAW flux of surface water area (m^2) = USIDT*DYP 1018. C**** FASI flux of sea ice area (m^2) = USIDT*DYP*RSIedge 1018.5 C**** FMSI flux of sea ice mass (kg) or heat (J) 1019. C**** 1020. C**** North-South Advection of Sea Ice 1021. C**** 1022. SFASI = 0. 1023. DO 5 K=1,2+LMSI 1024. 5 SFMSI(K) = 0. 1027. DO 340 I=1,IM 1028. C**** 1029. C**** Calculate south-north sea ice fluxes at grid box edges 1030. C**** 1031. DO 120 J=2,JM-2 1032. IF(VSIDT(I,J).eq.0.) GO TO 120 1033. FAW(J) = VSIDT(I,J)*DXV(J) 1034. IF(VSIDT(I,J).gt.0.) GO TO 110 1035. C**** Sea ice velocity is southward at grid box edge 1036. FASI(J) = FAW(J)*(RSI(I,J+1)-(1.+FAW(J)*BYDXYP(J+1))*RSIY(I,J+1)) 1037. FXSI(J) = FAW(J)*RSIX(I,J+1) 1038. FYSI(J) = FAW(J)* 1039. * (FAW(J)*BYDXYP(J+1)*FAW(J)*RSIY(I,J+1) - 3.*FASI(J)) 1040. DO 105 K=1,2+LMSI 1041. 105 FMSI(K,J) = FASI(J)*MSI(I,J+1,K) 1043. AIJ(I,J,31) = AIJ(I,J,31) + (FMSI(1,J)+FMSI(2,J)) 1044. GO TO 120 1045. C**** Sea ice velocity is northward at grid box edge 1046. 110 FASI(J) = FAW(J)*(RSI(I,J)+(1.-FAW(J)*BYDXYP(J))*RSIY(I,J)) 1047. FXSI(J) = FAW(J)*RSIX(I,J) 1048. FYSI(J) = FAW(J)*(FAW(J)*BYDXYP(J)*FAW(J)*RSIY(I,J) - 3.*FASI(J)) 1049. DO 115 K=1,2+LMSI 1050. 115 FMSI(K,J) = FASI(J)*MSI(I,J,K) 1052. AIJ(I,J,31) = AIJ(I,J,31) + (FMSI(1,J)+FMSI(2,J)) 1053. 120 continue 1054. C**** 1055. C**** Calculate south-north sea ice fluxes near North Pole 1056. C**** 1057. IF(VSIDT(I,JM-1).eq.0.) GO TO 200 1058. FAW(JM-1) = VSIDT(I,JM-1)*DXV(JM-1) 1059. IF(VSIDT(I,JM-1).gt.0.) GO TO 130 1060. C**** Sea ice velocity is southward from North Pole box 1061. FASI(JM-1) = FAW(JM-1)*RSI(1,JM) 1062. FXSI(JM-1) = 0. 1063. FYSI(JM-1) = -FAW(JM-1)*FASI(JM-1) 1064. DO 125 K=1,2+LMSI 1065. 125 FMSI(K,JM-1) = FASI(JM-1)*MSI(1,JM,K) 1067. AIJ(I,JM-1,31) = AIJ(I,JM-1,31) + (FMSI(1,JM-1)+FMSI(2,JM-1)) 1068. GO TO 140 1069. C**** Sea ice velocity is northward into North Pole box 1070. 130 FASI(JM-1) = FAW(JM-1)* 1071. * (RSI(I,JM-1)+(1.-FAW(JM-1)*BYDXYP(JM-1))*RSIY(I,JM-1)) 1072. FXSI(JM-1) = FAW(JM-1)*RSIX(I,JM-1) 1073. FYSI(JM-1) = FAW(JM-1)* 1074. * (FAW(JM-1)*BYDXYP(JM-1)*FAW(JM-1)*RSIY(I,JM-1) - 3.*FASI(JM-1)) 1075. DO 135 K=1,2+LMSI 1076. 135 FMSI(K,JM-1) = FASI(JM-1)*MSI(I,JM-1,K) 1078. AIJ(I,JM-1,31) = AIJ(I,JM-1,31) + (FMSI(1,JM-1)+FMSI(2,JM-1)) 1079. C**** Accumulate sea ice leaving and entering North Pole box 1080. 140 SFASI = SFASI + FASI(JM-1) 1081. DO 145 K=1,2+LMSI 1082. 145 SFMSI(K) = SFMSI(K) + FMSI(K,JM-1) 1085. C**** 1086. C**** Update sea ice variables due to south-north fluxes 1087. C**** 1088. 200 DO 330 J=2,JM-1 1089. IF(VSIDT(I,J-1)) 240,210,280 1090. C**** VSIDT(J-1)=0. 1091. 210 IF(VSIDT(I,J)) 220,330,230 1092. C**** VSIDT(J-1)=0, VSIDT(J)<0. 1093. 220 ASI = RSI(I,J)*DXYP(J) - FASI(J) 1094. DO 225 K=1,2+LMSI 1095. 225 AMSI(K) = RSI(I,J)*DXYP(J)*MSI(I,J,K) - FMSI(K,J) 1098. IF(ASI.gt.DXYP(J)) GO TO 320 1099. YSI = (RSIY(I,J)*DXYP(J)*DXYP(J) - FYSI(J) 1100. * + 3.*(FAW(J)*ASI-DXYP(J)*FASI(J))) / (DXYP(J)-FAW(J)) 1101. RSI(I,J) = ASI*BYDXYP(J) 1102. RSIY(I,J) = YSI*BYDXYP(J) 1103. RSIX(I,J) = RSIX(I,J) - FXSI(J)*BYDXYP(J) 1104. DO 226 K=1,2+LMSI 1105. 226 MSI(I,J,K) = AMSI(K)/ASI 1108. GO TO 310 1109. C**** VSIDT(J-1)=0, VSIDT(J)>0. 1110. 230 RSI(I,J) = RSI(I,J) - FASI(J)*BYDXYP(J) 1111. RSIX(I,J) = RSIX(I,J)*(1.-FAW(J)*BYDXYP(J)) 1112. RSIY(I,J) = RSIY(I,J)*(1.-FAW(J)*BYDXYP(J))**2 1112.1 IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J) 1112.2 IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J) 1113. GO TO 330 1114. C**** VSIDT(J-1)<0. 1115. 240 IF(VSIDT(I,J)) 260,250,270 1116. C**** VSIDT(J-1)<0, VSIDT(J)=0. 1117. 250 RSI(I,J) = RSI(I,J) + FASI(J-1)*BYDXYP(J) 1118. RSIX(I,J) = RSIX(I,J)*(1.+FAW(J-1)*BYDXYP(J)) 1119. RSIY(I,J) = RSIY(I,J)*(1.+FAW(J-1)*BYDXYP(J))**2 1119.1 IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J) 1119.2 IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J) 1120. GO TO 330 1121. C**** VSIDT(J-1)<0, VSIDT(J)<0 or VSIDT(J-1)>0, VSIDT(J)Ý0. 1122. 260 ASI = RSI(I,J)*DXYP(J) + ( FASI(J-1)- FASI(J)) 1123. DO 265 K=1,2+LMSI 1124. 265 AMSI(K) = RSI(I,J)*DXYP(J)*MSI(I,J,K) + (FMSI(K,J-1)-FMSI(K,J)) 1127. IF(ASI.gt.DXYP(J)) GO TO 320 1128. YSI = (RSIY(I,J)*DXYP(J)*DXYP(J) + (FYSI(J-1)-FYSI(J)) 1129. * + 3.*((FAW(J-1)+FAW(J))*ASI-DXYP(J)*(FASI(J-1)+FASI(J)))) 1130. * / (DXYP(J) + (FAW(J-1)-FAW(J))) 1131. RSI(I,J) = ASI*BYDXYP(J) 1132. RSIY(I,J) = YSI*BYDXYP(J) 1133. RSIX(I,J) = RSIX(I,J) + (FXSI(J-1)-FXSI(J))*BYDXYP(J) 1134. DO 266 K=1,2+LMSI 1135. 266 MSI(I,J,K) = AMSI(K)/ASI 1138. GO TO 310 1139. C**** VSIDT(J-1)<0, VSIDT(J)>0. 1140. 270 RSI(I,J) = RSI(I,J) + (FASI(J-1)-FASI(J))*BYDXYP(J) 1141. RSIX(I,J) = RSIX(I,J)*(1.+(FAW(J-1)-FAW(J))*BYDXYP(J)) 1142. RSIY(I,J) = RSIY(I,J)*(1.+(FAW(J-1)-FAW(J))*BYDXYP(J))**2 1142.1 IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J) 1142.2 IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J) 1143. GO TO 330 1144. C**** VSIDT(J-1)>0. 1145. 280 IF(VSIDT(I,J).ne.0.) GO TO 260 1146. C**** VSIDT(J-1)>0, VSIDT(J)=0. 1147. ASI = RSI(I,J)*DXYP(J) + FASI(J-1) 1148. DO 285 K=1,2+LMSI 1149. 285 AMSI(K) = RSI(I,J)*DXYP(J)*MSI(I,J,K) + FMSI(K,J-1) 1152. IF(ASI.gt.DXYP(J)) GO TO 320 1153. YSI = (RSIY(I,J)*DXYP(J)*DXYP(J) + FYSI(J-1) 1154. * + 3.*(FAW(J-1)*ASI-DXYP(J)*FASI(J-1))) / (DXYP(J)+FAW(J-1)) 1155. RSI(I,J) = ASI*BYDXYP(J) 1156. RSIY(I,J) = YSI*BYDXYP(J) 1157. RSIX(I,J) = RSIX(I,J) + FXSI(J-1)*BYDXYP(J) 1158. DO 286 K=1,2+LMSI 1159. 286 MSI(I,J,K) = AMSI(K)/ASI 1162. C**** Limit RSIX and RSIY so that sea ice is positive at the edges 1163. 310 IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J) 1164. IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J) 1165. IF(RSI(I,J)-RSIX(I,J).gt.1.) RSIX(I,J) = RSI(I,J)-1. 1166. IF(RSI(I,J)+RSIX(I,J).gt.1.) RSIX(I,J) = 1.-RSI(I,J) 1167. IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J) 1168. IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J) 1169. IF(RSI(I,J)-RSIY(I,J).gt.1.) RSIY(I,J) = RSI(I,J)-1. 1170. IF(RSI(I,J)+RSIY(I,J).gt.1.) RSIY(I,J) = 1.-RSI(I,J) 1171. GO TO 330 1172. C**** Sea ice crunches into itself and completely covers grid box 1173. 320 RSI(I,J) = 1. 1174. RSIX(I,J) = 0. 1175. RSIY(I,J) = 0. 1176. MSI(I,J,1) = AMSI(1)/ASI 1176.5 MSI(I,J,2) =(AMSI(1)+AMSI(2))*BYDXYP(J) - MSI(I,J,1) 1177. HSI(I,J,1) = AMSI(3)/ASI 1177.5 HSI(I,J,2) = AMSI(4)/ASI 1178. DHSI = (AMSI(3)+AMSI(4)+AMSI(5)+AMSI(6))*(BYDXYP(J) - 1./ASI) 1179. HSI(I,J,3) = AMSI(5)/ASI + XSI3*DHSI 1179.5 HSI(I,J,4)= AMSI(6)/ASI + XSI4*DHSI 1180. C**** End of loop over J 1181. 330 continue 1182. C**** End of loop over I 1183. 340 continue 1184. C**** 1185. C**** Advection of Sea Ice leaving and entering North Pole box 1186. C**** 1187. ASI = RSI(1,JM)*DXYP(JM) + SFASI/IM 1188. DO 345 K=1,2+LMSI 1189. 345 AMSI(K) = RSI(1,JM)*DXYP(JM)*MSI(1,JM,K) + SFMSI(K)/IM 1192. IF(ASI.gt.DXYP(JM)) GO TO 350 1193. RSI(1,JM) = ASI*BYDXYP(JM) 1194. DO 346 K=1,2+LMSI 1195. 346 MSI(1,JM,K) = AMSI(K)/ASI 1198. GO TO 400 1199. C**** Sea ice crunches into itself at North Pole box 1200. 350 RSI(1,JM) = 1. 1201. MSI(1,JM,1) = AMSI(1)/ASI 1201.5 MSI(1,JM,2) =(AMSI(1)+AMSI(2))*BYDXYP(JM) - MSI(1,JM,1) 1202. HSI(1,JM,1) = AMSI(3)/ASI 1202.5 HSI(1,JM,2) = AMSI(4)/ASI 1203. DHSI = (AMSI(3)+AMSI(4)+AMSI(5)+AMSI(6))*(BYDXYP(JM) - 1./ASI) 1204. HSI(1,JM,3) = AMSI(5)/ASI + XSI3*DHSI 1204.5 HSI(1,JM,4)= AMSI(6)/ASI + XSI4*DHSI 1205. C**** 1206. C**** East-West Advection of Sea Ice 1207. C**** 1208. 400 DO 640 J=2,JM-1 1209. C**** 1210. C**** Calculate west-east sea ice fluxes at grid box edges 1211. C**** 1212. I=IM 1213. DO 420 IP1=1,IM 1214. IF(USIDT(I,J).eq.0.) GO TO 420 1215. FAW(I) = USIDT(I,J)*DYP(J) 1216. IF(USIDT(I,J).gt.0.) GO TO 410 1217. C**** Sea ice velocity is westward at grid box edge 1218. FASI(I) = FAW(I)*(RSI(IP1,J)-(1.+FAW(I)*BYDXYP(J))*RSIX(IP1,J)) 1219. FXSI(I) = FAW(I)*(FAW(I)*BYDXYP(J)*FAW(I)*RSIX(IP1,J)-3.*FASI(I)) 1220. FYSI(I) = FAW(I)*RSIY(IP1,J) 1221. DO 405 K=1,2+LMSI 1222. 405 FMSI(K,I) = FASI(I)*MSI(IP1,J,K) 1224. AIJ(I,J,30) = AIJ(I,J,30) + (FMSI(1,I)+FMSI(2,I)) 1225. GO TO 420 1226. C**** Sea ice velocity is eastward at grid box edge 1227. 410 FASI(I) = FAW(I)*(RSI(I,J)+(1.-FAW(I)*BYDXYP(J))*RSIX(I,J)) 1228. FXSI(I) = FAW(I)*(FAW(I)*BYDXYP(J)*FAW(I)*RSIX(I,J) - 3.*FASI(I)) 1229. FYSI(I) = FAW(I)*RSIY(I,J) 1230. DO 415 K=1,2+LMSI 1231. 415 FMSI(K,I) = FASI(I)*MSI(I,J,K) 1233. AIJ(I,J,30) = AIJ(I,J,30) + (FMSI(1,I)+FMSI(2,I)) 1234. 420 I=IP1 1235. C**** 1236. C**** Update sea ice variables due to west-east fluxes 1237. C**** 1238. IM1=IM 1239. DO 630 I=1,IM 1240. IF(USIDT(IM1,J)) 540,510,580 1241. C**** USIDT(IM1)=0. 1242. 510 IF(USIDT(I,J)) 520,630,530 1243. C**** USIDT(IM1)=0, USIDT(I)<0. 1244. 520 ASI = RSI(I,J)*DXYP(J) - FASI(I) 1245. DO 525 K=1,2+LMSI 1246. 525 AMSI(K) = RSI(I,J)*DXYP(J)*MSI(I,J,K) - FMSI(K,I) 1249. IF(ASI.gt.DXYP(J)) GO TO 620 1250. XSI = (RSIX(I,J)*DXYP(J)*DXYP(J) - FXSI(I) 1251. * + 3.*(FAW(I)*ASI-DXYP(J)*FASI(I))) / (DXYP(J)-FAW(I)) 1252. RSI(I,J) = ASI*BYDXYP(J) 1253. RSIX(I,J) = XSI*BYDXYP(J) 1254. RSIY(I,J) = RSIY(I,J) - FYSI(I)*BYDXYP(J) 1255. DO 526 K=1,2+LMSI 1256. 526 MSI(I,J,K) = AMSI(K)/ASI 1259. GO TO 610 1260. C**** USIDT(IM1)=0, USIDT(I)>0. 1261. 530 RSI(I,J) = RSI(I,J) - FASI(I)*BYDXYP(J) 1262. RSIX(I,J) = RSIX(I,J)*(1.-FAW(I)*BYDXYP(J))**2 1263. RSIY(I,J) = RSIY(I,J)*(1.-FAW(I)*BYDXYP(J)) 1263.1 IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J) 1263.2 IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J) 1264. GO TO 630 1265. C**** USIDT(IM1)<0. 1266. 540 IF(USIDT(I,J)) 560,550,570 1267. C**** USIDT(IM1)<0, USIDT(I)=0. 1268. 550 RSI(I,J) = RSI(I,J) + FASI(IM1)*BYDXYP(J) 1269. RSIX(I,J) = RSIX(I,J)*(1.+FAW(IM1)*BYDXYP(J))**2 1270. RSIY(I,J) = RSIY(I,J)*(1.+FAW(IM1)*BYDXYP(J)) 1270.1 IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J) 1270.2 IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J) 1271. GO TO 630 1272. C**** USIDT(IM1)<0, USIDT(I)<0 or USIDT(IM1)>0, USIDT(I)Ý0. 1273. 560 ASI = RSI(I,J)*DXYP(J) + (FASI(IM1)- FASI(I)) 1274. DO 565 K=1,2+LMSI 1275. 565 AMSI(K) = RSI(I,J)*DXYP(J)*MSI(I,J,K) + (FMSI(K,IM1)-FMSI(K,I)) 1278. IF(ASI.gt.DXYP(J)) GO TO 620 1279. XSI = (RSIX(I,J)*DXYP(J)*DXYP(J) + (FXSI(IM1)-FXSI(I)) 1280. * + 3.*((FAW(IM1)+FAW(I))*ASI-DXYP(J)*(FASI(IM1)+FASI(I)))) 1281. * / (DXYP(J) + (FAW(IM1)-FAW(I))) 1282. RSI(I,J) = ASI*BYDXYP(J) 1283. RSIX(I,J) = XSI*BYDXYP(J) 1284. RSIY(I,J) = RSIY(I,J) + (FYSI(IM1)-FYSI(I))*BYDXYP(J) 1285. DO 566 K=1,2+LMSI 1286. 566 MSI(I,J,K) = AMSI(K)/ASI 1289. GO TO 610 1290. C**** USIDT(IM1)<0, USIDT(I)>0. 1291. 570 RSI(I,J) = RSI(I,J) + (FASI(IM1)-FASI(I))*BYDXYP(J) 1292. RSIX(I,J) = RSIX(I,J)*(1.+(FAW(IM1)-FAW(I))*BYDXYP(J))**2 1293. RSIY(I,J) = RSIY(I,J)*(1.+(FAW(IM1)-FAW(I))*BYDXYP(J)) 1293.1 IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J) 1293.2 IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J) 1294. GO TO 630 1295. C**** USIDT(IM1)>0. 1296. 580 IF(USIDT(I,J).ne.0.) GO TO 560 1297. C**** USIDT(IM1)>0, USIDT(I)=0. 1298. ASI = RSI(I,J)*DXYP(J) + FASI(IM1) 1299. DO 585 K=1,2+LMSI 1300. 585 AMSI(K) = RSI(I,J)*DXYP(J)*MSI(I,J,K) + FMSI(K,IM1) 1303. IF(ASI.gt.DXYP(J)) GO TO 620 1304. XSI = (RSIX(I,J)*DXYP(J)*DXYP(J) + FXSI(IM1) 1305. * + 3.*(FAW(IM1)*ASI-DXYP(J)*FASI(IM1))) / (DXYP(J)+FAW(IM1)) 1306. RSI(I,J) = ASI*BYDXYP(J) 1307. RSIX(I,J) = XSI*BYDXYP(J) 1308. RSIY(I,J) = RSIY(I,J) + FYSI(IM1)*BYDXYP(J) 1309. DO 586 K=1,2+LMSI 1310. 586 MSI(I,J,K) = AMSI(K)/ASI 1313. C**** Limit RSIX and RSIY so that sea ice is positive at the edges 1314. 610 IF(RSI(I,J)-RSIX(I,J).lt.0.) RSIX(I,J) = RSI(I,J) 1315. IF(RSI(I,J)+RSIX(I,J).lt.0.) RSIX(I,J) = -RSI(I,J) 1316. IF(RSI(I,J)-RSIX(I,J).gt.1.) RSIX(I,J) = RSI(I,J)-1. 1317. IF(RSI(I,J)+RSIX(I,J).gt.1.) RSIX(I,J) = 1.-RSI(I,J) 1318. IF(RSI(I,J)-RSIY(I,J).lt.0.) RSIY(I,J) = RSI(I,J) 1319. IF(RSI(I,J)+RSIY(I,J).lt.0.) RSIY(I,J) = -RSI(I,J) 1320. IF(RSI(I,J)-RSIY(I,J).gt.1.) RSIY(I,J) = RSI(I,J)-1. 1321. IF(RSI(I,J)+RSIY(I,J).gt.1.) RSIY(I,J) = 1.-RSI(I,J) 1322. GO TO 630 1323. C**** Sea ice crunches into itself and completely covers grid box 1324. 620 RSI(I,J) = 1. 1325. RSIX(I,J) = 0. 1326. RSIY(I,J) = 0. 1327. MSI(I,J,1) = AMSI(1)/ASI 1327.5 MSI(I,J,2) =(AMSI(1)+AMSI(2))*BYDXYP(J) - MSI(I,J,1) 1328. HSI(I,J,1) = AMSI(3)/ASI 1328.5 HSI(I,J,2) = AMSI(4)/ASI 1329. DHSI = (AMSI(3)+AMSI(4)+AMSI(5)+AMSI(6))*(BYDXYP(J) - 1./ASI) 1330. HSI(I,J,3) = AMSI(5)/ASI + XSI3*DHSI 1330.5 HSI(I,J,4) = AMSI(6)/ASI + XSI4*DHSI 1331. C**** End of loop over I 1332. 630 IM1=I 1333. C**** End of loop over J 1334. 640 continue 1335. RETURN 1336. END 2000. 2001. SUBROUTINE PRECSI 2002. C**** 2003. C**** PRECSI adds atmospheric precipitation to sea 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 ocean mass (kg^2) 2010. C**** 2011. C**** SEAICE RSI Horizontal ratio of sea ice to ocean (1) 2012. C**** MSI Mass of sea ice (kg/m^2) 2013. C**** HSI Enthalpy minus latent heat of sea ice (J/m^2) 2014. C**** 2015. C**** FIXDCB FOCEAN Ocean fraction (0 or 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. IF(PREC(I,J,1).le.0. .or. FOCEAN(I,J)*RSI(I,J).le.0.) GO TO 500 2027. FICE = FOCEAN(I,J)*RSI(I,J) 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) 2044. EJ(J,21) = EJ(J,21) + FICE*(RAIN+MELT1)*ZATMO(I,J) 2045. FJ(J,21) = FJ(J,21) - FICE*(RAIN+MELT1)*ZATMO(I,J) 2046. C EJ(J,33) = EJ(J,33) + FICE*(RAIN+MELT1) ! in DIAGJ 2047. FJ(J,33) = FJ(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) 2078. EJ(J,21) = EJ(J,21) + FICE*(RAIN-FREZ1)*ZATMO(I,J) 2079. FJ(J,21) = FJ(J,21) - FICE*(RAIN-FREZ1)*ZATMO(I,J) 2080. C EJ(J,33) = EJ(J,33) + FICE*(RAIN-FREZ1) ! in DIAGJ 2081. FJ(J,33) = FJ(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) 2108. EJ(J,21) = EJ(J,21) + FICE*(RAIN-FREZ1)*ZATMO(I,J) 2109. FJ(J,21) = FJ(J,21) - FICE*(RAIN-FREZ1)*ZATMO(I,J) 2110. C EJ(J,33) = EJ(J,33) + FICE*(RAIN-FREZ1) ! in DIAGJ 2111. FJ(J,33) = FJ(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 SEAICE 3002. C**** 3003. C**** SEAICE receives the atmospheric fluxes of heat and water and 3004. C**** applies them to sea ice, updating the heat, mass and snow 3005. C**** 3006. INCLUDE 'C070.COM' 3007. PARAMETER (RHOI=910., ALAMI=2.1762d0, ALPHA=1., 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 ocean mass (kg/m^2) 3013. C**** 3014. C**** SEAICB RSI Horizontal ratio of sea ice to ocean (1) 3015. C**** MSI Mass of sea ice (kg/m^2) 3016. C**** HSI Enthalpy minus latent heat of sea ice (J/m^2) 3017. C**** 3018. C**** FIXDCB FOCEAN Ocean fraction (0 or 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) 3022.5 C**** E1(1) Energy into ocean from bottom of ice (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. IF(FOCEAN(I,J)*RSI(I,J) .le. 0.) GO TO 400 3031. FICE = FOCEAN(I,J)*RSI(I,J) 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 3035.1 GO = G0M(I,J,1)/(MO(I,J,1)*DXYP(J)) 3035.2 SO = S0M(I,J,1)/(MO(I,J,1)*DXYP(J)) 3035.3 TCO = TEMGS(GO,SO) 3036. C**** Accumulate diagnostics 3037. FJ(J,30) = FJ(J,30) + FICE 3038. FJ(J,51) = FJ(J,51) + FICE*MSI(I,J,2) 3039. FJ(J,18) = FJ(J,18) + FICE*TSI1 3040. FJ(J,17) = FJ(J,17) + FICE*TSI2 3041. FJ(J,16) = FJ(J,16) + FICE*TSI3 3042. FJ(J,52) = FJ(J,52) + FICE*TSI4 3043. FJ(J,19) = FJ(J,19) + FICE*W0(I,J,2) 3044. C FJ(J,22) = FJ(J,22) + FICE*W0(I,J,2)*ZATMO(I,J) ! for Caspian 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. FJ(J,31) = FJ(J,31) + FICE 3050. FJ(J,53) = FJ(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. E1(I,J,1) = dF4dTI*(HCI4*(TSI4-TCO ) + ALPHA*F3) / 3068.1 / (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 - E1(I,J,1)) 3073. C G0M(I,J,1) = G0M(I,J,1) + FICE*E1(I,J,1)*DXYP(J) ! in OSOURC 3074. FJ(J,15) = FJ(J,15) - FICE*E1(I,J,1) 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 3081. C EJ(J,22) = EJ(J,22) + FICE*MELT1*ZATMO(I,J) ! for Caspian 3082. C FJ(J,22) = FJ(J,22) - FICE*MELT1*ZATMO(I,J) ! for Caspian 3083. C EJ(J,34) = EJ(J,34) + FICE*MELT1 ! in DIAGJ 3084. FJ(J,34) = FJ(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 3142. C EJ(J,22) = EJ(J,22) + FICE*MELT4*ZATMO(I,J) ! for Caspian 3143. C FJ(J,22) = FJ(J,22) - FICE*MELT4*ZATMO(I,J) ! for Caspian 3144. C EJ(J,38) = EJ(J,38) + FICE*MELT4 ! in DIAGJ 3145. FJ(J,38) = FJ(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 3170. C EJ(J,22) = EJ(J,22) + FICE*MELT4*ZATMO(I,J) ! for Caspian 3171. C FJ(J,22) = FJ(J,22) - FICE*MELT4*ZATMO(I,J) ! for Caspian 3172. C EJ(J,38) = EJ(J,38) + FICE*MELT4 ! in DIAGJ 3173. FJ(J,38) = FJ(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 3192. C EJ(J,22) = EJ(J,22) + FICE*MELT4*ZATMO(I,J) ! for Caspian 3193. C FJ(J,22) = FJ(J,22) - FICE*MELT4*ZATMO(I,J) ! for Caspian 3194. C EJ(J,38) = EJ(J,38) + FICE*MELT4 ! in DIAGJ 3195. FJ(J,38) = FJ(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 MOVEIB 4002. C**** 4003. C**** MOVEIB moves icebergs depending on ocean currents and presence 4004. C**** of sea ice 4005. C**** 4006. INCLUDE 'C070.COM' 4007. PARAMETER (LMIB=5, RHOA=1.27556d0, RHOI=910., MSIX=0.) 4008. REAL*8 OMF(LMIB) 4009. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 4010. * W0(IM,JM,4), E0(IM,JM,4),E1(IM,JM,4), UAS(IM,JM),VAS(IM,JM) 4011. DATA IFIRST /1/ 4012. C**** 4013. C**** SICECB BERGM mass of each ice berg (kg) 4014. C**** BERGH enthalpy minus latent energy of each ice berg (J) 4015. C**** BERGI longitude of each ice berg (1. = 175W) 4016. C**** BERGJ latitude of each ice berg (1. = 88S) 4017. C**** 4018. C**** FLUXCB UAS U atmospheric surface wind (m/s) on A grid 4019. C**** VAS V atmospheric surface wind (m/s) on A grid 4020. C**** 4021. PARAMETER (ACDxRHOAzMIB=.0000005d0*RHOA, 4022. * OCDxRHOOzMIB=.001667d0, DW=.25, OD=DW*OCDxRHOOzMIB) 4023. C**** 4024. C**** MIB (kg/m^2) = ice berg mass per unit area = RHOI*ZE(LMIB) 4025. C**** ACD = air - ice berg drag coefficient = .0000005*MIB 4026. C**** OCD = ocean - ice berg drag coefficient = .001667*MIB/RHOO 4027. C**** 4028. C**** AD (1/s) = ACD*RHOA*abs(WAS-WIB)/MIB 4029. C**** OD (1/s) = OCD*RHOO*abs(WOL-WIB)/MIB = OCD*RHOO*.25/MIB 4030. C**** SD (1/s) = OD*MSIX*MSI/MIB 4031. C**** 4032. C**** UIB = (AD*UAS + OD*UOB + SD*USI) / (AD + OD + SD) 4033. C**** VIB = (AD*VAS + OD*VOB + SD*VSI) / (AD + OD + SD) 4034. C**** 4035. IF(IFIRST.eq.0) GO TO 100 4036. IFIRST = 0 4037. DLON = TWOPI/IM 4038. DLAT = TWOPI*NINT(360./(JM-1))/720. 4039. zMIB = 1./(RHOI*ZE(LMIB)) 4040. C**** Calculate ocean mass fraction OMF of each of first LMIB layers 4041. DO 10 L=1,LMIB 4042. 10 OMF(L) = DSIGO(L) / (SIGEO(LMIB)-SIGEO(0)) 4043. C**** 4044. C**** Loop over Ice Bergs 4045. C**** 4046. 100 DO 410 N=1,NBERG 4047. I = BERGI(N) + 1. 4048. J = BERGJ(N) + 1. 4049. IF(LMM(I,J).lt.LMIB) GO TO 400 ! ice berg is grounded 4050. IF(J.eq.JM) GO TO 300 4051. C**** 4052. C**** Calculate ice berg movement away from the North Pole 4053. C**** 4054. AD = ACDxRHOAzMIB * SQRT(UAS(I,J)*UAS(I,J)+VAS(I,J)*VAS(I,J)) 4055. C**** Calculate vertically integrated ocean velocity 4056. Im1 = I-1 4057. IF(Im1.lt.1) Im1 = IM 4058. UOL = 0. 4059. VOL = 0. 4060. DO 110 L=1,LMIB 4061. UOL = UOL + OMF(L)*(UO(Im1,J,L)*(I-BERGI(N)) + 4062. + UO(I ,J,L)*(BERGI(N)+1.-I)) 4063. 110 VOL = VOL + OMF(L)*(VO(I,J-1,L)*(J-BERGJ(N)) + 4064. + VO(I,J ,L)*(BERGJ(N)+1.-J)) 4065. C**** Calculate sea ice velocity 4066. SD = 0. 4067. RSIB = RSI(I,J) + (BERGI(N)+.5-I)*RSIX(I,J) + 4068. + (BERGJ(N)+.5-J)*RSIY(I,J) 4069. IF(RSIB.le.0.) GO TO 120 4070. SD = OD*MSIX*RSIB*(ACE1I+MSI(I,J,2))*zMIB 4071. USIN = USI(Im1,J)*(I-BERGI(N)) + USI(I,J)*(BERGI(N)+1.-I) 4072. VSIN = VSI(I,J-1)*(J-BERGJ(N)) + VSI(I,J)*(BERGJ(N)+1.-J) 4073. C**** Calculate ice berg velocity 4074. 120 UIB = (UAS(I,J)*AD + UOL*OD + USIN*SD) / (AD + OD + SD) 4075. VIB = (VAS(I,J)*AD + VOL*OD + VSIN*SD) / (AD + OD + SD) 4076. C**** Move Ice Berg 4077. BERGJX = BERGJ(N) +.5*DTS*VIB/(DLAT*RADIUS) 4078. BERGJ(N) = BERGJ(N) + DTS*VIB/(DLAT*RADIUS) 4079. BERGI(N) = BERGI(N) + DTS*UIB/(DLON*RADIUS* 4080. * (COSV(J-1)*(J-BERGJX) + COSV(J)*(BERGJX+1.-J))) 4081. IF(BERGI(N).lt.0.) BERGI(N) = BERGI(N) + IM 4082. IF(BERGI(N).ge.IM) BERGI(N) = BERGI(N) - IM 4083. C**** 4084. C**** Check that ice berg did not leave the ocean 4085. C**** 4086. Inew = BERGI(N) + 1. 4087. Jnew = BERGJ(N) + 1. 4088. IF(LMM(Inew,Jnew).gt.0) GO TO 410 4089. C**** Ice berg left the ocean, bring it back 4090. IF(LMM(Inew,J).gt.0) GO TO 210 4091. IF(LMM(I,Jnew).gt.0) GO TO 220 4092. C**** Ice berg left the ocean at a corner, bring it back 4093. IF(UIB.lt.0.) then 4094. BERGI(N) = I - .999d0 4095. else 4096. BERGI(N) = I - .001d0 4097. endif 4098. C**** Bring ice berg back to ocean either north or south of current box 4099. 210 IF(VIB.lt.0.) then 4100. BERGJ(N) = J - .999d0 4101. else 4102. BERGJ(N) = J - .001d0 4103. endif 4104. GO TO 410 4105. C**** Bring ice berg back to ocean either east or west of current box 4106. 220 IF(UIB.lt.0.) then 4107. BERGI(N) = I - .999d0 4108. else 4109. BERGI(N) = I - .001d0 4110. endif 4111. GO TO 410 4112. C**** 4113. C**** Calculate ice berg movement near the North Pole (J=JM) 4114. C**** 4115. C**** Locate ice berg in unit polar circle with rectangular coordinates 4116. 300 X = (JM-BERGJ(N))*COS(TWOPI*(BERGI(N)/IM-.25)) 4117. Y = (JM-BERGJ(N))*SIN(TWOPI*(BERGI(N)/IM-.25)) 4118. AD = ACDxRHOAzMIB * SQRT(UAS(1,JM)*UAS(1,JM)+VAS(1,JM)*VAS(1,JM)) 4119. C**** Calculate vertically integrated ocean velocity 4120. I1 = 1 4121. I2 = IM/4 4122. I3 = IM/4 + 1 4123. I4 = IM/2 4124. I5 = IM/2 + 1 4125. I6 = 3*IM/4 4126. I7 = 3*IM/4 + 1 4127. I8 = IM 4128. UOL = 0. 4129. VOL = 0. 4130. DO 310 L=1,LMIB 4131. UOL = UOL + .25*OMF(L)*((VO(I6,JM-1,L)+VO(I7,JM-1,L))*(1.-X) - 4132. - (VO(I2,JM-1,L)+VO(I3,JM-1,L))*(1.+X)) 4133. 310 VOL = VOL + .25*OMF(L)*((VO(I8,JM-1,L)+VO(I1,JM-1,L))*(1.-Y) - 4134. - (VO(I4,JM-1,L)+VO(I5,JM-1,L))*(1.+Y)) 4135. C**** Calculate sea ice velocity 4136. SD = 0. 4137. IF(RSI(1,JM).le.0.) GO TO 320 4138. SD = OD*MSIX*RSI(1,JM)*(ACE1I+MSI(1,JM,2))*zMIB 4139. USIN = .25*((VSI(I6,JM-1)+VSI(I7,JM-1))*(1.-X) - 4140. - (VSI(I2,JM-1)+VSI(I3,JM-1))*(1.+X)) 4141. VSIN = .25*((VSI(I8,JM-1)+VSI(I1,JM-1))*(1.-Y) - 4142. - (VSI(I4,JM-1)+VSI(I5,JM-1))*(1.+Y)) 4143. C**** Calculate ice berg velocity 4144. 320 UIB = (UAS(1,JM)*AD + UOL*OD + USIN*SD) / (AD + OD + SD) 4145. VIB = (VAS(1,JM)*AD + VOL*OD + VSIN*SD) / (AD + OD + SD) 4146. C**** Move ice berg (change unit circle to circle of radius DYP(JM)) 4147. X = X*DYP(JM) + DTS*UIB 4148. Y = Y*DYP(JM) + DTS*VIB 4149. C**** Convert polar rectangular coordinates back to I and J 4150. BERGI(N) = IM*ATAN2(Y,X)/TWOPI + IM/4 4151. IF(BERGI(N).lt.0.) BERGI(N) = BERGI(N) + IM 4152. Z = SQRT(X*X + Y*Y) 4153. IF(Z.le.DYP(JM)) BERGJ(N) = JM - Z/DYP(JM) 4154. IF(Z.gt.DYP(JM)) BERGJ(N) = JM-1 - (Z-DYP(JM))/DYP(JM-1) 4155. GO TO 410 4156. C**** 4157. C**** Ice berg is grounded, check that it is still in the ocean 4158. C**** 4159. 400 IF(LMM(I,J).gt.0) GO TO 410 4160. WRITE (6,*) 'Ice berg ran into continent. N,I,J,BERG* =', 4161. * N,I,J,BERGM(N),BERGH(N),BERGI(N),BERGJ(N) 4162. STOP 'MOVEIB 400:' 4163. 410 continue 4164. END 5000. 5001. SUBROUTINE MELTSI 5002. C**** 5003. C**** MELTSI melts sea ice horizontally and vertically, each using half 5004. C**** the energy, if the ocean temperature is above 0øC 5005. C**** 5006. INCLUDE 'C070.COM' 5007. PARAMETER (ACE2SI=182., SHCO=4000., FDAILY=.3333333d0, R4=1/XSI4) 5008. C**** 5009. C**** OCENCB MO Mass of ocean (kg/m^2) 5010. C**** G0M Mean potential enthalpy of ocean (J) 5011. C**** S0M Mean salt of ocean (kg) 5012. C**** 5013. C**** SICECB RSI Horizontal ratio of sea ice to ocean (1) 5014. C**** MSI Mass of sea ice (kg/m^2) 5015. C**** HSI Enthalpy minus latent heat of sea ice (J/m^2) 5016. C**** 5017. C**** ENRG (J/m^2) = energy moved from sea ice to ocean < 0 5018. C**** DMO (kg/m^2) = mass moved from sea ice to ocean > 0 5019. C**** MELT(kg/m^2) = mass of sea ice melted vertically > 0 5020. C**** DRSI = loss of horizontal sea ice cover > 0 5021. C**** 5022. DO 20 J=2,JM 5023. IMAX=IM 5024. IF(J.eq.JM) IMAX=1 5025. DO 20 I=1,IMAX 5026. IF(RSI(I,J).le.0. .or. FOCEAN(I,J).le.0.) GO TO 20 5027. IF(RSI(I,J)*(MSI(I,J,1)+MSI(I,J,2)) .lt. 3600./ELHM) GO TO 10 5028. GO = G0M(I,J,1)/(MO(I,J,1)*DXYP(J)) 5029. SO = S0M(I,J,1)/(MO(I,J,1)*DXYP(J)) 5030. TO = TEMGS(GO,SO) 5031. IF(TO.le.0.) GO TO 20 5032. C**** Ocean is above 0, calculate minus energy available for ice melting 5033. ENRG = -(SHCO*FDAILY*DTS/SDAY)*TO*MO(I,J,1) 5034. IF(ENRG-RSI(I,J)*(HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4)) 5035. * .le. 0.) GO TO 10 5036. C**** Melt vertically when RSI = 1, melt horizontally when RSI is near 0 5037. MELT = XSI4*MSI(I,J,2)*ENRG/HSI(I,J,4) 5038. IF(MSI(I,J,2)-MELT .lt. ACE2SI) MELT = MSI(I,J,2) - ACE2SI 5039. C FMSI3 = XSI3*MELT ! > 0 5040. FHSI3 = HSI(I,J,3)*MELT / MSI(I,J,2) 5041. FHSI4 = HSI(I,J,4)*MELT*R4 / MSI(I,J,2) 5042. C**** Melt some sea ice horizontally 5043. DRSI = (ENRG - RSI(I,J)*FHSI4) / 5044. / (HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4)-FHSI4) 5045. DMO = RSI(I,J)*MELT + DRSI*(MSI(I,J,1)+MSI(I,J,2)-MELT) 5046. C GXM(I,J,1) = 5047. C GYM(I,J,1) = 5048. C GZM(I,J,1) = GZM(I,J,1) - 3.*ENRG*DXYP(J) + 5049. C + (3.*G0M(I,J,1)-GZM(I,J,1))*DMO/MO(I,J,1) 5050. C SXM(I,J,1) = 5051. C SYM(I,J,1) = 5052. C SZM(I,J,1) = SZM(I,J,1) + (3.*S0M(I,J,1)-SZM(I,J,1))*DMO/MO(I,J,1) 5053. G0M(I,J,1) = G0M(I,J,1) + ENRG*DXYP(J) 5054. MO(I,J,1) = MO(I,J,1) + DMO 5055. C UO(I,J,1) = 5056. C VO(I,J,1) = 5057. FJ(J,36) = FJ(J,36) - DMO 5058. FJ(J,45) = FJ(J,45) - ENRG 5059. AIJ(I,J,63) = AIJ(I,J,63) + DMO 5060. AIJ(I,J,64) = AIJ(I,J,64) + ENRG 5061. RSIX(I,J) = RSIX(I,J)*(1.-DRSI/RSI(I,J)) 5062. RSIY(I,J) = RSIY(I,J)*(1.-DRSI/RSI(I,J)) 5063. RSI(I,J) = RSI(I,J) - DRSI 5064. MSI(I,J,2) = MSI(I,J,2) - MELT 5065. HSI(I,J,3) = HSI(I,J,3) - FHSI3 5066. HSI(I,J,4) = HSI(I,J,4) + FHSI3 - FHSI4 5067. GO TO 20 5068. C**** All sea ice melts 5069. 10 DMO = RSI(I,J)*(MSI(I,J,1)+MSI(I,J,2)) 5070. ENRG = RSI(I,J)*(HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4)) 5071. C GXM(I,J,1) = 5072. C GYM(I,J,1) = 5073. C GZM(I,J,1) = GZM(I,J,1) - 3.*ENRG*DXYP(J) + 5074. C + (3.*G0M(I,J,1)-GZM(I,J,1))*DMO/MO(I,J,1) 5075. C SXM(I,J,1) = 5076. C SYM(I,J,1) = 5077. C SZM(I,J,1) = SZM(I,J,1) + (3.*S0M(I,J,1)-SZM(I,J,1))*DMO/MO(I,J,1) 5078. G0M(I,J,1) = G0M(I,J,1) + ENRG*DXYP(J) 5079. MO(I,J,1) = MO(I,J,1) + DMO 5080. C UO(I,J,1) = 5081. C VO(I,J,1) = 5082. FJ(J,36) = FJ(J,36) - DMO 5083. FJ(J,45) = FJ(J,45) - ENRG 5084. AIJ(I,J,63) = AIJ(I,J,63) + DMO 5085. AIJ(I,J,64) = AIJ(I,J,64) + ENRG 5086. RSI(I,J) = 0. 5087. RSIX(I,J) = 0. 5088. RSIY(I,J) = 0. 5089. PSI(I,J) = 0. 5090. 20 continue 5091. RETURN 5092. END 6000. 6001. SUBROUTINE MELTIB 6002. C**** 6003. C**** MELTIB melts icebergs depending on the ocean temperature 6004. C**** 6005. INCLUDE 'C070.COM' 6006. PARAMETER (LMIB=5, XMELT=2.) 6007. REAL*8 OMF(LMIB), BYMM(LMIB), EFLUX(LMIB), MMELT(LMIB) 6008. DATA IFIRST /1/ 6009. C**** 6010. C**** SICECB BERGM mass of each ice berg (kg) 6011. C**** BERGH enthalpy minus latent energy of each ice berg (J) 6012. C**** BERGI longitude of each ice berg on U-grid (1. = 175W) 6013. C**** BERGJ latitude of each ice berg on V-grid (1. = 88S) 6014. C**** 6015. IF(IFIRST.eq.0) GO TO 100 6016. IFIRST = 0 6017. C**** Calculate ocean mass fraction of each of first LMIB layers 6018. DO 10 L=1,LMIB 6019. 10 OMF(L) = DSIGO(L) / (SIGEO(LMIB)-SIGEO(0)) 6020. C**** 6021. C**** Melt all or part of ice berg 6022. C**** 6023. 100 N = 0 6024. 110 N = N+1 6025. 120 IF(N.gt.NBERG) GO TO 500 6026. SRM = SQRT(BERGM(N)) 6027. TCI = (BERGH(N)/BERGM(N)+ELHM)/SHCI 6028. I = BERGI(N) + 1. 6029. J = BERGJ(N) + 1. 6030. X = 2.*BERGI(N) - 2*I + 1. 6031. Y = 2.*BERGJ(N) - 2*J + 1. 6032. IF(J.eq.JM) I = 1 6033. LMAX = MIN (LMIB, LMM(I,J)) 6034. C**** Calculate energy flux EFLUX (J) from ice berg to ocean 6035. SEFLUX = 0. 6036. DO 130 L=1,LMAX 6037. BYMM(L) = BYDXYP(J) / MO(I,J,L) 6038. EFLUX(L) = -.125*ELHM*XMELT*DTS*SRM/LMIB 6039. GO = (G0M(I,J,L) + X*GXM(I,J,L) + Y*GYM(I,J,L)) * BYMM(L) 6040. SO = (S0M(I,J,L) + X*SXM(I,J,L) + Y*SYM(I,J,L)) * BYMM(L) 6041. TCO = TEMGS(GO,SO) 6042. IF(TCO.le.TCI) GO TO 130 6043. EFLUX(L) = EFLUX(L) - ELHM*XMELT*DTS*OMF(L)*SRM*(TCO-TCI) 6044. 130 SEFLUX = SEFLUX + EFLUX(L) 6045. C 130 continue 6046. C IF(SEFLUX.ge.0.) GO TO 110 6047. IF(BERGH(N)+ELHM*BERGM(N).gt.SEFLUX) GO TO 300 6048. C**** 6049. C**** No ice berg melting: add ice berg energy to ocean arrays 6050. C**** 6051. 200 BERGH(N) = BERGH(N) - SEFLUX 6052. AIJ(I,J,66) = AIJ(I,J,66) + SEFLUX*BYDXYP(J) 6053. IF(J.eq.JM) GO TO 220 6054. DO 210 L=1,LMAX 6055. G0M(I,J,L) = G0M(I,J,L) + EFLUX(L) 6056. GXM(I,J,L) = GXM(I,J,L) + 3.*X*EFLUX(L) 6057. 210 GYM(I,J,L) = GYM(I,J,L) + 3.*Y*EFLUX(L) 6058. GO TO 110 6059. C**** Add ice berg energy to ocean arrays at North Pole box 6060. 220 DO 230 L=1,LMAX 6061. 230 G0M(1,JM,L) = G0M(1,JM,L) + EFLUX(L)/IM 6062. GO TO 110 6063. C**** 6064. C**** Calculate ice berg melting 6065. C**** 6066. 300 IF(BERGH(N).lt.SEFLUX) GO TO 320 6067. C**** Entire ice berg will melt, recalculate flux to ocean 6068. DO 310 L=1,LMAX 6069. 310 EFLUX(L) = EFLUX(L)*(BERGH(N)/SEFLUX) 6070. SEFLUX = BERGH(N) 6071. C**** Calculate melt water flux MMELT (kg) from ice berg to ocean 6072. 320 SMMELT = BERGM(N) + (BERGH(N)-SEFLUX)/ELHM 6073. DO 330 L=1,LMAX 6074. 330 MMELT(L) = SMMELT*EFLUX(L)/SEFLUX 6075. C**** 6076. C**** Add melt water flux and energy flux to ocean arrays 6077. C**** 6078. IF(J.eq.JM) GO TO 420 6079. AIJ(I,J,65) = AIJ(I,J,65) + SMMELT*BYDXYP(J) 6080. AIJ(I,J,66) = AIJ(I,J,66) + SEFLUX*BYDXYP(J) 6081. DO 410 L=1,LMAX 6082. G0M(I,J,L) = G0M(I,J,L) + EFLUX(L) 6083. GXM(I,J,L) = GXM(I,J,L) + 3.*X*EFLUX(L) + MMELT(L)*BYMM(L)* 6084. * (GXM(I,J,L)*(.5-1.5*X*X)-3.*X*G0M(I,J,L)) 6085. GYM(I,J,L) = GYM(I,J,L) + 3.*Y*EFLUX(L) + MMELT(L)*BYMM(L)* 6086. * (GYM(I,J,L)*(.5-1.5*Y*Y)-3.*Y*G0M(I,J,L)) 6087. SXM(I,J,L) = SXM(I,J,L) + MMELT(L)*BYMM(L)* 6088. * (SXM(I,J,L)*(.5-1.5*X*X)-3.*X*S0M(I,J,L)) 6089. SYM(I,J,L) = SYM(I,J,L) + MMELT(L)*BYMM(L)* 6090. * (SYM(I,J,L)*(.5-1.5*Y*Y)-3.*Y*S0M(I,J,L)) 6091. 410 MO(I,J,L) = MO(I,J,L) + MMELT(L)*BYDXYP(J) 6092. GO TO 440 6093. C**** Add melt water flux and energy flux to ocean arrays at North Pole 6094. 420 AIJ(1,JM,65) = AIJ(1,JM,65) + SMMELT*BYDXYP(JM)/IM 6095. AIJ(1,JM,66) = AIJ(1,JM,66) + SEFLUX*BYDXYP(JM)/IM 6096. DO 430 L=1,LMAX 6097. G0M(1,JM,L) = G0M(1,JM,L) + EFLUX(L)/IM 6098. DO 430 I=1,IM 6099. 430 MO(I,JM,L) = MO(I,JM,L) + MMELT(L)*BYDXYP(JM)/IM 6100. 440 IF(BERGH(N).lt.SEFLUX) GO TO 460 6101. C**** Entire ice berg has melted 6102. NBERG = NBERG - 1 6103. DO 450 NN=N,NBERG 6104. BERGM(NN) = BERGM(NN+1) 6105. BERGH(NN) = BERGH(NN+1) 6106. BERGI(NN) = BERGI(NN+1) 6107. 450 BERGJ(NN) = BERGJ(NN+1) 6108. GO TO 120 6109. C**** Subtract melt water and energy from ice berg 6110. 460 BERGM(N) = BERGM(N) - SMMELT 6111. BERGH(N) = BERGH(N) - SEFLUX 6112. GO TO 110 6113. C**** 6114. C**** If there are too many ice bergs, melt the oldest 6115. C**** 6116. 500 IF(NBERG.le.128-4) RETURN 6117. I = BERGI(1) + 1. 6118. J = BERGJ(1) + 1. 6119. IF(J.eq.JM) I = 1 6120. LMAX = MIN (LMIB, LMM(I,J)) 6121. IF(LMAX.lt.LMIB) GO TO 520 6122. C**** LMAX = LMIB, calculate ice berg mass and energy for each layer 6123. DO 510 L=1,LMAX 6124. MMELT(L) = BERGM(1)*OMF(L) 6125. 510 EFLUX(L) = BERGH(1)*OMF(L) 6126. GO TO 540 6127. C**** LMAX < LMIB, calculate ice berg mass and energy for each layer 6128. 520 DO 530 L=1,LMAX 6129. MMELT(L) = BERGM(1)*DSIGO(L)/(SIGEO(LMAX)-SIGEO(0)) 6130. 530 EFLUX(L) = BERGH(1)*DSIGO(L)/(SIGEO(LMAX)-SIGEO(0)) 6131. C**** Add melt water flux and energy flux to ocean arrays 6132. 540 IF(J.eq.JM) GO TO 560 6133. AIJ(I,J,65) = AIJ(I,J,65) + BERGM(1)*BYDXYP(J) 6134. AIJ(I,J,66) = AIJ(I,J,66) + BERGH(1)*BYDXYP(J) 6135. DO 550 L=1,LMAX 6136. G0M(I,J,L) = G0M(I,J,L) + EFLUX(L) 6137. GXM(I,J,L) = GXM(I,J,L) + 3.*X*EFLUX(L) + MMELT(L)*BYMM(L)* 6138. * (GXM(I,J,L)*(.5-1.5*X*X)-3.*X*G0M(I,J,L)) 6139. GYM(I,J,L) = GYM(I,J,L) + 3.*Y*EFLUX(L) + MMELT(L)*BYMM(L)* 6140. * (GYM(I,J,L)*(.5-1.5*Y*Y)-3.*Y*G0M(I,J,L)) 6141. SXM(I,J,L) = SXM(I,J,L) + MMELT(L)*BYMM(L)* 6142. * (SXM(I,J,L)*(.5-1.5*X*X)-3.*X*S0M(I,J,L)) 6143. SYM(I,J,L) = SYM(I,J,L) + MMELT(L)*BYMM(L)* 6144. * (SYM(I,J,L)*(.5-1.5*Y*Y)-3.*Y*S0M(I,J,L)) 6145. 550 MO(I,J,L) = MO(I,J,L) + MMELT(L)*BYDXYP(J) 6146. GO TO 580 6147. C**** Add melt water flux and energy flux to ocean arrays at North Pole 6148. 560 AIJ(1,JM,65) = AIJ(1,JM,65) + BERGM(1)*BYDXYP(JM)/IM 6149. AIJ(1,JM,66) = AIJ(1,JM,66) + BERGH(1)*BYDXYP(JM)/IM 6150. DO 570 L=1,LMAX 6151. G0M(1,JM,L) = G0M(1,JM,L) + EFLUX(L)/IM 6152. DO 570 I=1,IM 6153. 570 MO(I,JM,L) = MO(I,JM,L) + MMELT(L)*BYDXYP(JM)/IM 6154. C**** Remove oldest ice berg from ice berg arrays 6155. 580 NBERG = NBERG - 1 6155.1 WRITE (6,*) 'Oldest ice berg melted.', 6155.2 * BERGM(1),BERGH(1),BERGI(1),BERGJ(1) 6156. DO 590 NN=1,NBERG 6157. BERGM(NN) = BERGM(NN+1) 6158. BERGH(NN) = BERGH(NN+1) 6159. BERGI(NN) = BERGI(NN+1) 6160. 590 BERGJ(NN) = BERGJ(NN+1) 6161. RETURN 6162. END 6500. 6501. SUBROUTINE CALVEI 6502. C**** 6503. C**** CALVEI causes glacial ice that has left the coast to calve into 6504. C**** the ocean, thus incresing the sea ice 6505. C**** 6506. INCLUDE 'C070.COM' 6507. PARAMETER (RHOI=910., MSICalve=2.*RHOI, SNOWCalve=45.5) 6508. LOGICAL*4 QSHELF 6509. COMMON /RIVRCB/ RATE(IM,JM), IFLOW(IM,JM),JFLOW(IM,JM), 6510. * KDIREC(IM,JM), QSHELF(IM,JM) 6511. C**** 6512. C**** GICECB MGI mass to be converted from glacial ice to sea ice (kg) 6513. C**** HGI heat to be converted from glacial ice to sea ice (J) 6514. C**** 6515. C**** SEAICB RSI horizontal ratio of sea ice to ocean (1) 6516. C**** MSI mass of sea ice (kg/m^2) 6517. C**** HSI enthalpy minus latent heat of sea ice (J/m^2) 6518. C**** 6519. DO 10 J=JM/12,JM/6+1 ! southern oceans only 6520. DO 10 I=1,IM 6521. IF(FOCEAN(I,J).le.0. .or. QSHELF(I,J) .or. MGI(I,J).le.0.) 6522. * GO TO 10 6523. DRSI = MGI(I,J)*BYDXYP(J)/MSICalve 6524. IF(DRSI .gt. 1.-RSI(I,J)) DRSI = 1.-RSI(I,J) 6525. DRSIY = - DRSI 6526. IF(DRSIY .lt. RSI(I,J)+DRSI-RSIY(I,J)-1.) 6527. * DRSIY = RSI(I,J)+DRSI-RSIY(I,J)-1. 6528. DRSIxMSI1C = DRSI*(ACE1I+SNOWCalve) 6529. DRSIxMSI2C = MGI(I,J)*BYDXYP(J) - DRSIxMSI1C 6530. DRSIxHSI1C = XSI1*DRSIxMSI1C*HGI(I,J,1) / MGI(I,J) 6531. DRSIxHSI2C = XSI2*DRSIxMSI1C*HGI(I,J,1) / MGI(I,J) 6532. DRSIxHSI3C = XSI3*DRSIxMSI2C*HGI(I,J,1) / MGI(I,J) 6533. DRSIxHSI4C = XSI4*DRSIxMSI2C*HGI(I,J,1) / MGI(I,J) 6534. MSI(I,J,1) = (RSI(I,J)*MSI(I,J,1) + DRSIxMSI1C) / (RSI(I,J)+DRSI) 6535. MSI(I,J,2) = (RSI(I,J)*MSI(I,J,2) + DRSIxMSI2C) / (RSI(I,J)+DRSI) 6536. HSI(I,J,1) = (RSI(I,J)*HSI(I,J,1) + DRSIxHSI1C) / (RSI(I,J)+DRSI) 6537. HSI(I,J,2) = (RSI(I,J)*HSI(I,J,2) + DRSIxHSI2C) / (RSI(I,J)+DRSI) 6538. HSI(I,J,3) = (RSI(I,J)*HSI(I,J,3) + DRSIxHSI3C) / (RSI(I,J)+DRSI) 6539. HSI(I,J,4) = (RSI(I,J)*HSI(I,J,4) + DRSIxHSI4C) / (RSI(I,J)+DRSI) 6540. RSI(I,J) = RSI(I,J) + DRSI 6541. RSIY(I,J) = RSIY(I,J) + DRSIY 6542. IF(RSI(I,J)+RSIX(I,J) .gt. 1.) RSIX(I,J) = 1.-RSI(I,J) 6543. IF(RSI(I,J)-RSIX(I,J) .gt. 1.) RSIX(I,J) = RSI(I,J)-1. 6544. MGI(I,J) = 0. 6545. HGI(I,J,1) = 0. 6546. 10 continue 6547. RETURN 6548. END 7000. 7001. SUBROUTINE CREAIB 7002. C**** 7003. C**** CREAIB creates an ice berg from the glacial ice coast line 7004. C**** accumulating arrays when those arrays contains too much ice 7005. C**** 7006. INCLUDE 'C070.COM' 7007. PARAMETER (MIBNEW=1.d13, GIMTSH=1.,GIMTNH=0.) 7008. REAL*8 GIMT(IM,JM) 7009. COMMON /RIVRCB/ RATE(IM,JM), IFLOW(IM,JM),JFLOW(IM,JM) 7010. C**** 7011. DATA IFIRST /1/ 7012. IF(IFIRST.eq.0) GO TO 100 7013. IFIRST = 0 7014. C**** 7015. C**** Read in constant glacial ice mass transport GIMT (10^12 kg/year) 7016. C**** Multiply GIMT times hemispheric constant factors 7017. C**** 7018. CALL READR4 (41,IM*JM,GIMT,GIMT) 7019. CLOSE (41) 7020. DO 10 I=1,IM*JM/2 7021. 10 GIMT(I,1) = GIMT(I,1) *GIMTSH 7022. C 10 GIMT(I,1+JM/2) = GIMT(I,1+JM/2)*GIMTNH 7023. C**** 7024. C**** Zero out GIMT unless it feeds into the ocean 7025. C**** 7026. DO 40 JU=2,JM/2 ! southern oceans only 7027. DO 40 IU=1,IM 7028. IF(FGICE(IU,JU).le.0) GO TO 30 7029. ID = IFLOW(IU,JU) 7030. JD = JFLOW(IU,JU) 7031. IF(FOCEAN(ID,JD).le.0) GO TO 30 7032. C**** GIMT(IU,JU) feeds into ocean, convert GIMT into unit of 7033. C**** (kg/source time step) 7034. GIMT(IU,JU) = GIMT(IU,JU)*1.d12*DTS/(SDAY*365.) 7035. GO TO 40 7036. C**** GIMT(IU,JU) does not feed into ocean, zero it out 7037. 30 GIMT(IU,JU) = 0. 7038. 40 continue 7039. C**** 7040. C**** Add glacial ice transport GIMT (kg) to glacial ice coast line 7041. C**** accumulating arrays MGI (kg) and HGI (J) 7042. C**** 7043. 100 DO 110 JU=2,JM/2 ! southern oceans only 7044. DO 110 IU=1,IM 7045. IF(GIMT(IU,JU).le.0) GO TO 110 7046. ID = IFLOW(IU,JU) 7047. JD = JFLOW(IU,JU) 7048. MGI(ID,JD) = MGI(ID,JD) + GIMT(IU,JU) 7049. HGI(ID,JD,1)= HGI(ID,JD,1)+ GIMT(IU,JU)*HGI(IU,JU,4)/(XGI4*ACE2GI) 7050. 110 continue 7051. C**** 7052. C**** If glacial ice coast line accumulating array MGI (kg) exceeds 7053. C**** MIBNEW, release it as an ice berg 7054. C**** 7055. DO 210 J=2,JM/2 ! southern oceans only 7056. DO 210 I=1,IM 7057. IF(MGI(I,J).lt.MIBNEW) GO TO 210 7058. IF(NBERG.ge.128) then 7059. WRITE (6,*) 'CREAIB: additional ice berg exceeds array size.' 7060. STOP 'CREAIB 200:' 7061. endif 7062. NBERG = NBERG+1 7063. BERGM(NBERG) = MGI(I,J) 7064. BERGH(NBERG) = HGI(I,J,1) 7065. BERGI(NBERG) = I - .5 7066. BERGJ(NBERG) = J - .5 7067. MGI(I,J) = 0. 7068. HGI(I,J,1) = 0. 7069. TIB = (BERGH(NBERG)/BERGM(NBERG)+ELHM)/SHCI 7070. WRITE (6,921) 'Ice berg create: N,IB,JB,T,M=', 7071. * NBERG,BERGI(NBERG),BERGJ(NBERG),TIB,BERGM(NBERG) 7072. 210 continue 7073. RETURN 7074. 921 FORMAT (1X,A29,I4,3F8.3,F16.0) 7075. END