1. C**** 2. C**** C089O.S Fortran source code for Ocean routines 98/11/12 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 OINPUT (ISTART) 12. C**** 13. C**** OINPUT calculates fixed arrays, and reads in the starting 14. C**** ocean prognostic variables when ISTART = 1 15. C**** 16. INCLUDE 'C070.COM' 17. PARAMETER (ZE1=12.,ZERAT=1.5, RHOI=910., MIBNEW=1.d13) 17.5 REAL*4 MO4,G0M4,GZM4,S0M4,SZM4 17.6 LOGICAL*4 QSHELF 18. CHARACTER*80 TITLE 19. COMMON /OFUNCB/ VGSP(-2:40,0:40,0:39),TGSP(-2:40,0:40,0:39), 20. * CGS(-2:40,0:40) ,HGSP(-2:40,0:40,0:39), 20.1 * AGSP(-2:40,0:40,0:39),BGSP(-2:40,0:40,0:39) 20.5 COMMON /RIVRCB/ RATE(IM,JM), IFLOW(IM,JM),JFLOW(IM,JM), 20.6 * KDIREC(IM,JM), QSHELF(IM,JM) 21. COMMON /WORK00/ MO4(IM,JM,LMO),G0M4(IM,JM,LMO),GZM4(IM,JM,LMO) 21.1 COMMON /WORK01/ S0M4(IM,JM,LMO),SZM4(IM,JM,LMO) 22. C**** 23. C**** Arrays needed each ocean model run 24. C**** 25. C**** Calculate ZE, SIGEO, DSIGO and SIGO 26. DO 110 L=0,LMO 27. 110 ZE(L) = ZE1*(ZERAT**L-1.)/(ZERAT-1.) 28. SIGEO(0) = 0. 29. DO 120 L=1,LMO 30. SIGEO(L) = ZE(L)/ZE(LMO) 31. DSIGO(L) = SIGEO(L) - SIGEO(L-1) 32. 120 SIGO(L) = (SIGEO(L) + SIGEO(L-1))*.5 33. C**** Read in table function for specific volume 34. READ (22) TITLE,VGSP 35. WRITE (6,*) 'VGSP read from unit 22: ',TITLE 36. READ (22) TITLE,TGSP 37. WRITE (6,*) 'TGSP read from unit 22: ',TITLE 38. READ (22) TITLE,CGS 39. WRITE (6,*) 'CGS read from unit 22: ',TITLE 39.1 READ (22) TITLE,HGSP 39.2 WRITE (6,*) 'HGSP read from unit 22: ',TITLE 39.3 READ (22) TITLE,AGSP 39.4 WRITE (6,*) 'AGSP read from unit 22: ',TITLE 39.5 READ (22) TITLE,BGSP 39.6 WRITE (6,*) 'BGSP read from unit 22: ',TITLE 40. CLOSE (22) 41. C**** Calculate LMM from ZE(L) and HOCEAN (stored in ZSOLID) 42. DO 170 I=1,IM*JM 43. L=0 44. IF(FOCEAN(I,1).le.0.) GO TO 160 45. DO 150 L=1,LMO-1 46. 150 IF(ZSOLID(I,1) .le. .5*(ZE(L)+ZE(L+1))) GO TO 160 47. C L=LMO 48. 160 LMM(I,1)=L 49. 170 ZSOLID(I,1) = ZATMO(I,1) - ZE(L) 50. C**** Calculate LMU 51. I=IM 52. DO 180 J=1,JM 53. DO 180 IP1=1,IM 54. LMU(I,J) = MIN(LMM(I,J),LMM(IP1,J)) 55. 180 I=IP1 56. C**** Calculate LMV 57. DO 190 J=1,JM-1 58. DO 190 I=1,IM 59. 190 LMV(I,J) = MIN(LMM(I,J),LMM(I,J+1)) 60. C**** 61. C**** Calculate distance of strait, distance between centers of 62. C**** ocean grid boxes through strait for pressure gradient force, 63. C**** and mass of water in the strait 64. C**** 64.5 DLON = TWOPI/IM 64.6 DLAT = TWOPI*NINT(360./(JM-1))/720. 65. FJEQ = (JM+1.)/2. 66. DO 230 N=1,NMST 67. DSLON = (IST(N,2)+.5*XST(N,2)-IST(N,1)-.5*XST(N,1))*DLON 68. DSLAT = (JST(N,2)+.5*YST(N,2)-JST(N,1)-.5*YST(N,1))*DLAT 69. SLAT = (JST(N,2)+.5*YST(N,2)+JST(N,1)+.5*YST(N,1)-2.*FJEQ)*DLAT/2 70. DIST(N) = RADIUS*SQRT((DSLON*COS(SLAT))**2 + DSLAT*DSLAT) 71. DFLON = .5*XST(N,1)*DLON 72. DFLAT = .5*YST(N,1)*DLAT 73. FLAT = (JST(N,1)+.25*YST(N,1)-FJEQ)*DLAT 74. DTLON = .5*XST(N,2)*DLON 75. DTLAT = .5*YST(N,2)*DLAT 76. TLAT = (JST(N,2)+.25*YST(N,2)-FJEQ)*DLAT 77. DISTPG(N) = DIST(N) + 78. + RADIUS*SQRT((DFLON*COS(FLAT))**2 + DFLAT*DFLAT) + 79. + RADIUS*SQRT((DTLON*COS(TLAT))**2 + DTLAT*DTLAT) 80. DO 210 L=1,LMST(N) 81. 210 MMST(L,N) = DIST(N)*WIST(N)*(ZE(L)-ZE(L-1))/.00097d0 82. DO 220 L=LMST(N)+1,LMO 83. 220 MMST(L,N) = 0. 84. 230 continue 85. C**** 86. IF(ISTART.ne.1) RETURN 87. C**** 88. C**** ISTART = 1: Initialize a run from ocean initial conditions 89. C**** 90. READ (20,ERR=820) TITLE,MO4,G0M4,GZM4,S0M4,SZM4 91. WRITE (6,*) 'OIC read from unit 10: ',TITLE 92. C**** Calculate layer mass from column mass and check for consistency 93. DO 313 J=1,JM 94. DO 313 I=1,IM 94.5 LMIJ=LMM(I,J) 95. DO 311 L=1,LMIJ 96. 311 MO(I,J,L) = MO4(I,J,L) 97. DO 312 L=LMIJ+1,LMO 98. 312 MO(I,J,L) = 0. 99. IF((LMM(I,J).GT.0).and.(ABS(MO(I,J,1)/ZE(1)-1000.).GT.50.)) 100. * WRITE (6,931) I,J,LMIJ,MO(I,J,1),ZE(1) 101. 313 continue 102. C**** Initialize velocity field to zero 103. DO 370 L=1,LMO 104. DO 320 I=1,IM*JM 105. UO(I,1,L) = 0. 106. 320 VO(I,1,L) = 0. 107. C**** Define mean value of mass, potential heat, and salinity at poles 108. J=JM 109. SM = 0. 110. SG0 = 0. 111. SGZ = 0. 112. SS0 = 0. 113. SSZ = 0. 114. DO 331 I=1,IM 115. SM = SM + MO(I,J,L) 116. SG0 = SG0 + G0M4(I,J,L) 117. SGZ = SGZ + GZM4(I,J,L) 118. SS0 = SS0 + S0M4(I,J,L) 119. 331 SSZ = SSZ + SZM4(I,J,L) 120. DO 332 I=1,IM 121. MO(I,J,L) = SM /IM 122. G0M4(I,J,L) = SG0/IM 123. GZM4(I,J,L) = SGZ/IM 124. S0M4(I,J,L) = SS0/IM 125. 332 SZM4(I,J,L) = SSZ/IM 126. C**** Define East-West horizontal gradients 127. DO 341 I=1,IM*JM 128. GXM(I,1,L) = 0. 129. GYM(I,1,L) = 0. 130. SXM(I,1,L) = 0. 131. 341 SYM(I,1,L) = 0. 132. IM1=IM-1 133. I=IM 134. DO 345 J=2,JM-1 135. DO 345 IP1=1,IM 136. IF(LMM(I ,J).LT.L) GO TO 344 137. IF(LMM(IM1,J).GE.L) GO TO 342 138. IF(LMM(IP1,J).LT.L) GO TO 344 139. GXM(I,J,L) = .5*(G0M4(IP1,J,L)-G0M4(I,J,L)) 140. SXM(I,J,L) = .5*(S0M4(IP1,J,L)-S0M4(I,J,L)) 141. GO TO 344 142. 342 IF(LMM(IP1,J).GE.L) GO TO 343 143. GXM(I,J,L) = .5*(G0M4(I,J,L)-G0M4(IM1,J,L)) 144. SXM(I,J,L) = .5*(S0M4(I,J,L)-S0M4(IM1,J,L)) 145. GO TO 344 146. 343 GXM(I,J,L) = .25*(G0M4(IP1,J,L)-G0M4(IM1,J,L)) 147. SXM(I,J,L) = .25*(S0M4(IP1,J,L)-S0M4(IM1,J,L)) 148. 344 IM1=I 149. 345 I=IP1 150. C**** Define North-South horizontal gradients 151. DO 354 J=2,JM-1 152. DO 354 I=1,IM 153. IF(LMM(I,J ).LT.L) GO TO 354 154. IF(LMM(I,J-1).GE.L) GO TO 352 155. IF(LMM(I,J+1).LT.L) GO TO 354 156. GYM(I,J,L) = .5*(G0M4(I,J+1,L)-G0M4(I,J,L)) 157. SYM(I,J,L) = .5*(S0M4(I,J+1,L)-S0M4(I,J,L)) 158. GO TO 354 159. 352 IF(LMM(I,J+1).GE.L) GO TO 353 160. GYM(I,J,L) = .5*(G0M4(I,J,L)-G0M4(I,J-1,L)) 161. SYM(I,J,L) = .5*(S0M4(I,J,L)-S0M4(I,J-1,L)) 162. GO TO 354 163. 353 GYM(I,J,L) = .25*(G0M4(I,J+1,L)-G0M4(I,J-1,L)) 164. SYM(I,J,L) = .25*(S0M4(I,J+1,L)-S0M4(I,J-1,L)) 165. 354 continue 166. C**** Multiply specific quantities by mass 167. DO 360 J=1,JM 168. DO 360 I=1,IM 169. G0M(I,J,L) = G0M4(I,J,L)*(MO(I,J,L)*DXYP(J)) 170. GXM(I,J,L) = GXM (I,J,L)*(MO(I,J,L)*DXYP(J)) 171. GYM(I,J,L) = GYM (I,J,L)*(MO(I,J,L)*DXYP(J)) 172. GZM(I,J,L) = GZM4(I,J,L)*(MO(I,J,L)*DXYP(J)) 173. S0M(I,J,L) = S0M4(I,J,L)*(MO(I,J,L)*DXYP(J)) 174. SXM(I,J,L) = SXM (I,J,L)*(MO(I,J,L)*DXYP(J)) 175. SYM(I,J,L) = SYM (I,J,L)*(MO(I,J,L)*DXYP(J)) 176. 360 SZM(I,J,L) = SZM4(I,J,L)*(MO(I,J,L)*DXYP(J)) 177. 370 continue 177.1 C**** 177.2 C**** Initialize gradients of ice cover and sea ice velocity 177.3 C**** 177.4 DO 410 I=1,IM*JM 177.5 RSIX(I,1) = 0. 177.6 PSI(I,1) = 0. 177.7 USI(I,1) = 0. 177.8 410 VSI(I,1) = 0. 177.9 DO 420 J=2,JM/2 178. DO 420 I=1,IM 178.1 RSIY(I,J) = -RSI(I,J) 178.2 IF(RSI(I,J).GT..5) RSIY(I,J) = RSI(I,J) - 1. 178.3 RSIY(I,JM+1-J) = RSI(I,JM+1-J) 178.4 420 IF(RSI(I,JM+1-J).GT..5) RSIY(I,JM+1-J) = 1. - RSI(I,JM+1-J) 178.9 C**** 179. C**** Initialize the mass flux, potential heat, and salinity in 180. C**** each strait 181. C**** 182. DO 530 N=1,NMST 183. I1=IST(N,1) 184. J1=JST(N,1) 185. I2=IST(N,2) 186. J2=JST(N,2) 187. DO 510 L=1,LMST(N) 188. MUST(L,N) = 0. 189. G01 = (G0M(I1,J1,L)+XST(N,1)*GXM(I1,J1,L)+YST(N,1)*GYM(I1,J1,L))/ 190. / (MO(I1,J1,L)*DXYP(J1)) 191. GZ1 = GZM(I1,J1,L)/(MO(I1,J1,L)*DXYP(J1)) 192. G02 = (G0M(I2,J2,L)+XST(N,2)*GXM(I2,J2,L)+YST(N,2)*GYM(I2,J2,L))/ 193. / (MO(I2,J2,L)*DXYP(J2)) 194. GZ2 = GZM(I2,J2,L)/(MO(I2,J2,L)*DXYP(J2)) 195. G0MST(L,N) = (G01+G02)*MMST(L,N)*.5 196. GXMST(L,N) = (G02-G01)*MMST(L,N)*.5 197. GZMST(L,N) = (GZ1+GZ2)*MMST(L,N)*.5 198. S01 = (S0M(I1,J1,L)+XST(N,1)*SXM(I1,J1,L)+YST(N,1)*SYM(I1,J1,L))/ 199. / (MO(I1,J1,L)*DXYP(J1)) 200. SZ1 = SZM(I1,J1,L)/(MO(I1,J1,L)*DXYP(J1)) 201. S02 = (S0M(I2,J2,L)+XST(N,2)*SXM(I2,J2,L)+YST(N,2)*SYM(I2,J2,L))/ 202. / (MO(I2,J2,L)*DXYP(J2)) 203. SZ2 = SZM(I2,J2,L)/(MO(I2,J2,L)*DXYP(J2)) 204. S0MST(L,N) = (S01+S02)*MMST(L,N)*.5 205. SXMST(L,N) = (S02-S01)*MMST(L,N)*.5 206. 510 SZMST(L,N) = (SZ1+SZ2)*MMST(L,N)*.5 207. DO 520 L=LMST(N)+1,LMO 208. MUST(L,N) = 0. 209. G0MST(L,N) = 0. 210. GXMST(L,N) = 0. 211. GZMST(L,N) = 0. 212. S0MST(L,N) = 0. 213. SXMST(L,N) = 0. 214. 520 SZMST(L,N) = 0. 215. 530 continue 216. C**** Initialize sea ice in straits 217. DO 540 N=1,NMST 218. RSIST(N) = 0. 219. RSIXST(N) = 0. 220. DO 540 K=1,2+LMSI 221. 540 MSIST(K,N) = 0. 222. C**** 223. C**** Initialize glacial ice coast line accumulating arrays 224. C**** 225. DO 610 J=2,JM 226. DO 610 I=1,IM 227. IF(FGICE(I,J).le.0.) MGI(I,J) = 0. 227.1 IF(FGICE(I,J).le.0.) HGI(I,J,1) = 0. 228. IF(QSHELF(I,J)) MGI(I,J) = MIBNEW*(.25*MOD(I+J,4)+.125) 229. IF(QSHELF(I,J)) HGI(I,J,1) = - ELHM*MGI(I,J) 230. 610 continue 231. RETURN 232. C**** Terminate because of improper start up 233. 820 STOP 'OINPUT 820: Error reading ocean IC on unit 10.' 234. C**** 235. 931 FORMAT ('0Inconsistency between LMM and M:',3I4,2F10.1) 236. END 1000. 1001. SUBROUTINE OVTOM (MM,UM,VM) 1002. C**** 1003. C**** OVTOM converts density and velocities in concentration units 1004. C**** to mass and momentum in mass units 1005. C**** Input: M (kg/m**2), U (m/s), V (m/s) 1006. C**** Output: MM (kg), UM (kg*m/s), VM (kg*m/s) 1007. C**** 1008. IMPLICIT REAL*8 (A-H,M,O-Z) 1009. PARAMETER (IM=72,JM=46,LMO=13) 1010. REAL*8 MM(IM,JM,LMO),UM(IM,JM,LMO),VM(IM,JM,LMO) 1011. COMMON /GEOMCB/ DXYP(JM) 1012. COMMON /OCENCB/ MO(IM,JM,LMO),UO(IM,JM,LMO),VO(IM,JM,LMO) 1013. DO 70 L=1,LMO 1014. C**** Convert density to mass 1015. DO 10 J=1,JM 1016. DO 10 I=1,IM 1017. 10 MM(I,J,L) = MO(I,J,L)*DXYP(J) 1018. C**** Convert U velocity to U momentum 1019. I=IM 1020. DO 40 J=2,JM-1 1021. DO 40 IP1=1,IM 1022. UM(I, J,L) = UO(I, J,L)*(MM(I,J,L)+MM(IP1,J,L))*.5 1023. 40 I=IP1 1024. C UM(IM,1,L) = UO(IM,1,L)*MM(IM,1,L)*IM 1025. UM(1,JM,L) = UO(1,JM,L)*MM(1,JM,L)*IM 1026. C**** Convert V velocity to V momentum 1027. DO 60 I=1,IM 1028. C VM(I, 1 ,L) = VO(I, 1 ,L)*(MM(1, 1,L)+.5*MM(I, 2 ,L)) 1029. VM(I,JM-1,L) = VO(I,JM-1,L)*(MM(1,JM,L)+.5*MM(I,JM-1,L)) 1030. DO 60 J=2,JM-2 1031. 60 VM(I,J,L) = VO(I,J,L)*(MM(I,J,L)+MM(I,J+1,L))*.5 1032. 70 continue 1033. RETURN 1034. END 1100. 1101. SUBROUTINE OMTOV (MM,UM,VM) 1102. C**** 1103. C**** OMTOV converts mass and momentum in mass units 1104. C**** to density and velocity in concentration units 1105. C**** Input: MM (kg), UM (kg*m/s), VM (kg*m/s) 1106. C**** Output: M (kg/m**2), U (m/s), V (m/s) 1107. C**** 1108. IMPLICIT REAL*8 (A-H,M,O-Z) 1109. PARAMETER (IM=72,JM=46,LMO=13) 1110. REAL*8 MM(IM,JM,LMO),UM(IM,JM,LMO),VM(IM,JM,LMO) 1111. COMMON /FIXDCB/ FANDZ(IM,JM,8),LMM(IM,JM),LMU(IM,JM),LMV(IM,JM) 1112. COMMON /GEOMCB/ DXYP(JM) 1113. COMMON /OCENCB/ MO(IM,JM,LMO),UO(IM,JM,LMO),VO(IM,JM,LMO) 1114. C**** Convert mass to density 1115. DO 10 J=1,JM 1116. DO 10 I=1,IM 1117. DO 10 L=1,LMM(I,J) 1118. 10 MO(I,J,L) = MM(I,J,L)/DXYP(J) 1119. C**** Convert U momentum to U velocity 1120. I=IM 1121. DO 41 J=2,JM-1 1122. DO 41 IP1=1,IM 1123. DO 40 L=1,LMU(I,J) 1124. 40 UO(I, J,L) = UM(I, J,L)*2./(MM(I,J,L)+MM(IP1,J,L)) 1125. 41 I=IP1 1126. C UO(IM,1,L) = UM(IM,1,L)/(MM(IM,1,L)*IM) 1127. DO 50 L=1,LMU(1,JM) 1128. UO(1,JM,L) = UM(1,JM,L)/(MM(1,JM,L)*IM) 1129. DO 50 I=2,IM 1130. 50 UO(I,JM,L) = UO(1,JM,L) 1131. C**** Convert V momentum to V velocity 1132. DO 70 I=1,IM 1133. DO 60 J=2,JM-2 1134. DO 60 L=1,LMV(I,J) 1135. 60 VO(I,J,L) = VM(I,J,L)*2./(MM(I,J,L)+MM(I,J+1,L)) 1136. C VO(I, 1 ,L) = VM(I, 1 ,L)/(MM(I, 1,L)+.5*MM(I, 2 ,L)) 1137. DO 70 L=1,LMV(I,JM-1) 1138. 70 VO(I,JM-1,L) = VM(I,JM-1,L)/(MM(I,JM,L)+.5*MM(I,JM-1,L)) 1139. RETURN 1140. END 1200. 1201. SUBROUTINE OFLUX (NS,MM,SMU) 1202. C**** 1203. C**** AFLUX calculates the fluid fluxes 1204. C**** Input: M (kg/m**2), U (m/s), V (m/s) 1205. C**** Output: MU (kg/s) = DY * U * M 1206. C**** MV (kg/s) = DX * V * M 1207. C**** MW (kg/s) = MW(L-1) + CONV-SCONV*DSIGO + (MM-SMM*DSIGO) 1208. C**** CONV (kg/s) = MU-MU+MV-MV 1209. C**** 1210. IMPLICIT REAL*8 (A-H,M,O-Z) 1211. PARAMETER (IM=72,JM=46,LMO=13, CHI=.125) 1212. REAL*8 MM(IM,JM,LMO),SMU(IM,JM,LMO,3) 1214. COMMON /PARMCB/ IPARM(100), DTS,DTA,DTO 1215. COMMON /FIXDCB/ FANDZ(IM,JM,8),LMM(IM,JM),LMU(IM,JM),LMV(IM,JM) 1216. COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM) 1217. COMMON /LAYOCB/ DSIGO(LMO),SIGEO(0:LMO) 1218. COMMON /OCENCB/ MO(IM,JM,LMO),UO(IM,JM,LMO),VO(IM,JM,LMO) 1219. COMMON /FLUXCB/ MU(IM,JM,LMO),MV(IM,JM,LMO),MW(IM,JM,LMO-1) 1220. COMMON /WORK03/ CONV(IM,JM,LMO),DMVS(IM),DMVN(IM) 1221. DATA SKIP /Z'FFEFFFFFFFFFFFFF'/ 1222. C**** 1223. C**** Compute fluid fluxes for the C grid 1224. C**** 1225. C**** Smooth the West-East velocity near the poles 1226. DO 110 L=1,LMO 1227. DO 110 I=IM+1,IM*(JM-1) 1228. 110 MU(I,1,L) = UO(I,1,L) 1229. CALL OPFIL (MU) 1230. C**** Compute MU, the West-East mass flux, at non-polar points 1231. DO 220 L=1,LMO 1232. I=IM 1233. DO 120 J=2,JM-1 1234. DO 120 IP1=1,IM 1235. MU(I,J,L) = .5*DYP(J)*MU(I,J,L)*(MO(I,J,L)+MO(IP1,J,L)) 1236. 120 I=IP1 1237. C**** Compute MV, the South-North mass flux 1238. DO 130 J=1,JM-1 1239. DO 130 I=1,IM 1240. 130 MV(I,J,L) = .5*DXV(J)*VO(I,J,L)*(MO(I,J,L)+MO(I,J+1,L)) 1241. C**** Compute MU*2/(1-2*CHI) at the poles 1242. MUS = DYP( 1)*UO(1, 1,L)*MO(1, 1,L) 1243. MUN = DYP(JM)*UO(1,JM,L)*MO(1,JM,L) 1244. MVS = 0. 1245. MVN = 0. 1246. DO 140 I=1,IM 1247. MVS = MVS + MV(I, 1 ,L) 1248. 140 MVN = MVN + MV(I,JM-1,L) 1249. MVS = MVS/IM 1250. MVN = MVN/IM 1251. DMVS(1) = 0. 1252. DMVN(1) = 0. 1253. DO 150 I=2,IM 1254. DMVS(I) = DMVS(I-1) + (MV(I, 1 ,L)-MVS) 1255. 150 DMVN(I) = DMVN(I-1) + (MV(I,JM-1,L)-MVN) 1256. ADMVS = 0. 1257. ADMVN = 0. 1258. DO 160 I=1,IM 1259. ADMVS = ADMVS + DMVS(I) 1260. 160 ADMVN = ADMVN + DMVN(I) 1261. ADMVS = ADMVS/IM 1262. ADMVN = ADMVN/IM 1263. DO 170 I=1,IM 1264. MU(I, 1,L) = (ADMVS-DMVS(I) + MUS)*2./(1.-2.*CHI) 1265. 170 MU(I,JM,L) = (DMVN(I)-ADMVN + MUN)*2./(1.-2.*CHI) 1266. C**** 1267. C**** Compute horizontal fluid convergence 1268. C**** 1269. IM1=IM 1270. DO 210 J=2,JM-1 1271. DO 210 I=1,IM 1272. CONV(I,J,L) = MU(IM1,J,L)-MU(I,J,L) + (MV(I,J-1,L)-MV(I,J,L)) 1273. 210 IM1=I 1274. CONV(IM,1,L) = -MVS 1275. CONV(1,JM,L) = MVN 1276. 220 continue 1277. C**** 1278. C**** Compute vertically integrated column convergence and mass 1279. C**** 1280. DO 440 I=IM,IM*(JM-1)+1 1281. LMIJ=1 1282. IF(LMM(I,1).le.1) GO TO 420 1283. LMIJ=LMM(I,1) 1284. SCONV = 0. 1285. SMM = 0. 1286. DO 310 L=1,LMIJ 1287. SCONV = SCONV + CONV(I,1,L) 1288. 310 SMM = SMM + MM(I,1,L) 1289. SCONV = SCONV/SIGEO(LMIJ) 1290. SMM = SMM /SIGEO(LMIJ) 1291. C**** 1292. C**** Compute MW, the downward fluid flux 1293. C**** 1294. MW(I,1,1) = CONV(I,1,1)-SCONV*DSIGO(1) + 1295. + ( MM(I,1,1)- SMM*DSIGO(1))/(NS*DTO) 1296. DO 410 L=2,LMIJ-1 1297. 410 MW(I,1,L) = CONV(I,1,L)-SCONV*DSIGO(L) + MW(I,1,L-1) + 1298. + ( MM(I,1,L)- SMM*DSIGO(L))/(NS*DTO) 1299. 420 DO 430 L=LMIJ,LMO-1 1300. 430 MW(I,1,L) = 0. 1301. 440 continue 1302. C**** 1303. C**** Sum mass fluxes to be used for advection of G and S 1304. C**** 1305. IF(SMU(1,1,1,1).eq.SKIP) RETURN 1306. DO 520 L=1,LMO 1307. DO 510 I=IM+1,IM*(JM-1) 1308. 510 SMU(I,1,L,1) = SMU(I,1,L,1) + MU(I,1,L) 1309. DO 520 I=1,IM*(JM-1) 1310. 520 SMU(I,1,L,2) = SMU(I,1,L,2) + MV(I,1,L) 1311. DO 530 L=1,LMO-1 1312. DO 530 I=IM,IM*(JM-1)+1 1313. 530 SMU(I,1,L,3) = SMU(I,1,L,3) + MW(I,1,L) 1314. RETURN 1315. END 1400. 1401. SUBROUTINE OPFIL (X) 1402. C**** 1403. C**** OPFIL smoothes X in zonal direction by reducing coefficients 1404. C**** of its Fourier series for high wave numbers near the poles. 1405. C**** INDM for polar ocean filter set to work with AVR4X5LD.Z12. 1405.1 C**** 1405.5 IMPLICIT REAL*8 (A-H,M,O-Z) 1406. PARAMETER (IM=72,JM=46,LMO=13, TWOPI=6.283185307179586477d0, 1406.5 * DLON=TWOPI/IM, NMAX=IM/2, NSEGM=7,INDM=90365) 1407. INTEGER*2 NSEG(LMO,4:25),IMIN(NSEGM,LMO,4:25),ILEN(NSEGM,LMO,4:25) 1408. INTEGER*4 INDEX(IM,2:13) 1409. REAL*4 REDUCO(INDM) 1410. REAL*8 X(IM,JM,LMO) 1411. CHARACTER*80 TITLE 1412. COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM) 1412.5 COMMON /AVRCOM/ AVCOS(IM,NMAX),AVSIN(IM,NMAX),BYSN(NMAX), 1413. * AN(0:NMAX),BN(0:NMAX), Y(IM) 1413.5 C**** 1414. DATA IFIRST /1/ 1415. IF(IFIRST.eq.0) GO TO 100 1416. IFIRST=0 1417. IF(JM.ne.46) STOP 'JM not equal to 46 in OPFIL.' 1418. C**** Calculate cos(TWOPI*N*I/IM) and sin(TWOPI*N*I/IM) 1420. DO 10 I=1,NMAX 1421. AVCOS(I,1) = COS(DLON*I) 1422. AVSIN(I,1) = SIN(DLON*I) 1423. AVCOS(I+NMAX,1) = -AVCOS(I,1) 1424. 10 AVSIN(I+NMAX,1) = -AVSIN(I,1) 1425. DO 20 N=2,NMAX 1426. DO 20 I=1,IM 1427. IN = I*N-((I*N-1)/IM)*IM 1428. AVCOS(I,N) = AVCOS(IN,1) 1429. 20 AVSIN(I,N) = AVSIN(IN,1) 1430. C**** Read in reduction contribution matrices from disk 1431. READ (23) TITLE,NSEG,IMIN,ILEN,INDEX,REDUCO 1432. CLOSE (23) 1433. WRITE (6,*) 'Read on unit 17: ',TITLE 1434. C**** 1435. C**** Loop over J and L. J = latitude, JA = absolute latitude. 1436. C**** 1437. 100 DO 330 JX=4,25 1438. J =JX 1439. JA=JX 1440. IF(JX.GT.13) J =20+JX 1441. IF(JX.GT.13) JA=27-JX 1442. DO 330 L=1,LMO 1443. IF(ILEN(1,L,JX).GE.IM) GO TO 300 1444. C**** 1445. C**** Land boxes exist at this latitude and layer, loop over ocean 1446. C**** segments. 1447. C**** 1448. DO 270 NS=1,NSEG(L,JX) 1449. I0 = IMIN(NS,L,JX) 1450. IL = ILEN(NS,L,JX) 1451. IND= INDEX(IL,JA) 1452. IF(IL+I0.gt.IM) GO TO 200 1453. C**** Ocean segment does not wrap around the IDL. 1454. C**** Copy X to temporary array and filter X in place. 1455. DO 110 I=1+I0,IL+I0 1456. 110 Y(I-I0) = X(I,J,L) 1457. DO 140 I=1+I0,IL+I0 1458. REDUC = 0. 1459. DO 130 K=1,IL 1460. IND=IND+1 1461. 130 REDUC = REDUC + REDUCO(IND)*Y(K) 1462. 140 X(I,J,L) = X(I,J,L) - REDUC 1463. GO TO 270 1464. C**** Ocean segment wraps around the IDL. 1465. C**** Copy X to temporary array and filter X in place. 1466. 200 DO 210 I=1+I0,IM 1467. 210 Y(I-I0) = X(I,J,L) 1468. DO 220 I=1,IL+I0-IM 1469. 220 Y(I-I0+IM) = X(I,J,L) 1470. DO 240 I=1+I0,IM 1471. REDUC = 0. 1472. DO 230 K=1,IL 1473. IND=IND+1 1474. 230 REDUC = REDUC + REDUCO(IND)*Y(K) 1475. 240 X(I,J,L) = X(I,J,L) - REDUC 1476. DO 260 I=1,IL+I0-IM 1477. REDUC = 0. 1478. DO 250 K=1,IL 1479. IND=IND+1 1480. 250 REDUC = REDUC + REDUCO(IND)*Y(K) 1481. 260 X(I,J,L) = X(I,J,L) - REDUC 1482. 270 continue 1483. GO TO 330 1484. C**** 1485. C**** No land boxes at this latitude and layer, 1486. C**** perform standard polar filter 1487. C**** 1488. 300 CALL FFT (X(1,J,L),AN,BN) 1489. DRATM = NMAX*DXP(J)/DYP(3) 1490. DO 320 N=NMAX,1,-1 1491. SM = 1. - DRATM/N 1492. IF(SM.le.0.) GO TO 330 1493. DO 320 I=1,IM 1494. 320 X(I,J,L) = X(I,J,L) - SM*(AN(N)*AVCOS(I,N)+BN(N)*AVSIN(I,N)) 1495. 330 continue 1496. RETURN 1497. END 1500. 1501. SUBROUTINE OADVM (MM2,MM0,DT) 1502. C**** 1503. C**** OADVM calculates the updated column mass 1504. C**** Input: MM0 (kg), CONV (kg/s), MW (kg/s), DT (s) 1505. C**** Output: MM2 (kg) = MM0 + DT*DM*DSIGO 1506. C**** 1506.5 IMPLICIT REAL*8 (A-H,M,O-Z) 1507. PARAMETER (IM=72,JM=46,LMO=13) 1509. REAL*8 MM2(IM,JM,LMO),MM0(IM,JM,LMO) 1510. COMMON /FLUXCB/ MU(IM,JM,LMO),MV(IM,JM,LMO),MW(IM,JM,LMO-1) 1511. COMMON /FIXDCB/ FANDZ(IM,JM,8),LMM(IM,JM) 1512. COMMON /WORK03/ CONV(IM,JM,LMO) 1513. C**** 1514. C**** Compute the new mass MM2 1515. C**** 1516. DO 40 I=IM,IM*(JM-1)+1 1517. LMIJ=LMM(I,1) 1518. IF(LMIJ-1) 40,10,20 1519. 10 MM2(I,1,1) = MM0(I,1,1) + DT*CONV(I,1,1) 1520. GO TO 40 1521. 20 MM2(I,1,1) = MM0(I,1,1) + DT*(CONV(I,1,1)-MW(I,1,1)) 1522. DO 30 L=2,LMIJ-1 1523. 30 MM2(I,1,L) = MM0(I,1,L) + DT*(CONV(I,1,L)+(MW(I,1,L-1)-MW(I,1,L))) 1524. MM2(I,1,LMIJ) = MM0(I,1,LMIJ) + DT*(CONV(I,1,LMIJ)+MW(I,1,LMIJ-1)) 1527. 40 continue 1528. C**** Fill in values at the poles 1529. DO 80 L=1,LMO 1530. DO 80 I=1,IM 1531. C MM2(I, 1,L) = MM2(IM,1,L) 1532. 80 MM2(I,JM,L) = MM2(1,JM,L) 1533. RETURN 1534. END 1600. 1601. SUBROUTINE OADVV (UM2,VM2,UM0,VM0,DT1) 1602. C**** 1603. C**** OADVV advects oceanic momentum (with coriolis force) 1604. C**** Input: MO (kg/m**2), UO (m/s), VO (m/s) = from odd solution 1605. C**** MU (kg/s), MV (kg/s), MW (kg/s) = fluid fluxes 1606. C**** Output: UM2 (kg*m/s) = UM0 + DT*(MU*U-MU*U+MV*U-MV*U+M*CM*V) 1607. C**** VM2 (kg*m/s) = VM0 + DT*(MU*V-MU*V+MV*V-MV*V-M*CM*U) 1608. C**** 1608.5 IMPLICIT REAL*8 (A-H,M,O-Z) 1609. PARAMETER (IM=72,JM=46,LMO=13, TWOPI=6.283185307179586477d0, 1609.5 * SDAY=86400., RADIUS=6375000., OMEGA=TWOPI*366./(365.*SDAY), 1609.6 * CHI=.125) 1611. REAL*8 UM2(IM,JM,LMO),VM2(IM,JM,LMO),UM0(IM,JM,LMO), 1612. * VM0(IM,JM,LMO), DCOSP(JM),TANP(JM) 1613. COMMON /FLUXCB/ MU(IM,JM,LMO),MV(IM,JM,LMO),MW(IM,JM,LMO-1) 1614. COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM), 1615. * COSV(0:JM),COSP(JM) 1616. COMMON /OCENCB/ MO(IM,JM,LMO),UO(IM,JM,LMO),VO(IM,JM,LMO) 1617. COMMON /WORK03/ DUM(IM,JM,LMO),DVM(IM,JM,LMO), 1618. * GX(0:JM),GXUY(0:JM),DGDUDN(0:JM),DGDUUP(0:JM), 1619. * UDCOSY(IM),UYUTAY(IM) 1620. LOGICAL*4 QFIRST/.TRUE./ 1621. C**** 1622. IF(.not.QFIRST) GO TO 20 1623. QFIRST = .FALSE. 1624. DO 10 J=1,JM 1625. DCOSP(J) = COSV(J-1)-COSV(J) 1626. 10 TANP(J) = DCOSP(J)/COSP(J) 1627. C**** 1628. 20 DT2 = DT1/2. 1629. DT4 = DT1/4. 1630. DTC2 = DT1*CHI/2. 1631. DTD4 = DT1*(1.-2.*CHI)/4. 1632. C**** Zero out momentum changes 1633. DO 30 I=1,IM*JM*LMO*2 1634. 30 DUM(I,1,1) = 0. 1635. C**** 1636. C**** Horizontal advection of momentum 1637. C**** 1638. DO 480 L=1,LMO 1639. IM1=IM-1 1640. I=IM 1641. DO 350 IP1=1,IM 1642. C**** Contribution of West-East flux to U wind 1643. DO 310 J=2,JM-1 1644. UMUX = DT4*(MU(IM1,J,L)+MU(I,J,L))*(UO(IM1,J,L)+UO(I,J,L)) 1645. DUM(IM1,J,L) = DUM(IM1,J,L) - UMUX 1646. 310 DUM(I ,J,L) = DUM(I ,J,L) + UMUX 1647. C**** Contribution of South-North and corner fluxes to U wind 1648. DO 320 J=1,JM-1 1649. UMVX = DTD4*(MV(I,J,L)+MV(IP1,J,L))*(UO(I,J,L)+UO(I,J+1,L)) 1650. UMVNE = DTC2* MV(I,J,L)*(UO(IM1,J,L)+UO(I,J+1,L)) 1651. UMVNW = DTC2* MV(I,J,L)*(UO(I,J,L)+UO(IM1,J+1,L)) 1652. DUM(IM1,J ,L) = DUM(IM1,J ,L) - UMVNE 1653. DUM(I ,J ,L) = DUM(I ,J ,L) - (UMVNW+UMVX) 1654. DUM(IM1,J+1,L) = DUM(IM1,J+1,L) + UMVNW 1655. 320 DUM(I ,J+1,L) = DUM(I ,J+1,L) + (UMVNE+UMVX) 1656. C**** Contribution of West-East flux to V wind 1657. DO 330 J=1,JM-1 1658. VMUY = DTD4*(MU(I,J,L)+MU(I,J+1,L))*(VO(I,J,L)+VO(IP1,J,L)) 1659. DVM(I ,J,L) = DVM(I ,J,L) - VMUY 1660. 330 DVM(IP1,J,L) = DVM(IP1,J,L) + VMUY 1661. C**** Contribution of South-North and corner fluxes to V wind 1662. DO 340 J=2,JM-1 1663. VMVY = DT4*(MV(I,J-1,L)+MV(I,J,L))*(VO(I,J-1,L)+VO(I,J,L)) 1664. VMUNE = DTC2*MU(I,J,L)*(VO(I,J-1,L)+VO(IP1,J,L)) 1665. VMUSE = DTC2*MU(I,J,L)*(VO(I,J,L)+VO(IP1,J-1,L)) 1666. DVM(I ,J-1,L) = DVM(I ,J-1,L) - (VMUNE+VMVY) 1667. DVM(IP1,J-1,L) = DVM(IP1,J-1,L) + VMUSE 1668. DVM(I ,J ,L) = DVM(I ,J ,L) - (VMUSE-VMVY) 1669. 340 DVM(IP1,J ,L) = DVM(IP1,J ,L) + VMUNE 1670. IM1=I 1671. 350 I=IP1 1672. C**** 1673. C**** Coriolis force and metric term 1674. C**** 1675. C**** U component 1676. GX( 0) = 0. 1677. GX(JM) = 0. 1678. GXUY( 0) = 0. 1679. GXUY(JM) = 0. 1680. DGDUDN( 0) = 0. 1681. DGDUDN(JM) = 0. 1682. DGDUUP( 0) = 0. 1683. DGDUUP(JM) = 0. 1684. IM1=IM-1 1685. I=IM 1686. DO 430 IP1=1,IM 1687. DO 410 J=1,JM-1 1688. GX(J) = MV(I,J,L)+MV(IP1,J,L) 1689. GXUY(J) = (MV(I,J,L)+MV(IP1,J,L))*(UO(I,J,L)+UO(I,J+1,L)) 1690. DGDUDN(J) = MV(I ,J,L)*(UO(IM1,J ,L)-UO(I ,J ,L)) - 1691. * MV(IP1,J,L)*(UO(I ,J ,L)-UO(IP1,J ,L)) 1692. 410 DGDUUP(J) = MV(I ,J,L)*(UO(IM1,J+1,L)-UO(I ,J+1,L)) - 1693. * MV(IP1,J,L)*(UO(I ,J+1,L)-UO(IP1,J+1,L)) 1694. DO 420 J=1,JM 1695. 420 DUM(I,J,L) = DUM(I,J,L) + 1696. * DT2*(OMEGA*RADIUS*(GX(J-1)+GX(J))*DCOSP(J) + TANP(J)* 1697. * (.25*(GXUY(J-1)+GXUY(J)) + .5*CHI*(DGDUDN(J-1)+DGDUUP(J)))) 1698. IM1=I 1699. 430 I=IP1 1700. C**** V component 1701. IM1=IM 1702. DO 450 J=1,JM-1 1703. DO 440 I=1,IM 1704. UDCOSY(I) = UO(I,J,L)*DCOSP(J) + UO(I,J+1,L)*DCOSP(J+1) 1705. 440 UYUTAY(I) = (UO(I,J,L)* TANP(J) + UO(I,J+1,L)* TANP(J+1))* 1706. * (UO(I,J,L)+UO(I,J+1,L)) 1707. DO 450 I=1,IM 1708. DVM(I,J,L) = DVM(I,J,L) - DT4*DXV(J)*(MO(I,J,L)+MO(I,J+1,L))* 1709. * (OMEGA*RADIUS*(UDCOSY(IM1)+UDCOSY(I)) + 1710. * .25*(UYUTAY(IM1)+UYUTAY(I)) - .5*CHI*(TANP(J)+TANP(J+1))* 1711. * (UO(IM1,J,L)-UO(I,J,L))*(UO(IM1,J+1,L)-UO(I,J+1,L))) 1712. 450 IM1=I 1713. 480 continue 1714. C**** 1715. C**** Vertical advection of momentum 1716. C**** 1717. DO 560 L=1,LMO-1 1718. C**** U component 1719. I=IM 1720. DO 530 J=2,JM-1 1721. DO 530 IP1=1,IM 1722. FLUX = DT4*(MW(I,J,L)+MW(IP1,J,L))*(UO(I,J,L)+UO(I,J,L+1)) 1723. DUM(I,J,L) = DUM(I,J,L) - FLUX 1724. DUM(I,J,L+1) = DUM(I,J,L+1) + FLUX 1725. 530 I=IP1 1726. C FLUX = DT2*MW(IM,1,L)*IM*(UO(IM,1,L)+UO(IM,1,L+1)) 1727. C DUM(IM,1,L) = DUM(IM,1,L) - FLUX 1728. C DUM(IM,1,L+1) = DUM(IM,1,L+1) + FLUX 1729. FLUX = DT2*MW(1,JM,L)*IM*(UO(1,JM,L)+UO(1,JM,L+1)) 1730. DUM(1,JM,L) = DUM(1,JM,L) - FLUX 1731. DUM(1,JM,L+1) = DUM(1,JM,L+1) + FLUX 1732. C**** V component 1733. DO 550 I=1,IM 1734. FLUX = DT4*(MW(I,JM-1,L)+2.*MW(1,JM,L))* 1734.1 * (VO(I,JM-1,L)+VO(I,JM-1,L+1)) 1735. DVM(I,JM-1,L) = DVM(I,JM-1,L) - FLUX 1736. DVM(I,JM-1,L+1) = DVM(I,JM-1,L+1) + FLUX 1737. DO 550 J=2,JM-2 1738. FLUX = DT4*(MW(I,J,L)+MW(I,J+1,L))*(VO(I,J,L)+VO(I,J,L+1)) 1739. DVM(I,J,L) = DVM(I,J,L) - FLUX 1740. 550 DVM(I,J,L+1) = DVM(I,J,L+1) + FLUX 1741. 560 continue 1742. C**** 1743. C**** Add changes to momentum 1744. C**** 1745. DO 640 L=1,LMO 1746. C**** U component 1747. DO 610 I=IM+1,IM*(JM-1) 1748. 610 UM2(I,1,L) = UM0(I,1,L) + DUM(I,1,L) 1749. DUMS = 0. 1750. DUMN = 0. 1751. DO 620 I=1,IM 1752. DUMS = DUMS + DUM(I, 1,L) 1753. 620 DUMN = DUMN + DUM(I,JM,L) 1754. UM2(IM,1,L) = UM0(IM,1,L) + DUMS 1755. UM2(1,JM,L) = UM0(1,JM,L) + DUMN 1756. C**** V component 1757. DO 630 I=1,IM*(JM-1) 1758. 630 VM2(I,1,L) = VM0(I,1,L) + DVM(I,1,L) 1759. 640 continue 1760. RETURN 1761. END 1800. 1801. SUBROUTINE OPGF (UM,VM,DT1) 1802. C**** 1803. C**** OPGF adds the pressure gradient force to the momentum 1804. C**** Input: G0M (J), GZM, S0M (kg), SZM, DT (s), MO (kg/m**2) 1805. C**** Output: UM (kg*m/s) = UM - DT*(DH*D(P)+MO*D(PHI))*DYP 1806. C**** VM (kg*m/s) = VM - DT*(DH*D(P)+MO*D(PHI))*DXV 1807. C**** 1808. INCLUDE 'C070.COM' 1809. PARAMETER (RRT12=.28867513d0) ! RRT12 = 1/SQRT(12) 1810. REAL*8 UM(IM,JM,LMO),VM(IM,JM,LMO) 1813. COMMON /WORK02/ MMI(IM,JM,LMO) 1814. COMMON /WORK03/ DUM(IM,JM,LMO),DVM(IM,JM,LMO) 1815. COMMON /WORK04/ P(IM,JM,LMO),PHI(IM,JM,LMO),DH(IM,JM,LMO) 1816. COMMON /WORK05/ DZGDP(IM,JM,LMO),VBAR(IM,JM,LMO) 1822. C**** 1822.1 C**** SURFCB BLDATA(6) Atmospheric and sea ice pressure (Pa-101325) 1822.2 C**** BLDATA(8) Ocean surface geopotential (m^2/s^2) 1822.3 C**** 1824. DT2 = DT1/2. 1825. C**** 1826. C**** Calculate the mass weighted pressure P (Pa), 1827. C**** geopotential PHI (m**2/s**2), and layer thickness DH (m) 1828. C**** 1829. DO 130 I=IM+1,IM*(JM-1)+1 1830. IF(FOCEAN(I,1).le.0.) GO TO 130 1831. C**** Calculate pressure by integrating from the top down 1832. P(I,1,1) = BLDATA(I,1,6) + MO(I,1,1)*.5*GRAV 1833. DO 110 L=2,LMM(I,1) 1834. 110 P(I,1,L) = P(I,1,L-1) + (MO(I,1,L-1)+MO(I,1,L))*.5*GRAV 1835. C**** Calculate geopotential by integrating from the bottom up 1837. PHIE = ZSOLID(I,1)*GRAV 1838. DO 120 L=LMM(I,1),1,-1 1839. PHI(I,1,L) = PHIE + DZGDP(I,1,L)*MO(I,1,L)*GRAV 1840. DH(I,1,L) = VBAR(I,1,L)*MO(I,1,L) 1841. 120 PHIE = PHIE + VBAR(I,1,L)*MO(I,1,L)*GRAV 1842. BLDATA(I,1,8) = PHIE 1843. 130 continue 1844. C**** Define polar values at all latitudes 1845. C P(I, 1,L) = P(1,IM,L) 1846. C PHI(I, 1,L) = PHI(1,IM,L) 1847. C DH(I, 1,L) = DH(1,IM,L) 1848. DO 140 L=1,LMM(1,JM) 1849. DO 140 I=2,IM 1850. P(I,JM,L) = P(1,JM,L) 1851. PHI(I,JM,L) = PHI(1,JM,L) 1852. 140 DH(I,JM,L) = DH(1,JM,L) 1853. C**** 1854. C**** Calculate smoothed East-West Pressure Gradient Force 1855. C**** 1856. I=IM 1857. DO 230 J=2,JM-1 1858. DO 230 IP1=1,IM 1859. DO 210 L=1,LMU(I,J) 1860. 210 DUM(I,J,L) = -DT2*DYP(J)* 1861. * (( DH(IP1,J,L)+ DH(I,J,L))*( P(IP1,J,L)- P(I,J,L)) + 1862. + (PHI(IP1,J,L)-PHI(I,J,L))*(MO(IP1,J,L)+MO(I,J,L))) 1863. DO 220 L=LMU(I,J)+1,LMO 1864. 220 DUM(I,J,L) = 0. 1865. 230 I=IP1 1866. CALL OPFIL (DUM) 1867. C**** 1868. C**** Calculate North-South Pressure Gradient Force 1869. C**** 1870. DO 310 J=1,JM-1 1871. DO 310 I=1,IM 1872. DO 310 L=1,LMV(I,J) 1873. 310 DVM(I,J,L) = -DT2*DXV(J)* 1874. * (( DH(I,J+1,L)+ DH(I,J,L))*( P(I,J+1,L)- P(I,J,L)) + 1875. + (PHI(I,J+1,L)-PHI(I,J,L))*(MO(I,J+1,L)+MO(I,J,L))) 1876. C**** 1877. C**** Add pressure gradient force to momentum 1878. C**** 1879. DO 610 I=IM+1,IM*(JM-1) 1880. DO 610 L=1,LMU(I,1) 1881. 610 UM(I,1,L) = UM(I,1,L) + DUM(I,1,L) 1882. DO 630 I=1,IM*(JM-1) 1883. DO 630 L=1,LMV(I,1) 1884. 630 VM(I,1,L) = VM(I,1,L) + DVM(I,1,L) 1885. RETURN 1886. C**** 1887. C**** 1888. ENTRY OPGF0 1889. C**** 1890. C**** OPGF0 calculates VOLGSP at two locations inside each grid box, 1891. C**** which will be kept constant during the dynamics of momentum, 1892. C**** and calculates two weighted averages from those values 1893. C**** 1894. DO 710 I=IM+1,IM*(JM-1)+1 1895. PE = BLDATA(I,1,6) 1896. DO 710 L=1,LMM(I,1) 1897. PBAR= PE + MO(I,1,L)*(GRAV*.5) 1898. PE = PE + MO(I,1,L)* GRAV 1899. PUP = PBAR - MO(I,1,L)*(GRAV*RRT12) 1900. PDN = PBAR + MO(I,1,L)*(GRAV*RRT12) 1901. GUP = (G0M(I,1,L)-2.*RRT12*GZM(I,1,L))/MMI(I,1,L) 1902. GDN = (G0M(I,1,L)+2.*RRT12*GZM(I,1,L))/MMI(I,1,L) 1903. SUP = (S0M(I,1,L)-2.*RRT12*SZM(I,1,L))/MMI(I,1,L) 1904. SDN = (S0M(I,1,L)+2.*RRT12*SZM(I,1,L))/MMI(I,1,L) 1905. VUP = VOLGSP(GUP,SUP,PUP) 1906. VDN = VOLGSP(GDN,SDN,PDN) 1907. DZGDP(I,1,L) = (VUP*(.5-RRT12) + VDN*(.5+RRT12))*.5 1908. 710 VBAR(I,1,L) = (VUP + VDN)*.5 1909. RETURN 1910. END 2000. 2001. SUBROUTINE OADVT (RM,RX,RY,RZ,DT,QLIMIT,OIJL) 2002. C**** 2003. C**** OADVT advects tracers using the linear upstream scheme. 2004. C**** 2005. C**** Input: MB (kg) = mass before advection 2006. C**** DT (s) = time step 2007. C**** MU (kg/s) = west to east mass flux 2008. C**** MV (kg/s) = south to north mass flux 2009. C**** MW (kg/s) = downward vertical mass flux 2010. C**** QLIMIT = whether slope limitations should be used 2011. C**** Output: RM (kg) = tracer mass 2012. C**** RX,RY,RZ (kg) = first moments of tracer mass 2013. C**** OIJL (kg) = diagnostic accumulation of tracer mass flux 2014. C**** 2014.5 IMPLICIT REAL*8 (A-H,M,O-Z) 2015. PARAMETER (IM=72,JM=46,LMO=13) 2016. REAL*8 RM(IM,JM,LMO),RX(IM,JM,LMO),RY(IM,JM,LMO),RZ(IM,JM,LMO), 2017. * OIJL(IM,JM,LMO,3) 2018. COMMON /WORK02/ MB(IM,JM,LMO), 2019. * SMU(IM,JM,LMO),SMV(IM,JM,LMO),SMW(IM,JM,LMO-1) 2020. COMMON /WORK03/ MA(IM,JM,LMO) 2021. C**** 2022. C**** Load mass after advection from mass before advection 2023. C**** 2024. DO 10 L=1,LMO 2025. DO 10 I=IM,IM*(JM-1)+1 2026. 10 MA(I,1,L) = MB(I,1,L) 2027. C**** 2028. C**** Advect the tracer using the Slopes Scheme 2029. C**** 2030. CALL OADVTX (RM,RX,RY,RZ,MA,SMU,.5*DT,QLIMIT,OIJL) 2031. CALL OADVTY (RM,RX,RY,RZ,MA,SMV, DT,QLIMIT,OIJL(1,1,1,2)) 2032. CALL OADVTZ (RM,RX,RY,RZ,MA,SMW, DT,QLIMIT,OIJL(1,1,1,3)) 2033. CALL OADVTX (RM,RX,RY,RZ,MA,SMU,.5*DT,QLIMIT,OIJL) 2034. C**** Fill in values at the poles 2035. DO 20 L=1,LMO 2036. DO 20 I=1,IM 2037. C RM(I, 1,L) = RM(IM,1,L) 2038. C RZ(I, 1,L) = RZ(IM,1,L) 2039. RM(I,JM,L) = RM(1,JM,L) 2040. 20 RZ(I,JM,L) = RZ(1,JM,L) 2041. RETURN 2042. END 2200. 2201. SUBROUTINE OADVTX (RM,RX,RY,RZ,MO,MU,DT,QLIMIT,OIJL) 2202. C**** 2203. C**** OADVTX advects tracers in the west to east direction using the 2204. C**** linear upstream scheme. If QLIMIT is true, the gradients are 2205. C**** limited to prevent the mean tracer from becoming non-negative. 2206. C**** 2207. C**** Input: DT (s) = time step 2208. C**** MU (kg/s) = west to east mass flux 2209. C**** QLIMIT = whether slope limitations should be used 2210. C**** Input and Output: RM (kg) = tracer mass 2211. C**** RX,RY,RZ (kg) = first moments of tracer mass 2212. C**** M (kg) = water mass 2213. C**** 2213.5 IMPLICIT REAL*8 (A-H,M,O-Z) 2214. PARAMETER (IM=72,JM=46,LMO=13) 2215. REAL*8 RM(IM,JM,LMO),RX(IM,JM,LMO),RY(IM,JM,LMO),RZ(IM,JM,LMO), 2217. * MO(IM,JM,LMO),MU(IM,JM,LMO), OIJL(IM,JM,LMO) 2218. LOGICAL*4 QLIMIT 2219. COMMON /FIXDCB/ FANDZ(IM,JM,8),LMM(IM,JM) 2220. COMMON /WORK04/ AM(IM),A(IM),FM(IM),FX(IM),FY(IM),FZ(IM) 2221. C**** Loop over layers and latitudes 2222. DO 320 L=1,LMO 2223. DO 320 J=2,JM-1 2224. C**** 2225. C**** Calculate FM (kg), FX (kg**2), FY (kg) and FZ (kg) 2226. C**** 2227. I=IM 2228. DO 140 IP1=1,IM 2229. AM(I) = DT*MU(I,J,L) 2230. IF(AM(I)) 110,120,130 2231. C**** Water mass flux is negative 2232. 110 A(I) = AM(I)/MO(IP1,J,L) 2233. IF(A(I).LT.-1.) WRITE (6,*) 'A<-1:',I,J,L,A(I),MO(IP1,J,L) 2234. FM(I) = A(I)*(RM(IP1,J,L)-(1.+A(I))*RX(IP1,J,L)) 2235. FX(I) = AM(I)*(A(I)*A(I)*RX(IP1,J,L)-3.*FM(I)) 2236. FY(I) = A(I)*RY(IP1,J,L) 2237. FZ(I) = A(I)*RZ(IP1,J,L) 2238. GO TO 140 2239. C**** Water mass flux is zero 2240. 120 A(I) = 0. 2241. FM(I) = 0. 2242. FX(I) = 0. 2243. FY(I) = 0. 2244. FZ(I) = 0. 2245. GO TO 140 2246. C**** Water mass flux is positive 2247. 130 A(I) = AM(I)/MO(I,J,L) 2248. IF(A(I).GT.1.) WRITE (6,*) 'A>1:',I,J,L,A(I),MO(I,J,L) 2249. FM(I) = A(I)*(RM(I,J,L)+(1.-A(I))*RX(I,J,L)) 2250. FX(I) = AM(I)*(A(I)*A(I)*RX(I,J,L)-3.*FM(I)) 2251. FY(I) = A(I)*RY(I,J,L) 2252. FZ(I) = A(I)*RZ(I,J,L) 2253. 140 I=IP1 2254. C**** 2255. C**** Modify the tracer moments so that the tracer mass in each 2256. C**** division is non-negative 2257. C**** 2258. IF(.not.QLIMIT) GO TO 300 2259. IM1=IM 2260. DO 290 I=1,IM 2261. IF(A(IM1).GE.0.) GO TO 240 2262. C**** Water is leaving through the left edge: 2 or 3 divisions 2263. IF(FM(IM1).le.0.) GO TO 210 2264. C**** Left most division is negative, RML = -FM(I-1) < 0: Case 2 or 4 2265. RX(I,J,L) = RM(I,J,L)/(1.+A(IM1)) 2266. FM(IM1) = 0. 2267. FX(IM1) = AM(IM1)*A(IM1)*A(IM1)*RX(I,J,L) 2268. IF(A(I).le.0.) GO TO 290 2269. FM(I) = A(I)*(RM(I,J,L)+(1.-A(I))*RX(I,J,L)) 2270. FX(I) = AM(I)*(A(I)*A(I)*RX(I,J,L)-3.*FM(I)) 2271. GO TO 290 2272. C**** Left most division is non-negative, RML = -FM(I-() > 0: 2273. C**** Case 1, 3 or 5 2274. 210 IF(A(I).le.0.) GO TO 230 2275. C**** Water is leaving through the right edge: 3 divisions 2276. IF(FM(I).GE.0.) GO TO 290 2277. C**** Right most division is negative, RMR = FM(I) < 0: Case 3 or 5 2278. 220 RX(I,J,L) = -RM(I,J,L)/(1.-A(I)) 2279. FM(I) = 0. 2280. FX(I) = AM(I)*A(I)*A(I)*RX(I,J,L) 2281. FM(IM1) = A(IM1)*(RM(I,J,L)-(1.+A(IM1))*RX(I,J,L)) 2282. FX(IM1) = AM(IM1)*(A(IM1)*A(IM1)*RX(I,J,L)-3.*FM(IM1)) 2283. GO TO 290 2284. C**** No water is leaving through the right edge: 2 divisions 2285. 230 IF(RM(I,J,L)+FM(IM1).GE.0.) GO TO 290 2286. C**** Right most division is negative, RMR = RM(I,J)+FM(I-1) < 0: Case 3 2287. RX(I,J,L) = RM(I,J,L)/A(IM1) 2288. FM(IM1) = -RM(I,J,L) 2289. FX(IM1) = AM(IM1)*(A(IM1)+3.)*RM(I,J,L) 2290. GO TO 290 2291. C**** No water is leaving through the left edge: 1 or 2 divisions 2292. 240 IF(A(I).le.0.) GO TO 290 2293. C**** Water is leaving through the right edge: 2 divisions 2294. IF(FM(I).GE.0.) GO TO 250 2295. C**** Right most division is negative, RMR = FM(I) < 0: Case 3 2296. RX(I,J,L) = -RM(I,J,L)/(1.-A(I)) 2297. FM(I) = 0. 2298. FX(I) = AM(I)*A(I)*A(I)*RX(I,J,L) 2299. GO TO 290 2300. C**** Right most division is non-negative, RMR = FM(I) > 0: Case 1 or 2 2301. 250 IF(RM(I,J,L)-FM(I).GE.0.) GO TO 290 2302. C**** Left most division is negative, RML = RM(I,J)-FM(I) < 0: Case 2 2303. RX(I,J,L) = RM(I,J,L)/A(I) 2304. FM(I) = RM(I,J,L) 2305. FX(I) = AM(I)*(A(I)-3.)*RM(I,J,L) 2306. C**** 2307. 290 IM1=I 2308. C**** 2309. C**** Calculate new tracer mass and first moments of tracer mass 2310. C**** 2311. 300 IM1=IM 2312. DO 310 I=1,IM 2313. IF(L.GT.LMM(I,J)) GO TO 310 2314. RM(I,J,L) = RM(I,J,L) + FM(IM1)-FM(I) 2315. RX(I,J,L) = (RX(I,J,L)*MO(I,J,L) + (FX(IM1)-FX(I)) 2316. * + 3.*((AM(IM1)+AM(I))*RM(I,J,L)-MO(I,J,L)*(FM(IM1)+FM(I)))) 2317. * / (MO(I,J,L)+AM(IM1)-AM(I)) 2318. RY(I,J,L) = RY(I,J,L) + (FY(IM1)-FY(I)) 2319. RZ(I,J,L) = RZ(I,J,L) + (FZ(IM1)-FZ(I)) 2320. MO(I,J,L) = MO(I,J,L) + AM(IM1)-AM(I) 2321. IF(MO(I,J,L).le.0.) GO TO 800 2322. IF(QLIMIT .and. RM(I,J,L).LT.0.) GO TO 810 2323. OIJL(I,J,L) = OIJL(I,J,L) + FM(I) 2324. 310 IM1=I 2325. 320 continue 2326. RETURN 2327. 800 WRITE (6,*) 'MO<0 in OADVTX:',I,J,L,MO(I,J,L) 2328. 810 WRITE (6,*) 'RM in OADVTX:',I,J,L,RM(I,J,L) 2329. WRITE (6,*) 'A=',(I,A(I),I=1,IM) 2330. STOP 2331. END 2400. 2401. SUBROUTINE OADVTY (RM,RX,RY,RZ,MO,MV,DT,QLIMIT,OIJL) 2402. C**** 2403. C**** OADVTY advects tracers in the south to north direction using the 2404. C**** linear upstream scheme. If QLIMIT is true, the gradients are 2405. C**** limited to prevent the mean tracer from becoming non-negative. 2406. C**** 2407. C**** Input: DT (s) = time step 2408. C**** MV (kg/s) = south to north mass flux 2409. C**** QLIMIT = whether slope limitations should be used 2410. C**** Input and Output: RM (kg) = tracer mass 2411. C**** RX,RY,RZ (kg) = first moments of tracer mass 2412. C**** M (kg) = water mass 2413. C**** 2413.5 IMPLICIT REAL*8 (A-H,M,O-Z) 2414. PARAMETER (IM=72,JM=46,LMO=13) 2415. REAL*8 RM(IM,JM,LMO),RX(IM,JM,LMO),RY(IM,JM,LMO),RZ(IM,JM,LMO), 2417. * MO(IM,JM,LMO),MV(IM,JM,LMO), OIJL(IM,JM,LMO) 2418. LOGICAL*4 QLIMIT 2419. COMMON /FIXDCB/ FANDZ(IM,JM,8),LMM(IM,JM) 2420. COMMON /WORK04/ BM(JM),B(JM),FM(JM),FX(JM),FY(JM),FZ(JM) 2421. C**** Loop over layers and longitudes 2422. DO 350 L=1,LMO 2423. C SBMS = 0. 2424. C SFMS = 0. 2425. C SFZS = 0. 2426. SBMN = 0. 2427. SFMN = 0. 2428. SFZN = 0. 2429. DO 330 I=1,IM 2430. C**** 2431. C**** Calculate FM (kg), FX (kg), FY (kg**2) and FZ (kg) near South Pole 2432. C**** 2433. BM(1) = DT*MV(I,1,L) 2434. C IF(BM(1)) 101,102,103 2435. C**** Water mass flux is negative 2436. C 101 B(1) = BM(1)/MO(I,2,L) 2437. C IF(B(1).LT.-1.) WRITE (6,*) 'B<-1:',I,1,L,B(1),MO(1,2,L) 2438. C FM(1) = B(1)*(RM(I,2,L)-(1.+B(1))*RY(I,2,L)) 2439. C FX(1) = B(1)*RX(I,2,L) 2440. C FY(1) = BM(1)*(B(1)*B(1)*RY(I,2,L)-3.*FM(1)) 2441. C FZ(1) = B(1)*RZ(I,2,L) 2442. C GO TO 110 2443. C**** Water mass flux is zero 2444. 102 B(1) = 0. 2445. FM(1) = 0. 2446. FX(1) = 0. 2447. FY(1) = 0. 2448. FZ(1) = 0. 2449. C GO TO 110 2450. C**** Water mass flux is positive 2451. C 103 B(1) = BM(1)/MO(IM,1,L) 2452. C IF(B(1).GT.1.) WRITE (6,*) 'B>1:',I,1,L,B(1),MO(IM,1,L) 2453. C FM(1) = B(1)*RM(IM,1,L) 2454. C FX(1) = 0. 2455. C FY(1) = -3.*BM(1)*FM(1) 2456. C FZ(1) = B(1)*RZ(IM,1,L) 2457. C**** 2458. C**** Calculate FM (kg), FX (kg), FY (kg**2) and FZ (kg) in the interior 2459. C**** 2460. 110 DO 114 J=2,JM-2 2461. BM(J) = DT*MV(I,J,L) 2462. IF(BM(J)) 111,112,113 2463. C**** Water mass flux is negative 2464. 111 B(J) = BM(J)/MO(I,J+1,L) 2465. IF(B(J).LT.-1.) WRITE (6,*) 'B<-1:',I,J,L,B(J),MO(I,J+1,L) 2466. FM(J) = B(J)*(RM(I,J+1,L)-(1.+B(J))*RY(I,J+1,L)) 2467. FX(J) = B(J)*RX(I,J+1,L) 2468. FY(J) = BM(J)*(B(J)*B(J)*RY(I,J+1,L)-3.*FM(J)) 2469. FZ(J) = B(J)*RZ(I,J+1,L) 2470. GO TO 114 2471. C**** Water mass flux is zero 2472. 112 B(J) = 0. 2473. FM(J) = 0. 2474. FX(J) = 0. 2475. FY(J) = 0. 2476. FZ(J) = 0. 2477. GO TO 114 2478. C**** Water mass flux is positive 2479. 113 B(J) = BM(J)/MO(I,J,L) 2480. IF(B(J).GT.1.) WRITE (6,*) 'B>1:',I,J,L,B(J),MO(I,J,L) 2481. FM(J) = B(J)*(RM(I,J,L)+(1.-B(J))*RY(I,J,L)) 2482. FX(J) = B(J)*RX(I,J,L) 2483. FY(J) = BM(J)*(B(J)*B(J)*RY(I,J,L)-3.*FM(J)) 2484. FZ(J) = B(J)*RZ(I,J,L) 2485. 114 continue 2486. C**** 2487. C**** Calculate FM (kg), FX (kg), FY (kg**2) and FZ (kg) near North Pole 2488. C**** 2489. BM(JM-1) = DT*MV(I,JM-1,L) 2490. IF(BM(JM-1)) 121,122,123 2491. C**** Water mass flux is negative 2492. 121 B(JM-1) = BM(JM-1)/MO(1,JM,L) 2493. IF(B(JM-1).LT.-1.) WRITE (6,*) 'B<-1:',I,JM-1,L,B(JM-1),MO(1,JM,L) 2494. FM(JM-1) = B(JM-1)*RM(1,JM,L) 2495. FX(JM-1) = 0. 2496. FY(JM-1) = -3.*BM(JM-1)*FM(JM-1) 2497. FZ(JM-1) = B(JM-1)*RZ(1,JM,L) 2498. GO TO 200 2499. C**** Water mass flux is zero 2500. 122 B(JM-1) = 0. 2501. FM(JM-1) = 0. 2502. FX(JM-1) = 0. 2503. FY(JM-1) = 0. 2504. FZ(JM-1) = 0. 2505. GO TO 200 2506. C**** Water mass flux is positive 2507. 123 B(JM-1) = BM(JM-1)/MO(I,JM-1,L) 2508. IF(B(JM-1).GT.1.) WRITE(6,*) 'B>1:',I,JM-1,L,B(JM-1),MO(1,JM-1,L) 2509. FM(JM-1) = B(JM-1)*(RM(I,JM-1,L)+(1.-B(JM-1))*RY(I,JM-1,L)) 2510. FX(JM-1) = B(JM-1)*RX(I,JM-1,L) 2511. FY(JM-1) = BM(JM-1)*(B(JM-1)*B(JM-1)*RY(I,JM-1,L)-3.*FM(JM-1)) 2512. FZ(JM-1) = B(JM-1)*RZ(I,JM-1,L) 2513. C**** 2514. C**** Modify the tracer moments so that the tracer mass in each 2515. C**** division is non-negative 2516. C**** 2517. 200 IF(.not.QLIMIT) GO TO 300 2518. DO 290 J=2,JM-1 2519. IF(B(J-1).GE.0.) GO TO 240 2520. C**** Water is leaving through the south edge: 2 or 3 divisions 2521. IF(FM(J-1).le.0.) GO TO 210 2522. C**** South most division is negative, RMS = -FM(J-1) < 0: Case 2 or 4 2523. RY(I,J,L) = RM(I,J,L)/(1.+B(J-1)) 2524. FM(J-1) = 0. 2525. FY(J-1) = BM(J-1)*B(J-1)*B(J-1)*RY(I,J,L) 2526. IF(B(J).le.0.) GO TO 290 2527. FM(J) = B(J)*(RM(I,J,L)+(1.-B(J))*RY(I,J,L)) 2528. FY(J) = BM(J)*(B(J)*B(J)*RY(I,J,L)-3.*FM(J)) 2529. GO TO 290 2530. C**** South most division is non-negative, RMS = -FM(J-1) > 0: 2531. C**** Case 1, 3 or 5 2532. 210 IF(B(J).le.0.) GO TO 230 2533. C**** Water is leaving through the north edge: 3 divisions 2534. IF(FM(J).GE.0.) GO TO 290 2535. C**** North most division is negative, RMN = FM(J) < 0: Case 3 or 5 2536. 220 RY(I,J,L) = -RM(I,J,L)/(1.-B(J)) 2537. FM(J) = 0. 2538. FY(J) = BM(J)*B(J)*B(J)*RY(I,J,L) 2539. FM(J-1) = B(J-1)*(RM(I,J,L)-(1.+B(J-1))*RY(I,J,L)) 2540. FY(J-1) = BM(J-1)*(B(J-1)*B(J-1)*RY(I,J,L)-3.*FM(J-1)) 2541. GO TO 290 2542. C**** No water is leaving through the north edge: 2 divisions 2543. 230 IF(RM(I,J,L)+FM(J-1).GE.0.) GO TO 290 2544. C**** North most division is negative, RMN = RM(I,J)+FM(J-1) < 0: Case 3 2545. RY(I,J,L) = RM(I,J,L)/B(J-1) 2546. FM(J-1) = -RM(I,J,L) 2547. FY(J-1) = BM(J-1)*(B(J-1)+3.)*RM(I,J,L) 2548. GO TO 290 2549. C**** No water is leaving through the south edge: 1 or 2 divisions 2550. 240 IF(B(J).le.0.) GO TO 290 2551. C**** Water is leaving through the north edge: 2 divisions 2552. IF(FM(J).GE.0.) GO TO 250 2553. C**** North most division is negative, RMN = FM(J) < 0: Case 3 2554. RY(I,J,L) = -RM(I,J,L)/(1.-B(J)) 2555. FM(J) = 0. 2556. FY(J) = BM(J)*B(J)*B(J)*RY(I,J,L) 2557. GO TO 290 2558. C**** North most division is non-negative, RMN = FM(J) > 0: Case 1 or 2 2559. 250 IF(RM(I,J,L)-FM(J).GE.0.) GO TO 290 2560. C**** South most division is negative, RMS = RM(I,J)-FM(J) < 0: Case 2 2561. RY(I,J,L) = RM(I,J,L)/B(J) 2562. FM(J) = RM(I,J,L) 2563. FY(J) = BM(J)*(B(J)-3.)*RM(I,J,L) 2564. C**** 2565. 290 continue 2566. C 300 SBMS = SBMS + BM(1) 2567. C SFMS = SFMS + FM(1) 2568. C SFZS = SFZS + FZ(1) 2569. 300 SBMN = SBMN + BM(JM-1) 2570. SFMN = SFMN + FM(JM-1) 2571. SFZN = SFZN + FZ(JM-1) 2572. C**** 2573. C**** Calculate new tracer mass and first moments of tracer mass 2574. C**** 2575. DO 310 J=2,JM-1 2576. IF(L.GT.LMM(I,J)) GO TO 310 2577. RM(I,J,L) = RM(I,J,L) + FM(J-1)-FM(J) 2578. RX(I,J,L) = RX(I,J,L) + (FX(J-1)-FX(J)) 2579. RY(I,J,L) = (RY(I,J,L)*MO(I,J,L) + (FY(J-1)-FY(J)) 2580. * + 3.*((BM(J-1)+BM(J))*RM(I,J,L)-MO(I,J,L)*(FM(J-1)+FM(J)))) 2581. * / (MO(I,J,L)+BM(J-1)-BM(J)) 2582. RZ(I,J,L) = RZ(I,J,L) + (FZ(J-1)-FZ(J)) 2583. MO(I,J,L) = MO(I,J,L) + BM(J-1)-BM(J) 2584. IF(MO(I,J,L).le.0.) GO TO 800 2585. IF(QLIMIT.and.RM(I,J,L).LT.0.) GO TO 810 2586. 310 continue 2587. DO 320 J=1,JM-1 2588. 320 OIJL(I,J,L) = OIJL(I,J,L) + FM(J) 2589. 330 continue 2590. C IF(L.GT.LMM(IM,1)) GO TO 340 2591. C RM(IM,1,L) = RM(IM,1,L) - SFMS/IM 2592. C RZ(IM,1,L) = RZ(IM,1,L) - SFZS/IM 2593. C MO(IM,1,L) = MO(IM,1,L) - SBMS/IM 2594. C IF(MO(IM,1,L).le.0.) GO TO 800 2595. C IF(QLIMIT.and.RM(IM,1,L).LT.0.) GO TO 810 2596. 340 IF(L.GT.LMM(1,JM)) GO TO 350 2597. RM(1,JM,L) = RM(1,JM,L) + SFMN/IM 2598. RZ(1,JM,L) = RZ(1,JM,L) + SFZN/IM 2599. MO(1,JM,L) = MO(1,JM,L) + SBMN/IM 2600. IF(MO(1,JM,L).le.0.) GO TO 800 2601. IF(QLIMIT.and.RM(1,JM,L).LT.0.) GO TO 810 2602. 350 continue 2603. RETURN 2604. 800 WRITE (6,*) 'MO<0 in OADVTY:',I,J,L,MO(I,J,L) 2605. 810 WRITE (6,*) 'RM in ADVSSY:',I,J,L,RM(I,J,L) 2606. WRITE (6,*) 'B=',(J,B(J),J=1,JM-1) 2607. STOP 2608. END 2700. 2701. SUBROUTINE OADVTZ (RM,RX,RY,RZ,MO,MW,DT,QLIMIT,OIJL) 2702. C**** 2703. C**** OADVTZ advects tracers in the vertical direction using the 2704. C**** linear upstream scheme. If QLIMIT is true, the gradients are 2705. C**** limited to prevent the mean tracer from becoming non-negative. 2706. C**** 2707. C**** Input: DT (s) = time step 2708. C**** MW (kg/s) = downward vertical mass flux 2709. C**** QLIMIT = whether slope limitations should be used 2710. C**** Input and Output: RM (kg) = tracer mass 2711. C**** RX,RY,RZ (kg) = first moments of tracer mass 2712. C**** M (kg) = water mass 2713. C**** 2713.5 IMPLICIT REAL*8 (A-H,M,O-Z) 2714. PARAMETER (IM=72,JM=46,LMO=13) 2715. REAL*8 RM(IM,JM,LMO),RX(IM,JM,LMO),RY(IM,JM,LMO),RZ(IM,JM,LMO), 2717. * MO(IM,JM,LMO),MW(IM,JM,LMO-1), OIJL(IM,JM,LMO) 2718. LOGICAL*4 QLIMIT 2719. COMMON /FIXDCB/ FANDZ(IM,JM,8),LMM(IM,JM) 2720. COMMON /WORK04/ CM(0:LMO),C(0:LMO), 2721. * FM(0:LMO),FX(0:LMO),FY(0:LMO),FZ(0:LMO) 2722. C**** 2723. CM(0) = 0. 2724. C(0) = 0. 2725. FM(0) = 0. 2726. FX(0) = 0. 2727. FY(0) = 0. 2728. FZ(0) = 0. 2729. C**** Loop over latitudes and longitudes 2730. DO 330 I=IM,IM*(JM-1)+1 2731. LMIJ=LMM(I,1) 2732. IF(LMIJ.le.1) GO TO 330 2733. CM(LMIJ) = 0. 2734. C(LMIJ) = 0. 2735. FM(LMIJ) = 0. 2736. FX(LMIJ) = 0. 2737. FY(LMIJ) = 0. 2738. FZ(LMIJ) = 0. 2739. C**** 2740. C**** Calculate FM (kg), FX (kg), FY (kg) and FZ (kg**2) 2741. C**** 2742. DO 120 L=1,LMIJ-1 2743. CM(L) = DT*MW(I,1,L) 2744. IF(CM(L).LT.0.) GO TO 110 2745. C**** Water mass flux is positive 2746. C(L) = CM(L)/MO(I,1,L) 2747. IF(C(L).GT.1.) WRITE (6,*) 'C>1:',I,L,C(L),MO(I,1,L) 2748. FM(L) = C(L)*(RM(I,1,L)+(1.-C(L))*RZ(I,1,L)) 2749. FX(L) = C(L)*RX(I,1,L) 2750. FY(L) = C(L)*RY(I,1,L) 2751. FZ(L) = CM(L)*(C(L)*C(L)*RZ(I,1,L)-3.*FM(L)) 2752. GO TO 120 2753. C**** Water mass flux is negative 2754. 110 C(L) = CM(L)/MO(I,1,L+1) 2755. IF(C(L).LT.-1.) WRITE (6,*) 'C<-1:',I,L,C(L),MO(I,1,L+1) 2756. FM(L) = C(L)*(RM(I,1,L+1)-(1.+C(L))*RZ(I,1,L+1)) 2757. FX(L) = C(L)*RX(I,1,L+1) 2758. FY(L) = C(L)*RY(I,1,L+1) 2759. FZ(L) = CM(L)*(C(L)*C(L)*RZ(I,1,L+1)-3.*FM(L)) 2760. 120 continue 2761. C**** 2762. C**** Modify the tracer moments so that the tracer mass in each 2763. C**** division is non-negative 2764. C**** 2765. IF(.not.QLIMIT) GO TO 300 2766. DO 290 L=1,LMIJ 2767. IF(C(L-1).GE.0.) GO TO 240 2768. C**** Water is leaving through the bottom edge: 2 or 3 divisions 2769. IF(FM(L-1).le.0.) GO TO 210 2770. C**** Bottom most division is negative, RMB = -FM(L-1) < 0: Case 2 or 4 2771. RZ(I,1,L) = RM(I,1,L)/(1.+C(L-1)) 2772. FM(L-1) = 0. 2773. FZ(L-1) = CM(L-1)*C(L-1)*C(L-1)*RZ(I,1,L) 2774. IF(C(L).le.0.) GO TO 290 2775. FM(L) = C(L)*(RM(I,1,L)+(1.-C(L))*RZ(I,1,L)) 2776. FZ(L) = CM(L)*(C(L)*C(L)*RZ(I,1,L)-3.*FM(L)) 2777. GO TO 290 2778. C**** Bottom most division is non-negative, RMB = -FM(L-1) > 0: 2779. C**** Case 1, 3 or 5 2780. 210 IF(C(L).le.0.) GO TO 230 2781. C**** Water is leaving through the top edge: 3 divisions 2782. IF(FM(L).GE.0.) GO TO 290 2783. C**** Top most division is negative, RMT = FM(L) < 0: Case 3 or 5 2784. 220 RZ(I,1,L) = -RM(I,1,L)/(1.-C(L)) 2785. FM(L) = 0. 2786. FZ(L) = CM(L)*C(L)*C(L)*RZ(I,1,L) 2787. FM(L-1) = C(L-1)*(RM(I,1,L)-(1.+C(L-1))*RZ(I,1,L)) 2788. FZ(L-1) = CM(L-1)*(C(L-1)*C(L-1)*RZ(I,1,L)-3.*FM(L-1)) 2789. GO TO 290 2790. C**** No water is leaving through the top edge: 2 divisions 2791. 230 IF(RM(I,1,L)+FM(L-1).GE.0.) GO TO 290 2792. C**** Top most division is negative, RMT = RM(I,1,L)+FM(L-1) < 0: Case 3 2793. RZ(I,1,L) = RM(I,1,L)/C(L-1) 2794. FM(L-1) = -RM(I,1,L) 2795. FZ(L-1) = CM(L-1)*(C(L-1)+3.)*RM(I,1,L) 2796. GO TO 290 2797. C**** No water is leaving through the bottom edge: 1 or 2 divisions 2798. 240 IF(C(L).le.0.) GO TO 290 2799. C**** Water is leaving through the top edge: 2 divisions 2800. IF(FM(L).GE.0.) GO TO 250 2801. C**** Top most division is negative, RMT = FM(L) < 0: Case 3 2802. RZ(I,1,L) = -RM(I,1,L)/(1.-C(L)) 2803. FM(L) = 0. 2804. FZ(L) = CM(L)*C(L)*C(L)*RZ(I,1,L) 2805. GO TO 290 2806. C**** Top most division is non-negative, RMT = FM(L) > 0: Case 1 or 2 2807. 250 IF(RM(I,1,L)-FM(L).GE.0.) GO TO 290 2808. C**** Bottom most division is negative, RMB = RM(I,1,L)-FM(L) < 0: Cas 2 2809. RZ(I,1,L) = RM(I,1,L)/C(L) 2810. FM(L) = RM(I,1,L) 2811. FZ(L) = CM(L)*(C(L)-3.)*RM(I,1,L) 2812. C**** 2813. 290 continue 2814. C**** 2815. C**** Calculate new tracer mass and first moments of tracer mass 2816. C**** 2817. 300 DO 310 L=1,LMIJ 2818. RM(I,1,L) = RM(I,1,L) + FM(L-1)-FM(L) 2819. RX(I,1,L) = RX(I,1,L) + (FX(L-1)-FX(L)) 2820. RY(I,1,L) = RY(I,1,L) + (FY(L-1)-FY(L)) 2821. RZ(I,1,L) = (RZ(I,1,L)*MO(I,1,L) + (FZ(L-1)-FZ(L)) 2822. * + 3.*((CM(L-1)+CM(L))*RM(I,1,L)-MO(I,1,L)*(FM(L-1)+FM(L)))) 2823. * / (MO(I,1,L)+CM(L-1)-CM(L)) 2824. MO(I,1,L) = MO(I,1,L) + CM(L-1)-CM(L) 2825. IF(MO(I,1,L).le.0.) GO TO 800 2826. IF(QLIMIT.and.RM(I,1,L).LT.0.) GO TO 810 2827. 310 continue 2828. DO 320 L=1,LMIJ-1 2829. 320 OIJL(I,1,L) = OIJL(I,1,L) + FM(L) 2830. 330 continue 2831. RETURN 2832. 800 WRITE (6,*) 'MO<0 in OADVTZ:',I,L,MO(I,1,L) 2833. 810 WRITE (6,*) 'RM in OADVTZ:',I,L,RM(I,1,L) 2834. WRITE (6,*) 'C=',(L,C(L),L=0,LMIJ) 2835. STOP 2836. END 3000. 3001. SUBROUTINE STPGF 3002. C**** 3003. C**** STPGF accelerates the water through the strait due to the 3004. C**** pressure gradient force between the GCM grid boxes at the 3005. C**** opposite ends of the strait. 3006. C**** 3007. INCLUDE 'C070.COM' 3008. PARAMETER (RRT12=.28867513d0) ! RRT12 = 1/SQRT(12) 3009. COMMON /WORK01/ MEND(LMO,2),PEND(LMO,2),PHI(LMO,2),DH(LMO,2) 3010. C**** 3012. DO 50 N=1,NMST 3013. C**** 3014. C**** Accelerate the channel speed due to the pressure gradient 3015. C**** force, reduce the channel speed due to friction 3016. C**** 3017. DO 30 K=1,2 3018. I=IST(N,K) 3019. J=JST(N,K) 3020. C**** Calculate pressure by integrating from the top down 3021. PE = BLDATA(I,J,6) 3022. DO 10 L=1,LMM(I,J) 3023. MEND(L,K) = MO(I,J,L) 3024. PEND(L,K) = PE + GRAV*MO(I,J,L)*.5 3025. 10 PE = PE + GRAV*MO(I,J,L) 3026. C**** Calculate geopotential by integrating from the bottom up, 3027. C**** also calculate the specific volume 3028. PHIE = ZSOLID(I,J)*GRAV 3029. DO 20 L=LMM(I,J),1,-1 3030. GUP = (G0M(I,J,L)-2.*RRT12*GZM(I,J,L))/(MO(I,J,L)*DXYP(J)) 3031. GDN = (G0M(I,J,L)+2.*RRT12*GZM(I,J,L))/(MO(I,J,L)*DXYP(J)) 3032. SUP = (S0M(I,J,L)-2.*RRT12*SZM(I,J,L))/(MO(I,J,L)*DXYP(J)) 3033. SDN = (S0M(I,J,L)+2.*RRT12*SZM(I,J,L))/(MO(I,J,L)*DXYP(J)) 3034. DP = GRAV*MEND(L,K) 3035. PUP = PEND(L,K) - RRT12*DP 3036. PDN = PEND(L,K) + RRT12*DP 3037. VUP = VOLGSP(GUP,SUP,PUP) 3038. VDN = VOLGSP(GDN,SDN,PDN) 3039. PHI(L,K) = PHIE + (VUP*(.5-RRT12)+VDN*(.5+RRT12))*.5*DP 3040. DH(L,K) = MEND(L,K)*(VUP+VDN)*.5 3041. C**** Calculate PHI at top of current layer (or bottom of next layer) 3042. 20 PHIE = PHIE + (VDN+VUP)*.5*DP 3043. 30 continue 3044. C**** Subtract Pressure Gradient Force from the mass flux 3045. DO 40 L=1,LMST(N) 3046. 40 MUST(L,N) = MUST(L,N) - .5*DTS*WIST(N)* 3047. * (( DH(L,2)+ DH(L,1))*(PEND(L,2)-PEND(L,1)) + 3048. + (PHI(L,2)-PHI(L,1))*(MEND(L,2)+MEND(L,1)))/DISTPG(N) 3049. 50 continue 3050. RETURN 3051. END 3500. 3501. SUBROUTINE STADV 3502. C**** 3503. C**** STADV advects water mass and tracers through the straits 3504. C**** 3505. INCLUDE 'C070.COM' 3506. C**** 3507. C**** STRACB MMST = Sea water mass in strait (kg) 3508. C**** MUST = Mass flux through strait (kg/s) 3509. C**** GOMST = mean potential enthalpy in strait (J) 3510. C**** GXMST = horizontal gradient of potential enthalpy (J) 3511. C**** GZMST = vertical gradient of potential enthalpy (J) 3512. C**** S0MST = mean salt in strait (kg) 3513. C**** SXMST = horizontal gradient of salt (kg) 3514. C**** SZMST = vertical gradient of salt (kg) 3515. C**** 3516. DO 310 N=1,NMST 3517. I1 = IST(N,1) 3518. J1 = JST(N,1) 3519. X1 = XST(N,1) 3520. Y1 = YST(N,1) 3521. I2 = IST(N,2) 3522. J2 = JST(N,2) 3523. X2 = XST(N,2) 3524. Y2 = YST(N,2) 3525. DO 310 L=1,LMST(N) 3526. AM = DTS*MUST(L,N) 3527. MM1 = MO(I1,J1,L)*DXYP(J1) 3528. MM2 = MO(I2,J2,L)*DXYP(J2) 3529. IF(AM.lt.0.) GO TO 200 3530. C**** 3531. C**** AM > 0: Water flux is moving from grid box 1 to grid box 2 3532. C**** 3533. A1 = AM/MM1 3534. FGM1 = A1*(G0M(I1,J1,L) + GXM(I1,J1,L)*X1 + GYM(I1,J1,L)*Y1) 3535. FSM1 = A1*(S0M(I1,J1,L) + SXM(I1,J1,L)*X1 + SYM(I1,J1,L)*Y1) 3536. FGZ1 = A1* GZM(I1,J1,L) 3537. FSZ1 = A1* SZM(I1,J1,L) 3538. A2 = AM/MMST(L,N) 3539. FGM2 = A2*(G0MST(L,N) + (1.-A2)*GXMST(L,N)) 3540. FSM2 = A2*(S0MST(L,N) + (1.-A2)*SXMST(L,N)) 3541. FGZ2 = A2* GZMST(L,N) 3542. FSZ2 = A2* SZMST(L,N) 3543. C**** Calculate first moments of tracers 3544. GXM(I1,J1,L) = GXM(I1,J1,L)*(1. - A1*(.5+1.5*X1*X1)) 3545. SXM(I1,J1,L) = SXM(I1,J1,L)*(1. - A1*(.5+1.5*X1*X1)) 3546. GYM(I1,J1,L) = GYM(I1,J1,L)*(1. - A1*(.5+1.5*Y1*Y1)) 3547. SYM(I1,J1,L) = SYM(I1,J1,L)*(1. - A1*(.5+1.5*Y1*Y1)) 3548. GXM(I2,J2,L) = GXM(I2,J2,L) + 3.*X2*FGM2 + 3549. + (AM/MM2)*(GXM(I2,J2,L)*(.5-1.5*X2*X2)-3.*X2*G0M(I2,J2,L)) 3550. SXM(I2,J2,L) = SXM(I2,J2,L) + 3.*X2*FSM2 + 3551. + (AM/MM2)*(SXM(I2,J2,L)*(.5-1.5*X2*X2)-3.*X2*S0M(I2,J2,L)) 3552. GYM(I2,J2,L) = GYM(I2,J2,L) + 3.*Y2*FGM2 + 3553. + (AM/MM2)*(GYM(I2,J2,L)*(.5-1.5*Y2*Y2)-3.*Y2*G0M(I2,J2,L)) 3554. SYM(I2,J2,L) = SYM(I2,J2,L) + 3.*Y2*FSM2 + 3555. + (AM/MM2)*(SYM(I2,J2,L)*(.5-1.5*Y2*Y2)-3.*Y2*S0M(I2,J2,L)) 3556. G0MST(L,N) = G0MST(L,N) + FGM1 - FGM2 3557. S0MST(L,N) = S0MST(L,N) + FSM1 - FSM2 3558. GXMST(L,N) = GXMST(L,N)*(1.-A2)**3 - 3.*FGM1 + 3.*A2*G0MST(L,N) 3559. SXMST(L,N) = SXMST(L,N)*(1.-A2)**3 - 3.*FSM1 + 3.*A2*S0MST(L,N) 3560. GO TO 300 3561. C**** 3562. C**** AM < 0: Water flux is moving from grid box 2 to grid box 1 3563. C**** 3564. 200 A1 = AM/MMST(L,N) 3565. FGM1 = A1*(G0MST(L,N) - (1.+A1)*GXMST(L,N)) 3566. FSM1 = A1*(S0MST(L,N) - (1.+A1)*SXMST(L,N)) 3567. FGZ1 = A1* GZMST(L,N) 3568. FSZ1 = A1* SZMST(L,N) 3569. A2 = AM/MM2 3570. FGM2 = A2*(G0M(I2,J2,L) + GXM(I2,J2,L)*X2 + GYM(I2,J2,L)*Y2) 3571. FSM2 = A2*(S0M(I2,J2,L) + SXM(I2,J2,L)*X2 + SYM(I2,J2,L)*Y2) 3572. FGZ2 = A2* GZM(I2,J2,L) 3573. FSZ2 = A2* SZM(I2,J2,L) 3574. C**** Calculate first moments of tracers 3575. GXM(I1,J1,L) = GXM(I1,J1,L) - 3.*X1*FGM1 - 3576. - (AM/MM1)*(GXM(I1,J1,L)*(.5-1.5*X1*X1)-3.*X1*G0M(I1,J1,L)) 3577. SXM(I1,J1,L) = SXM(I1,J1,L) - 3.*X1*FSM1 - 3578. - (AM/MM1)*(SXM(I1,J1,L)*(.5-1.5*X1*X1)-3.*X1*S0M(I1,J1,L)) 3579. GYM(I1,J1,L) = GYM(I1,J1,L) - 3.*Y1*FGM1 - 3580. - (AM/MM1)*(GYM(I1,J1,L)*(.5-1.5*Y1*Y1)-3.*Y1*G0M(I1,J1,L)) 3581. SYM(I1,J1,L) = SYM(I1,J1,L) - 3.*Y1*FSM1 - 3582. - (AM/MM1)*(SYM(I1,J1,L)*(.5-1.5*Y1*Y1)-3.*Y1*S0M(I1,J1,L)) 3583. GXM(I2,J2,L) = GXM(I2,J2,L)*(1. + A2*(.5+1.5*X2*X2)) 3584. SXM(I2,J2,L) = SXM(I2,J2,L)*(1. + A2*(.5+1.5*X2*X2)) 3585. GYM(I2,J2,L) = GYM(I2,J2,L)*(1. + A2*(.5+1.5*Y2*Y2)) 3586. SYM(I2,J2,L) = SYM(I2,J2,L)*(1. + A2*(.5+1.5*Y2*Y2)) 3587. G0MST(L,N) = G0MST(L,N) + FGM1 - FGM2 3588. S0MST(L,N) = S0MST(L,N) + FSM1 - FSM2 3589. GXMST(L,N) = GXMST(L,N)*(1.+A1)**3 - 3.*FGM2 + 3.*A1*G0MST(L,N) 3590. SXMST(L,N) = SXMST(L,N)*(1.+A1)**3 - 3.*FSM2 + 3.*A1*S0MST(L,N) 3591. C**** 3592. C**** Calculate new tracer mass, vertical moment of tracer mass, 3593. C**** and horizontal moment of tracer mass for the straits 3594. C**** 3595. 300 MO(I1,J1,L) = MO(I1,J1,L) - AM*BYDXYP(J1) 3596. MO(I2,J2,L) = MO(I2,J2,L) + AM*BYDXYP(J2) 3597. G0M(I1,J1,L) = G0M(I1,J1,L) - FGM1 3598. S0M(I1,J1,L) = S0M(I1,J1,L) - FSM1 3599. GZM(I1,J1,L) = GZM(I1,J1,L) - FGZ1 3600. SZM(I1,J1,L) = SZM(I1,J1,L) - FSZ1 3601. G0M(I2,J2,L) = G0M(I2,J2,L) + FGM2 3602. S0M(I2,J2,L) = S0M(I2,J2,L) + FSM2 3603. GZM(I2,J2,L) = GZM(I2,J2,L) + FGZ2 3604. SZM(I2,J2,L) = SZM(I2,J2,L) + FSZ2 3605. GZMST(L,N) = GZMST(L,N) + (FGZ1-FGZ2) 3606. SZMST(L,N) = SZMST(L,N) + (FSZ1-FSZ2) 3607. OLNST(L,N,1) = OLNST(L,N,1) + AM 3608. OLNST(L,N,2) = OLNST(L,N,2) + (FGM1+FGM2) 3609. OLNST(L,N,3) = OLNST(L,N,3) + (FSM1+FSM2) 3610. 310 continue 3611. RETURN 3612. END 5000. 5001. SUBROUTINE OBDRAG 5002. C**** 5003. C**** OBDRAG exerts a drag on the Ocean Model's bottom layer 5004. C**** 5005. C**** TAU (kg/m*s^2) = vector bottom drag stress = BDRAGC*RHOO*|WO|*WO 5006. C**** WO = WO - DTS*TAU/MO = 5007. C**** = WO - DTS*BDRAGC*RHOO*|WO|*WO/MO = 5008. C**** = WO * (1 - DTS*BDRAGC*|WO|/dZ) ~ 5009. C**** ~ WO / (1 + DTS*BDRAGC*|WO|/dZ) 5010. C**** 5011. INCLUDE 'C070.COM' 5012. REAL*8 DTSxBDRAGCzDZ(LMO), DTSxBDRAGCxDISTzMMSTxDZ(LMO,NMST), 5013. * DTSxSDRAGCxDISTzMMSTxWIST(LMO,NMST) 5014. COMMON /WORK01/ UT(IM,JM,LMO),VT(IM,JM,LMO) 5015. PARAMETER (IQ1=1+IM/4, IQ2=IQ1+IM/4, IQ3=IQ2+IM/4, 5016. * BDRAGC=.001d0, SDRAGC=.1d0) 5017. DATA IFIRST /1/ 5018. C**** 5019. IF(IFIRST.eq.0) GO TO 100 5020. IFIRST = 0 5021. DO 10 L=1,LMO 5022. 10 DTSxBDRAGCzDZ(L) = DTS*BDRAGC / (ZE(L)-ZE(L-1)) 5023. DO 20 N=1,NMST 5024. DO 20 L=1,LMST(N) 5025. DTSxBDRAGCxDISTzMMSTxDZ(L,N) = 5026. = DTS*BDRAGC*DIST(N) / (MMST(L,N)*(ZE(L)-ZE(L-1))) 5027. 20 DTSxSDRAGCxDISTzMMSTxWIST(L,N) = 5028. = DTS*SDRAGC*DIST(N) / (MMST(L,N)*WIST(N)) 5029. C**** 5030. C**** Save UO,VO into UT,VT which will be unchanged 5031. C**** 5032. 100 DO 110 I=1,IM*JM*LMO*2 5033. 110 UT(I,1,1) = UO(I,1,1) 5034. C**** 5035. C**** Reduce West-East ocean current 5036. C**** 5037. C**** Bottom drag in the interior 5038. I=IM 5039. DO 220 J=2,JM-1 5040. DO 210 IP1=1,IM 5041. IF(LMU(I,J).le.0) GO TO 210 5042. L=LMU(I,J) 5043. WSQ = UT(I,J,L)*UT(I,J,L) + 5044. * .25*(VT(I,J ,L)*VT(I,J ,L) + VT(IP1,J ,L)*VT(IP1,J ,L) 5045. * + VT(I,J-1,L)*VT(I,J-1,L) + VT(IP1,J-1,L)*VT(IP1,J-1,L)) 5046. UO(I,J,L) = UO(I,J,L) / (1. + DTSxBDRAGCzDZ(L)*SQRT(WSQ)) 5047. 210 I=IP1 5048. 220 continue 5049. C**** Bottom drag at the poles 5050. J=JM 5051. JV=JM-1 5052. C IF(LMU(1,J).le.0) GO TO 5053. L=LMU(1,J) 5054. WSQ = UT(1,J,L)*UT(1,J,L) + 5055. * .25*(VT( 1,JV,L)*VT( 1,JV,L) + VT(IQ1,JV,L)*VT(IQ1,JV,L) 5056. * + VT(IQ2,JV,L)*VT(IQ2,JV,L) + VT(IQ3,JV,L)*VT(IQ3,JV,L)) 5057. UO(1,J,L) = UO(1,J,L) / (1. + DTSxBDRAGCzDZ(L)*SQRT(WSQ)) 5058. DO 230 I=2,IM 5059. 230 UO(I,J,L) = UO(1,J,L) 5060. C**** 5061. C**** Reduce South-North ocean current 5062. C**** 5063. IM1=IM 5064. DO 320 J=1,JM-1 5065. DO 310 I=1,IM 5066. IF(LMV(I,J).le.0) GO TO 310 5067. L=LMV(I,J) 5068. WSQ = VT(I,J,L)*VT(I,J,L) + 5069. * .25*(UT(IM1,J+1,L)*UT(IM1,J+1,L) + UT(I,J+1,L)*UT(I,J+1,L) 5070. * + UT(IM1,J ,L)*UT(IM1,J ,L) + UT(I,J ,L)*UT(I,J ,L)) 5071. VO(I,J,L) = VO(I,J,L) / (1. + DTSxBDRAGCzDZ(L)*SQRT(WSQ)) 5072. 310 IM1=I 5073. 320 continue 5074. RETURN 5075. C**** 5076. C**** 5077. ENTRY STDRAG 5078. C**** 5079. C**** STDRAG exerts a bottom drag on the lowest layer and a side drag 5080. C**** on each layer of each strait. 5081. C**** UST (m/s) = strait velocity = MUST*DIST/MMST 5082. C**** 5083. C**** STRACB MMST Mass of sea water in strait (kg) 5084. C**** MUST Mass flow through strait (kg/s) 5085. C**** WIST Width of strait (m) 5086. C**** DIST Length of strait (m) 5087. C**** 5088. DO 410 N=1,NMST 5089. C**** 5090. C**** TAU (kg/m*s^2) = bottom drag stress = BDRAGC*RHOO*|UST|*UST 5091. C**** MUST = MUST - DTS*TAU*WIST = 5092. C**** = MUST - DTS*BDRAGC*RHOO*|UST|*UST*WIST = 5093. C**** = MUST - DTS*BDRAGC*RHOO*|MUST|*MUST*WIST*DIST^2/MMST^2 = 5094. C**** = MUST * (1 - DTS*BDRAGC*|MUST|*DIST/MMST*dZ) ~ 5095. C**** ~ MUST / (1 + DTS*BDRAGC*|MUST|*DIST/MMST*dZ) 5096. C**** 5097. L=LMST(N) 5098. MUST(L,N) = MUST(L,N) / 5099. * (1. + DTSxBDRAGCxDISTzMMSTxDZ(L,N)*ABS(MUST(L,N))) 5100. C**** 5101. C**** TAU (kg/m*s^2) = side drag stress = SDRAGC*RHOO*|UST|*UST 5102. C**** MUST = MUST - DTS*TAU*dZ = 5103. C**** = MUST - DTS*SDRAGC*RHOO*|UST|*UST*dZ = 5104. C**** = MUST - DTS*SDRAGC*RHOO*|MUST|*MUST*dZ*DIST^2/MMST^2 = 5105. C**** = MUST * (1 - DTS*SDRAGC*|MUST|*DIST/MMST*WIST) ~ 5106. C**** ~ MUST / (1 + DTS*SDRAGC*|MUST|*DIST/MMST*WIST) 5107. C**** 5108. DO 410 L=1,LMST(N) 5109. 410 MUST(L,N) = MUST(L,N) / 5110. * (1. + DTSxSDRAGCxDISTzMMSTxWIST(L,N)*ABS(MUST(L,N))) 5111. RETURN 5112. END 5500. 5501. SUBROUTINE OCOAST 5502. C**** 5503. C**** OCOAST reduces the horizontal gradients of tracers that are 5504. C**** perpendicular to coastline ocean grid boxes. 5505. C**** OCOAST also applies a side drag to ocean currents that are 5506. C**** parallel to the coast. The drag is multiplied by 0, .5, 1, 1.5 5507. C**** or 2 depending on how many of the four grid boxes adjacent to the 5508. C**** the ocean current boxes are land. Side drag is currently unused. 5509. C**** 5510. INCLUDE 'C070.COM' 5511. C PARAMETER (SDRAGC=1d0) 5512. C SAVE :: REDUCE, IFIRST 5513. COMMON /OCOSCB/ LMINX(IM,JM),LMINY(IM,JM) 5514. C COMMON /OCOSCB/ SDRAGU(IM,JM,LMO),SDRAGV(IM,JM,LMO), 5515. C * LMINX(IM,JM),LMINY(IM,JM),LMINU(IM,JM),LMINV(IM,JM) 5516. DATA IFIRST /1/ 5517. C**** 5518. IF(IFIRST.eq.0) GO TO 300 5519. IFIRST = 0 5520. C**** 5521. C**** Calculate first ocean grid box in vertical that is adjacent to 5522. C**** a coast line 5523. C**** 5524. REDUCE = 1. - DTS/(SDAY*20.) 5525. IM1=IM-1 5526. I=IM 5527. DO 110 J=2,JM-1 5528. DO 110 IP1=1,IM 5529. LMINX(I,J) = MIN(LMM(IM1,J),LMM(IP1,J)) + 1 5530. IM1=I 5531. 110 I=IP1 5532. DO 120 J=2,JM-1 5533. DO 120 I=1,IM 5534. 120 LMINY(I,J) = MIN(LMM(I,J-1),LMM(I,J+1)) + 1 5535. C**** 5536. C**** Calculate side drag factors that will be fixed in time 5537. C**** 5538. C**** TAUU = east-west side drag stress = SDRAGC*RHOO*|UO|*UO 5539. C**** UO = UO - DTS*TAUU/RHOO*DYP = 5540. C**** = UO - DTS*SDRAGC*|UO|*UO/DYP = 5541. C**** = UO * (1 - DTS*SDRAGC*|UO|/DYP) ~ 5542. C**** ~ UO / (1 + DTS*SDRAGC*|UO|/DYP) 5543. C I=IM 5544. C DO 220 J=2,JM-1 5545. C DO 220 IP1=1,IM 5546. C LMINU(I,J) = LMU(I,J) + 1 5547. C DO 210 L=LMU(I,J),1,-1 5548. C SDRAGU(I,J,L) = 0. 5549. C IF(LMM(I,J-1).lt.L) then 5550. C SDRAGU(I,J,L) = SDRAGU(I,J,L) + .5*DTS*SDRAGC/DYP(J) 5551. C LMINU(I,J) = L 5552. C endif 5553. C IF(LMM(IP1,J-1).lt.L) then 5554. C SDRAGU(I,J,L) = SDRAGU(I,J,L) + .5*DTS*SDRAGC/DYP(J) 5555. C LMINU(I,J) = L 5556. C endif 5557. C IF(LMM(I,J+1).lt.L) then 5558. C SDRAGU(I,J,L) = SDRAGU(I,J,L) + .5*DTS*SDRAGC/DYP(J) 5559. C LMINU(I,J) = L 5560. C endif 5561. C IF(LMM(IP1,J+1).lt.L) then 5562. C SDRAGU(I,J,L) = SDRAGU(I,J,L) + .5*DTS*SDRAGC/DYP(J) 5563. C LMINU(I,J) = L 5564. C endif 5565. C 210 continue 5566. C 220 I=IP1 5567. C**** TAUV = north-south side drag stress = SDRAGC*RHOO*|VO|*VO 5568. C**** VO = VO - DTS*TAUV/RHOO*DXV = 5569. C**** = VO - DTS*SDRAGC*|VO|*VO/DXV = 5570. C**** = VO * (1 - DTS*SDRAGC*|VO|/DXV) ~ 5571. C**** ~ VO / (1 + DTS*SDRAGC*|VO|/DXV) 5572. C IM1=IM-1 5573. C I=IM 5574. C DO 240 J=1,JM-1 5575. C DO 240 IP1=1,IM 5576. C LMINV(I,J) = LMV(I,J) + 1 5577. C DO 230 L=LMV(I,J),1,-1 5578. C SDRAGV(I,J,L) = 0. 5579. C IF(LMM(IM1,J).lt.L) then 5580. C SDRAGV(I,J,L) = SDRAGV(I,J,L) + .5*DTS*SDRAGC/DXV(J) 5581. C LMINV(I,J) = L 5582. C endif 5583. C IF(LMM(IM1,J+1).lt.L) then 5584. C SDRAGV(I,J,L) = SDRAGV(I,J,L) + .5*DTS*SDRAGC/DXV(J) 5585. C LMINV(I,J) = L 5586. C endif 5587. C IF(LMM(IP1,J).lt.L) then 5588. C SDRAGV(I,J,L) = SDRAGV(I,J,L) + .5*DTS*SDRAGC/DXV(J) 5589. C LMINV(I,J) = L 5590. C endif 5591. C IF(LMM(IP1,J+1).lt.L) then 5592. C SDRAGV(I,J,L) = SDRAGV(I,J,L) + .5*DTS*SDRAGC/DXV(J) 5593. C LMINV(I,J) = L 5594. C endif 5595. C 230 continue 5596. C IM1=I 5597. C 240 I=IP1 5598. C**** 5599. C**** Reduce horizontal gradients of tracers 5600. C**** 5601. 300 DO 320 I=IM+1,IM*(JM-1) 5602. DO 310 L=LMINX(I,1),LMM(I,1) 5603. GXM(I,1,L) = GXM(I,1,L)*REDUCE 5604. 310 SXM(I,1,L) = SXM(I,1,L)*REDUCE 5605. DO 320 L=LMINY(I,1),LMM(I,1) 5606. GYM(I,1,L) = GYM(I,1,L)*REDUCE 5607. 320 SYM(I,1,L) = SYM(I,1,L)*REDUCE 5608. C**** 5609. C**** Apply side drag to ocean currents 5610. C**** 5611. C DO 410 I=IM+1,IM*(JM-1) 5612. C DO 410 L=LMINU(I,1),LMU(I,1) 5613. C 410 UO(I,1,L) = UO(I,1,L) / (1. + SDRAGU(I,1,L)*ABS(UO(I,1,L))) 5614. C DO 420 I=1,IM*(JM-1) 5615. C DO 420 L=LMINV(I,1),LMV(I,1) 5616. C 420 VO(I,1,L) = VO(I,1,L) / (1. + SDRAGV(I,1,L)*ABS(VO(I,1,L))) 5617. RETURN 5618. END 6000. 6001. SUBROUTINE OABFIL 6002. C**** 6003. C**** OABFIL applies 8-th order Alternating Binomial Filter to UO and VO 6004. C**** fields. The frequency and strength depends on NFILXO and NFILYO. 6005. C**** The south pole is assumed to be land. 6006. C**** 6007. INCLUDE 'C070.COM' 6008. REAL*8 X(IM),Y(IM) 6008.5 LOGICAL*4 QFILX,QFILY 6009. C**** 6010. C**** Filtering in the east-west direction 6011. C**** 6011.5 QFILX = .FALSE. 6012. IF(NFILXO.le.0) GO TO 300 6013. IF(MOD(IDT-IDT0,NFILXO).ne.0) GO TO 300 6013.5 QFILX = .TRUE. 6014. X4TO8 = 1./(NFILXO*4.d0**8) 6015. IM1=IM 6016. DO 250 L=1,LMO 6017. C**** Filter the U component 6018. DO 140 J=2,JM-1 6019. DO 110 I=1,IM 6020. X(I) = UO(I,J,L) 6021. 110 Y(I) = 0. 6022. DO 130 N=1,8 6023. DO 120 I=1,IM 6024. IF(L.le.LMU(IM1,J) .and. L.le.LMU(I,J)) Y(I) = X(I)-X(IM1) 6025. 120 IM1=I 6026. DO 130 I=1,IM 6027. X(IM1) = Y(I)-Y(IM1) 6028. 130 IM1=I 6029. DO 140 I=1,IM 6030. 140 UO(I,J,L) = UO(I,J,L) - X(I)*X4TO8 6031. C**** Filter the V component 6032. DO 240 J=1,JM-1 6033. DO 210 I=1,IM 6034. X(I) = VO(I,J,L) 6035. 210 Y(I) = 0. 6036. DO 230 N=1,8 6037. DO 220 I=1,IM 6038. IF(L.le.LMV(IM1,J) .and. L.le.LMV(I,J)) Y(IM1) = X(I)-X(IM1) 6039. 220 IM1=I 6040. DO 230 I=1,IM 6041. X(I) = Y(I)-Y(IM1) 6042. 230 IM1=I 6043. DO 240 I=1,IM 6044. 240 VO(I,J,L) = VO(I,J,L) - X(I)*X4TO8 6045. 250 continue 6046. C**** 6047. C**** Filtering in the north-south direction 6048. C**** 6048.5 300 QFILY = .FALSE. 6049. IF(NFILYO.le.0) GO TO 500 6050. IF(MOD(IDT-IDT0,NFILYO).ne.0) GO TO 500 6050.5 QFILY = .TRUE. 6051. Y4TO8 = 1./(NFILYO*4.d0**8) 6052. DO 450 L=1,LMO 6053. C**** Filter the U component 6054. DO 340 I=1,IM 6055. DO 310 J=1,JM 6056. X(J) = UO(I,J,L) 6057. 310 Y(J) = 0. 6058. DO 330 N=1,8 6059. DO 320 J=2,JM-1 6060. IF(L.le.LMU(I,J) .and. L.le.LMU(I,J+1)) Y(J) = X(J+1)-X(J) 6061. 320 continue 6062. DO 330 J=2,JM 6063. 330 X(J) = Y(J)-Y(J-1) 6064. DO 340 J=2,JM 6065. 340 UO(I,J,L) = UO(I,J,L) - X(J)*Y4TO8 6066. C**** Make U winds at north pole to be uniform 6067. C Y(1) = 0. 6068. DO 350 I=1,IM 6069. 350 Y(1) = Y(1) + UO(I,JM,L) 6070. DO 360 I=1,IM 6071. 360 UO(I,JM,L) = Y(1)/IM 6072. C**** Filter the V component 6073. DO 440 I=1,IM 6074. DO 410 J=2,JM 6075. X(J) = VO(I,J,L) 6076. 410 Y(J) = 0. 6077. DO 430 N=1,8 6078. DO 420 J=3,JM-1 6079. IF(L.le.LMV(I,J-1) .and. L.le.LMV(I,J)) Y(J) = X(J)-X(J-1) 6080. 420 continue 6081. DO 430 J=2,JM-1 6082. 430 X(J) = Y(J+1)-Y(J) 6083. DO 440 J=2,JM-1 6084. 440 VO(I,J,L) = VO(I,J,L) - X(J)*Y4TO8 6085. 450 continue 6092. 500 RETURN 6093. END 7000. 7001. SUBROUTINE OSTRES 7002. C**** 7003. C**** OSTRES applies the atmospheric surface stress over open ocean 7004. C**** and the sea ice stress to the layer 1 ocean velocities 7005. C**** 7006. INCLUDE 'C070.COM' 7007. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 7008. * W0(IM,JM,4), E0(IM,JM,4),E1(IM,JM,4), UAS(IM,JM),VAS(IM,JM), 7009. * DMUI(IM,JM),DMVI(IM,JM) 7011. C**** 7012. C**** FLUXCB DMUA(1) U momentum downward into open ocean (kg/m*s) 7013. C**** DMVA(1) V momentum downward into open ocean (kg/m*s) 7013.1 C**** DMUA(2,JM,1) polar atmo. mass slowed to zero (kg/m**2) 7014. C**** DMUI U momentum downward from sea ice (kg/m*s) 7015. C**** DMVI V momentum downward from sea ice (kg/m*s) 7015.1 C**** DMUI(2,JM) polar sea ice mass slowed to zero (kg/m**2) 7016. C**** 7017. C**** Surface stress is applied to U component 7018. C**** 7019. I=IM 7020. DO 10 J=2,JM-1 7021. DO 10 IP1=1,IM 7022. IF(LMU(I,J).gt.0.) UO(I,J,1) = UO(I,J,1) + 7023. * (DMUA(I ,J,1)*(1.-RSI(I ,J)) + DMUI(I,J)*RSI(I ,J) + 7024. * DMUA(IP1,J,1)*(1.-RSI(IP1,J)) + DMUI(I,J)*RSI(IP1,J)) / 7025. * (MO(I,J,1) + MO(IP1,J,1)) 7026. 10 I=IP1 7027. UO(1,JM,1) = UO(1,JM,1)*(1. - DMUA(2,JM,1)/MO(1,JM,1)) 7029. C**** 7030. C**** Surface stress is applied to V component 7031. C**** 7032. DO 20 J=2,JM-2 7033. DO 20 I=1,IM 7034. 20 IF(LMV(I,J).GT.0.) VO(I,J,1) = VO(I,J,1) + 7035. *((DMVA(I,J ,1)*(1.-RSI(I,J ))+DMVI(I,J)*RSI(I,J ))*DXYN(J) + 7036. * (DMVA(I,J+1,1)*(1.-RSI(I,J+1))+DMVI(I,J)*RSI(I,J+1))*DXYS(J+1)) 7037. * / (MO(I,J,1)*DXYN(J) + MO(I,J+1,1)*DXYS(J+1)) 7038. C**** Surface stress is applied to V component at the North Pole 7041. DO 30 I=1,IM 7042. 30 VO(I,JM-1,1) = VO(I,JM-1,1) + 7043. * ((DMVA(I,JM-1,1)*(1.-RSI(I,JM-1)) + DMVI(I,JM-1)*RSI(I,JM-1))* 7044. * DXYN(JM-1) + 7045. * (DMVA(1,JM,1)*COSI(I) - DMUA(1,JM,1)*SINI(I))*DXYS(JM)) / 7046. * (MO(I,JM-1,1)*DXYN(JM-1) + MO(I,JM,1)*DXYS(JM)) 7047. RETURN 7048. END 8000. 8001. SUBROUTINE OSOURC 8002. C**** 8003. C**** OSOURC adds the vertical fluxes into the ocean 8004. C**** 8005. INCLUDE 'C070.COM' 8006. COMMON /FLUXCB/ SROX(IM,JM), DMUA(IM,JM,2),DMVA(IM,JM,2), 8007. * W0(IM,JM,4), E0(IM,JM,4), E1(IM,JM,4) 8008. PARAMETER (RHOI=910., ACE2SI=182., FLEAD=.06d0, 8009. * RFRAC=.62d0, ZETA1=1.5, ZETA2=20., 8010. * YSI1=XSI1*ACE1I/(ACE1I+ACE2SI), YSI3=XSI3*ACE2SI/(ACE1I+ACE2SI), 8011. * YSI2=XSI2*ACE1I/(ACE1I+ACE2SI), YSI4=XSI4*ACE2SI/(ACE1I+ACE2SI)) 8012. DATA IFIRST /1/ 8013. C**** 8014. C**** FIXDCB FOCEAN Ocean fraction (1) 8015. C**** 8016. C**** OCENCB MO Mass of ocean (kg/m^2) 8017. C**** G0M Mean potential enthalpy of ocean (J) 8018. C**** S0M Mean salt of ocean (kg) 8019. C**** 8020. C**** SEAICB RSI Horizontal ratio of sea ice to ocean (1) 8021. C**** MSI Mass of sea ice (kg/m**2) 8022. C**** HSI Enthalpy minus latent heat of sea ice (J/m^2) 8023. C**** 8026. C**** FLUXCB SROX Solar radiation absorbed into open ocean (J/m**2) 8027. C**** DMUA U surface momentum stress from atmosphere (kg/m*s) 8028. C**** DMVA V surface momentum stress from atmosphere (kg/m*s) 8029. C**** W0(1) Dew received over open ocean (kg/m**2) 8030. C**** E0(1) Vertical flux received over open ocean (J/m**2) 8031. C**** E1(1) Vertical flux from sea ice (J/m**2) 8032. C**** 8032.9 EF(Z) = RFRAC*EXP(-Z/ZETA1) + (1.-RFRAC)*EXP(-Z/ZETA2) 8033.2 IF(IFIRST.eq.0) GO TO 100 8033.4 IFIRST=0 8033.5 C**** 8033.6 C**** Calculate the fraction of solar energy absorbed in each layer 8033.7 C**** 8033.8 E0TO1 = 1. - EF(ZE(1)) 8033.9 C E0TOI = 1. 8034. E1TO2 = EF(ZE(1)) - EF(ZE(2)) 8034.1 E1TOI = EF(ZE(1)) 8034.2 C E2TO3 = EF(ZE(2)) - EF(ZE(3)) 8034.3 E2TOI = EF(ZE(2)) 8034.6 C**** Calculate gradient of solar energy absorbed in each layer 8034.9 C ZCEN = ZE(1)*.5 8035. C EZ0TO1 = -6.*(ZCEN - (ZCEN-ZE(1))*EF(ZE(1)) 8035.1 C * - ZETA1*(1. - EXP(-ZE(1)/ZETA1))*RFRAC 8035.2 C * - ZETA2*(1. - EXP(-ZE(1)/ZETA2))*(1.-RFRAC)) / ZE(1) 8035.3 C EZ0TOI = EZ0TO1 + 2.*EF(ZE(1)) 8035.5 ZCEN = (ZE(1)+ZE(2))*.5 8035.6 EZ1TO2 = -6.*((ZCEN-ZE(1))*EF(ZE(1)) - (ZCEN-ZE(2))*EF(ZE(2)) 8035.7 * - ZETA1*(EXP(-ZE(1)/ZETA1) - EXP(-ZE(2)/ZETA1))*RFRAC 8035.8 * - ZETA2*(EXP(-ZE(1)/ZETA2) - EXP(-ZE(2)/ZETA2))*(1.-RFRAC)) 8035.9 * / (ZE(2)-ZE(1)) 8036. EZ1TOI = EZ1TO2 + 2.*EF(ZE(2)) 8036.2 ZCEN = (ZE(2)+ZE(3))*.5 8036.3 EZ2TO3 = -6.*((ZCEN-ZE(2))*EF(ZE(2)) - (ZCEN-ZE(3))*EF(ZE(3)) 8036.4 * - ZETA1*(EXP(-ZE(2)/ZETA1) - EXP(-ZE(3)/ZETA1))*RFRAC 8036.5 * - ZETA2*(EXP(-ZE(2)/ZETA2) - EXP(-ZE(3)/ZETA2))*(1.-RFRAC)) 8036.6 * / (ZE(3)-ZE(2)) 8036.7 EZ2TOI = EZ2TO3 + 2.*EF(ZE(3)) 8050. C**** 8051. C**** Add surface source of fresh water and heat 8052. C**** 8053. 100 DO 620 J=2,JM 8054. IMAX=IM 8055. IF(J.eq.JM) IMAX=1 8063. DO 610 I=1,IMAX 8064. IF(FOCEAN(I,J).lt.1.) GO TO 610 8065. FOOCN = 1. - RSI(I,J) 8066. FSICE = RSI(I,J) 8066.05 GZM(I,J,1) = 0. 8066.06 SZM(I,J,1) = 0. 8066.07 IF(FSICE.eq.1.) GO TO 300 8066.08 MOO = MO(I,J,1) + W0(I,J,1) 8066.09 GMOO = G0M(I,J,1)*BYDXYP(J) + E0(I,J,1) 8066.1 C**** 8066.11 C**** Add insolation to the lower layers of the ocean 8066.12 C**** 8066.13 GO TO (200,120,130),LMM(I,J) 8066.14 C**** More than 3 layers in the ocean 8066.18 GZM(I,J,3) = GZM(I,J,3) + SROX(I,J)*DXYP(J)*FOOCN*EZ2TO3 8066.19 GO TO 140 8066.2 C**** Exactly 2 layers in the ocean 8066.22 120 GZM(I,J,2) = GZM(I,J,2) + SROX(I,J)*DXYP(J)*FOOCN*EZ1TOI 8066.23 G0M(I,J,2) = G0M(I,J,2) + SROX(I,J)*DXYP(J)*FOOCN*E1TOI 8066.24 GO TO 150 8066.3 C**** Exactly 3 layers in the ocean 8066.34 130 GZM(I,J,3) = GZM(I,J,3) + SROX(I,J)*DXYP(J)*FOOCN*EZ2TOI 8066.35 140 G0M(I,J,3) = G0M(I,J,3) + SROX(I,J)*DXYP(J)*FOOCN*E2TOI 8066.36 GZM(I,J,2) = GZM(I,J,2) + SROX(I,J)*DXYP(J)*FOOCN*EZ1TO2 8066.37 G0M(I,J,2) = G0M(I,J,2) + SROX(I,J)*DXYP(J)*FOOCN*E1TO2 8066.38 150 GMOO = GMOO + SROX(I,J)*(E0TO1-1.) 8067. C**** 8068. C**** Open Ocean 8069. C**** 8072. 200 GOO = GMOO/MOO 8073. SOO = S0M(I,J,1)*BYDXYP(J)/MOO 8074. GFOO = GFREZS(SOO) 8075. IF(GOO.lt.GFOO) GO TO 210 8076. IF(FSICE.gt.0.) GO TO 300 8077. C**** Update ocean prognostic variables 8078. MO(I,J,1) = MOO 8079. G0M(I,J,1) = GMOO*DXYP(J) 8080. GO TO 600 8081. C**** Open ocean is below freezing, calculate 8082. C**** DMOO = mass of ocean that freezes over open fraction from 8083. C**** GOO*MOO = GFOO*(MOO-DMOO) + (TFOO*SHCI-ELHM)*DMOO 8084. 210 TFOO = TFREZS(SOO) 8085. DMOO = MOO*(GOO-GFOO)/(TFOO*SHCI-ELHM-GFOO) ! > 0 8086. DEOO = (TFOO*SHCI-ELHM)*DMOO 8087. C EJ(J,37) = EJ(J,37) - DMOO*FOOCN ! in DIAGJ 8088. C EJ(J,46) = EJ(J,46) - DEOO*FOOCN ! in DIAGJ 8088.1 FJ(J,37) = FJ(J,37) + DMOO*FOOCN 8088.2 FJ(J,46) = FJ(J,46) + DEOO*FOOCN 8089. IF(FSICE.gt.0.) GO TO 310 8090. C**** Update ocean prognostic variables 8091. MO(I,J,1) = MOO-DMOO 8092. G0M(I,J,1) = (GMOO-DEOO)*DXYP(J) 8093. C**** Form new ice at surface 8094. RSI(I,J) = DMOO/(ACE1I+ACE2SI) 8094.4 RSIX(I,J) = 0. 8094.5 RSIY(I,J) = 0. 8095. MSI(I,J,1) = ACE1I 8096. MSI(I,J,2) = ACE2SI 8097. HSI(I,J,1) = ACE1I *XSI1*(DEOO/DMOO) 8098. HSI(I,J,2) = ACE1I *XSI2*(DEOO/DMOO) 8098.1 HSI(I,J,3) = ACE2SI*XSI3*(DEOO/DMOO) 8098.2 HSI(I,J,4) = ACE2SI*XSI4*(DEOO/DMOO) 8098.3 PSI(I,J) = 0. 8098.4 AIJ(I,J,61) = AIJ(I,J,61) + DMOO 8098.5 AIJ(I,J,62) = AIJ(I,J,62) + DEOO 8099. GO TO 600 8100. C**** 8101. C**** Ocean underneath the ice 8102. C**** 8103. 300 DMOO = 0. 8103.1 DEOO = 0. 8104. 310 MOI = MO(I,J,1) 8105. GMOI = G0M(I,J,1)*BYDXYP(J) + E1(I,J,1) 8106. GOI = GMOI/MOI 8107. SOI = S0M(I,J,1)*BYDXYP(J)/MOI 8108. GFOI = GFREZS(SOI) 8109. IF(GOI.GE.GFOI) GO TO 420 8110. C**** Ocean underneath the ice is below freezing, calculate 8111. C**** DMOI = mass of ocean that freezes under sea ice fraction from 8112. C**** GOI*MOI = GFOI*(MOI-DMOI) + (TFOI*SHCI-ELHM)*DMOI 8113. TFOI = TFREZS(SOI) 8114. DMOI = MOI*(GOI-GFOI)/(TFOI*SHCI-ELHM-GFOI) 8115. DEOI = (TFOI*SHCI-ELHM)*DMOI 8116. C EJ(J,37) = EJ(J,37) - DMOI*FSICE ! in DIAGJ 8117. C EJ(J,46) = EJ(J,46) - DEOI*FSICE ! in DIAGJ 8117.1 FJ(J,37) = FJ(J,37) + DMOI*FSICE 8117.2 FJ(J,46) = FJ(J,46) + DEOI*FSICE 8117.3 AIJ(I,J,61) = AIJ(I,J,61) + DMOO*FOOCN + DMOI*FSICE 8117.4 AIJ(I,J,62) = AIJ(I,J,62) + DEOO*FOOCN + DEOI*FSICE 8117.5 C**** Calculate advective heat flux from layer 3 to layer 4 of sea ice 8117.6 C FMSI3 = - XSI3*DMOI ! < 0 8117.7 FHSI3 = - HSI(I,J,4)*DMOI*(XSI3/XSI4) / MSI(I,J,2) 8118. C**** 8119. C**** Combine open ocean and sea ice fractions to form new prognostic 8120. C**** variables 8121. C**** 8122. IF(DMOO.gt.0.) GO TO 410 8123. C**** New sea ice is formed below old sea ice 8124. MO(I,J,1) = MOO*FOOCN + (MOI-DMOI)*FSICE 8125. G0M(I,J,1) = (GMOO*FOOCN + (GMOI-DEOI)*FSICE)*DXYP(J) 8126. MSI(I,J,2) = MSI(I,J,2) + DMOI 8127. HSI(I,J,3) = HSI(I,J,3) - FHSI3 8127.1 HSI(I,J,4) = HSI(I,J,4) + (FHSI3 + DEOI) 8128. GO TO 500 8129. C**** New sea ice is formed below old sea ice and on open ocean 8130. 410 MO(I,J,1) = (MOO-DMOO)*FOOCN + (MOI-DMOI)*FSICE 8131. G0M(I,J,1) = ((GMOO-DEOO)*FOOCN + (GMOI-DEOI)*FSICE)*DXYP(J) 8133. DRSI = DMOO*FOOCN / (ACE1I+ACE2SI) 8134. MSI(I,J,1) = (DRSI*ACE1I + FSICE* MSI(I,J,1)) / (FSICE+DRSI) 8135. MSI(I,J,2) = (DRSI*ACE2SI + FSICE*(MSI(I,J,2)+DMOI))/(FSICE+DRSI) 8137. HSI(I,J,1) = (FOOCN*DEOO*YSI1 + FSICE* HSI(I,J,1)) / (FSICE+DRSI) 8138. HSI(I,J,2) = (FOOCN*DEOO*YSI2 + FSICE* HSI(I,J,2)) / (FSICE+DRSI) 8139. HSI(I,J,3) = (FOOCN*DEOO*YSI3 + FSICE*(HSI(I,J,3)-FHSI3)) / 8140. / (FSICE+DRSI) 8140.1 HSI(I,J,4) = (FOOCN*DEOO*YSI4 + FSICE*(HSI(I,J,4)+FHSI3+DEOI)) / 8140.2 / (FSICE+DRSI) 8140.4 RSIX(I,J) = RSIX(I,J)*((FOOCN-DRSI)/FOOCN) 8140.5 RSIY(I,J) = RSIY(I,J)*((FOOCN-DRSI)/FOOCN) 8141. RSI(I,J) = RSI(I,J) + DRSI 8142. GO TO 500 8143. 420 IF(DMOO.gt.0.) GO TO 430 8144. C**** No new sea ice is formed 8145. MO(I,J,1) = MOO*FOOCN + MOI*FSICE 8146. G0M(I,J,1) = (GMOO*FOOCN + GMOI*FSICE)*DXYP(J) 8147. GO TO 500 8148. C**** New sea ice is formed on open ocean 8149. 430 MO(I,J,1) = (MOO-DMOO)*FOOCN + MOI*FSICE 8150. G0M(I,J,1) = ((GMOO-DEOO)*FOOCN + GMOI*FSICE)*DXYP(J) 8151. DRSI = DMOO*FOOCN / (ACE1I+ACE2SI) 8152. MSI(I,J,1) = (FSICE*MSI(I,J,1) + DRSI*ACE1I ) / (FSICE+DRSI) 8153. MSI(I,J,2) = (FSICE*MSI(I,J,2) + DRSI*ACE2SI) / (FSICE+DRSI) 8154. HSI(I,J,1) = (FSICE*HSI(I,J,1) + FOOCN*DEOO*YSI1) / (FSICE+DRSI) 8155. HSI(I,J,2) = (FSICE*HSI(I,J,2) + FOOCN*DEOO*YSI2) / (FSICE+DRSI) 8156. HSI(I,J,3) = (FSICE*HSI(I,J,3) + FOOCN*DEOO*YSI3) / (FSICE+DRSI) 8157. HSI(I,J,4) = (FSICE*HSI(I,J,4) + FOOCN*DEOO*YSI4) / (FSICE+DRSI) 8157.4 RSIX(I,J) = RSIX(I,J)*((FOOCN-DRSI)/FOOCN) 8157.5 RSIY(I,J) = RSIY(I,J)*((FOOCN-DRSI)/FOOCN) 8158. RSI(I,J) = RSI(I,J) + DRSI 8159. C**** 8160. C**** Compress the sea ice horizontally 8161. C**** 8162. 500 IF(MSI(I,J,2).ge.ACE2SI) GO TO 510 8163. C**** Sea ice is too thin 8164. RSIN = RSI(I,J)*(ACE1I+MSI(I,J,2))/(ACE1I+ACE2SI) 8165. GO TO 520 8166. 510 IF((1.-RSI(I,J))*(ACE1I+MSI(I,J,2)) .ge. FLEAD*RHOI) GO TO 600 8167. C**** Too little open ocean 8168. RSIN = 1. - FLEAD*RHOI/(ACE1I+MSI(I,J,2)) 8169. 520 DRSI = RSIN - RSI(I,J) ! < 0 8170. FMSI4 = (MSI(I,J,1)+MSI(I,J,2))*(DRSI/RSIN) ! < 0 8170.1 C FMSI3 = XSI3*FMSI4 ! < 0 8170.2 FHSI3 = HSI(I,J,4)*FMSI4*(XSI3/XSI4) / MSI(I,J,2) 8170.3 FHSI4 = (HSI(I,J,1)+HSI(I,J,2)+HSI(I,J,3)+HSI(I,J,4))*(DRSI/RSIN) 8171. HSI(I,J,3) = HSI(I,J,3) - FHSI3 8171.1 HSI(I,J,4) = HSI(I,J,4) + (FHSI3 - FHSI4) 8171.2 MSI(I,J,2) = MSI(I,J,2) - FMSI4 8172. RSI(I,J) = RSIN 8172.1 IF(RSIX(I,J).gt. RSIN) RSIX(I,J) = RSIN 8172.2 IF(RSIX(I,J).lt.-RSIN) RSIX(I,J) = -RSIN 8172.3 IF(RSIY(I,J).gt. RSIN) RSIY(I,J) = RSIN 8172.4 IF(RSIY(I,J).lt.-RSIN) RSIY(I,J) = -RSIN 8175. C**** Accumulate diagnostics 8175.5 600 EJ(J,18) = EJ(J,18) + FOOCN*TEMGS(GOO,SOO) 8176. EJ(J,19) = EJ(J,19) + FOOCN*W0(I,J,1) 8176.1 C EJ(J,22) = EJ(J,22) + FOOCN*W0(I,J,1)*ZATMO(I,J) ! Caspian 8178. AIJ(I,J,53) = AIJ(I,J,53) - FOOCN*W0(I,J,1) 8179. AIJ(I,J,57) = AIJ(I,J,57) + FOOCN*E0(I,J,1) 8180. 610 continue 8188. 620 continue 8207. RETURN 8208. END 9000. 9001. SUBROUTINE PRECOO 9002. C**** 9003. C**** PRECOO adds precipitation to the open ocean 9004. C**** 9005. INCLUDE 'C070.COM' 9006. COMMON /PRECCB/ PREC(IM,JM,2),EPRE(IM,JM,2) 9007. C**** 9008. C**** FIXDCB FOCEAN Ocean fraction (0 or 1) 9009. C**** 9010. C**** OCENCB MO Mass of ocean (kg/m^2) 9011. C**** G0M Mean potential enthalpy of ocean (J) 9012. C**** 9013. C**** SEAICB RSI Horizontal ratio of sea ice to ocean (1) 9014. C**** 9015. C**** FLUXCB PREC Precipitation from atmosphere (kg/m^2) 9016. C**** EPRE Energy of precipitation (J/m^2) 9017. C**** 9018. DO 10 J=2,JM 9019. IMAX=IM 9020. IF(J.eq.JM) IMAX=1 9021. DO 10 I=1,IMAX 9022. IF(FOCEAN(I,J).le.0. .or. PREC(I,J,1).le.0.) GO TO 10 9023. MO(I,J,1) = MO(I,J,1) + PREC(I,J,1)*(1.-RSI(I,J)) 9024. G0M(I,J,1) = G0M(I,J,1) + EPRE(I,J,1)*(1.-RSI(I,J))*DXYP(J) 9025. 10 continue 9026. RETURN 9027. END