1. C**** 2. C**** C088K.S Fortran source code for KPP mixing scheme 98/09/10 3. C**** 4. C**** NASA/GISS Climate Model III Gavin A. Schmidt 5. C**** 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**** 1001. C**** KPPMIX.F contains self contained code for one dimensional 1002. C**** calculation of Large et al (1996) KPP scheme (no slab calcs) 1003. C**** Gavin Schmidt Feb 1998 1004. 1005. SUBROUTINE KVINIT 1006. 1007. C**** Initialise KMIX and save pre-source term surface values of 1008. C**** enthalpy, mass and horizontal gradients for kppmix calculation 1009. 1010. INCLUDE 'C070.COM' 1011. COMMON /DIFVCB/ G0M1(IM,JM,3),MO1(IM,JM),GXM1(IM,JM),GYM1(IM,JM) 1012. * ,SXM1(IM,JM),SYM1(IM,JM),UO1(IM,JM),VO1(IM,JM) 1013. 1014. C**** Save surface values 1015. DO J=1,JM 1016. DO I=1,IM 1017. G0M1(I,J,1:3) = G0M(I,J,1:3) 1018. GXM1(I,J) = GXM(I,J,1) 1019. GYM1(I,J) = GYM(I,J,1) 1020. SXM1(I,J) = SXM(I,J,1) 1021. SYM1(I,J) = SYM(I,J,1) 1022. MO1(I,J) = MO(I,J,1) 1023. UO1(I,J) = UO(I,J,1) 1024. VO1(I,J) = VO(I,J,1) 1025. END DO 1026. END DO 1027. RETURN 1028. END 1029. 1030. subroutine kppmix(LDD , km , 1031. $ zgrid , hwide , kmtj , Shsq , dVsq , 1032. $ ustar , Bo , Bosol , alphaDT, betaDS, 1033. $ dbloc , Ritop , Coriol, byhwide, 1034. $ visc , difs , dift , ghats , hbl 1035. $ , kbl ) 1036. 1037. c Main driver subroutine for kpp vertical mixing scheme and 1038. c interface to greater ocean model 1039. c 1040. c written by: bill large, june 6, 1994 1041. c modified by: jan morzel, june 30, 1994 1042. c bill large, august 11, 1994 1043. c bill large, january 25, 1995 : "dVsq" and 1d code 1044. c modified for GISS by Gavin Schmidt, march 1998 1045. c 1046. c 1047. INCLUDE "KPP_1D.COM" 1048. c 1049. parameter (mdiff = 3) ! number of diffusivities for local arrays 1050. c input 1051. real*8 zgrid(0:km+1) ! vertical grid (<= 0) (m) 1052. real*8 hwide(0:km+1) ! layer thicknesses (m) 1053. real*8 byhwide(0:km+1) ! 1/layer thicknesses (1/m) 1054. integer kmtj ! number of vertical layers on this row 1055. real*8 Shsq(km) ! (local velocity shear)^2 (m/s)^2 1056. real*8 dVsq(km) ! (velocity shear re sfc)^2 (m/s)^2 1057. real*8 ustar ! surface friction velocity (m/s) 1058. real*8 Bo ! surface turbulent buoy. forcing (m^2/s^3) 1059. real*8 Bosol ! radiative buoyancy forcing (m^2/s^3) 1060. real*8 alphaDT(km) ! alpha * DT across interfaces (kg/m^3) 1061. real*8 betaDS(km) ! beta * DS across interfaces (kg/m^3) 1062. real*8 dbloc(km) ! local delta buoyancy (m/s^2) 1063. c ! across interfaces 1064. real*8 Ritop(km) ! numerator of bulk Richardson 1065. c Number: (m/s)^2 1066. c Ritop = (-z - -zref)* delta buoyancy w/ respect to sfc 1067. real*8 Coriol ! Coriolis parameter (1/s) 1068. logical LDD ! = TRUE for double diffusive parameterisation 1069. c output 1070. real*8 visc(0:km+1) ! vertical viscosity coefficient (m^2/s) 1071. real*8 difs(0:km+1) ! vertical scalar diffusivity (m^2/s) 1072. real*8 dift(0:km+1) ! vertical temperature diffusivity(m^2/s) 1073. real*8 ghats(km) ! nonlocal transport (s/m^2) 1074. real*8 hbl ! boundary layer depth (m) 1075. real*8 byhbl ! 1/boundary layer depth (1/m) 1076. c local 1077. real*8 bfsfc ! surface buoyancy forcing (m^2/s^3) 1078. real*8 ws ! momentum velocity scale 1079. real*8 wm ! scalar velocity scale 1080. real*8 caseA ! = 1 in case A; =0 in case B 1081. real*8 stable ! = 1 in stable forcing; =0 in unstable 1082. real*8 dkm1(mdiff) ! boundary layer difs at kbl-1 level 1083. real*8 gat1(mdiff) ! shape function at sigma=1 1084. real*8 dat1(mdiff) ! derivative of shape function at sigma=1 1085. real*8 blmc(km,mdiff)! boundary layer mixing coefficients 1086. real*8 sigma ! normalized depth (d / hbl) 1087. real*8 Rib(2) ! bulk Richardson number 1088. integer kbl ! index of first grid level below hbl 1089. integer kmax ! minimum of ksrpd and kmtj, used in swfrac 1090. c local ri_iwmix 1091. real*8 Rigg ! local richardson number 1092. real*8 fri,fcon ! function of Rig 1093. c local enhance 1094. real*8 delta ! fraction hbl lies beteen zgrid neighbors 1095. c local wscale 1096. real*8 zehat ! = zeta * ustar**3 1097. INTENT (IN) LDD,km,zgrid,hwide,kmtj,Shsq,dVsq,ustar,Bo,Bosol 1098. * ,alphaDT,betaDS,dbloc,Ritop,Coriol,byhwide 1099. INTENT (OUT) visc, difs, dift, ghats, hbl, kbl 1100. 1101. kmax = MIN(ksrpd,kmtj) 1102. c compute interior mixing coefficients everywhere, due to constant 1103. c internal wave activity, static instability, and local shear 1104. c instability. 1105. c Zero surface values 1106. 1107. c call ri_iwmix ( km, km+1, kmtj, Shsq, dbloc, zgrid, 1108. c * visc, difs , dift ) 1109. c compute interior viscosity diffusivity coefficients due 1110. c to shear instability (dependent on a local richardson number), 1111. c to background internal wave activity, and 1112. c to static instability (local richardson number < 0). 1113. 1114. c compute interior gradient Ri at all interfaces ki=1,km, (not surfa 1115. ce) 1116. c use visc(ki=1,km) as temporary storage to be smoothed 1117. c use dift(ki=1,km) as temporary storage of N^2 = BVSQ 1118. c set values at bottom and below to nearest value above bottom 1119. 1120. do ki = 1, kmtj 1121. visc(ki) = dbloc(ki) * (zgrid(ki)-zgrid(ki+1)) / 1122. $ ( Shsq(ki) + epsl) 1123. dift(ki) = dbloc(ki) * byhwide(ki+1) 1124. end do 1125. 1126. c vertically smooth Ri num_v_smooth_Ri times 1127. 1128. do mr = 1,num_v_smooth_Ri 1129. call z121(visc,kmtj,km) 1130. end do 1131. 1132. do ki = 1, kmtj 1133. 1134. c evaluate f of Vaisala squared for convection store in fcon 1135. c evaluate f of smooth Ri (fri) for shear instability store in fri 1136. 1137. Rigg = DMAX1( dift(ki) , BVSQcon ) 1138. ratio = DMIN1( (BVSQcon-Rigg) * rBVSQcon , c1 ) 1139. fcon = (c1 - ratio*ratio) 1140. fcon = fcon * fcon * fcon 1141. 1142. Rigg = DMAX1( visc(ki) , c0 ) 1143. ratio = DMIN1( Rigg * rRiinfty , c1 ) 1144. fri = (c1 - ratio*ratio) 1145. fri = fri * fri * fri 1146. 1147. c evaluate diffusivities and viscosity 1148. c mixing due to internal waves, and shear and static instability 1149. 1150. visc(ki) = (difmiw + fcon * difmcon + fri * difm0) 1151. difs(ki) = (difsiw + fcon * difscon + fri * difs0) 1152. dift(ki) = difs(ki) 1153. end do 1154. 1155. c set surface values to 0.0 1156. 1157. visc(0) = c0 1158. dift(0) = c0 1159. difs(0) = c0 1160. 1161. c add double diffusion if desired 1162. 1163. if (LDD) call ddmix (km,alphaDT,betaDS,visc,difs,dift,kmtj) 1164. 1165. c Zero values at seafloor and below for blmix 1166. 1167. visc(kmtj:km+1) = c0 1168. difs(kmtj:km+1) = c0 1169. dift(kmtj:km+1) = c0 1170. 1171. c compute boundary layer mixing coefficients: 1172. c diagnose the new boundary layer depth 1173. 1174. c call bldepth (km, km+1, zgrid, hwide, kmtj, dVsq, 1175. c $ dbloc, Ritop, ustar, Bo, Bosol, Coriol, 1176. c $ hbl, bfsfc, stable, caseA, kbl, 1177. c $ Rib, sigma, wm, ws) 1178. 1179. c the oceanic planetray boundary layer depth, hbl, is determined as 1180. c the shallowest depth where the bulk richardson number is 1181. c equal to the critical value, Ricr. 1182. c Bulk richardson numbers are evaluated by computing velocity and 1183. c buoyancy differences between values at zgrid(kl) < 0 and surface 1184. c reference values. 1185. c in this configuration, the reference values are equal to the 1186. c values in the surface layer. 1187. c when using a very fine vertical grid, these values should be 1188. c computed as the vertical average of velocity and buoyancy from 1189. c the surface down to epsilon*zgrid(kl). 1190. c When the bulk richardson number at k exceeds Ricr, hbl is 1191. c linearly interpolated between grid levels zgrid(k) and zgrid(k-1). 1192. c The water column and the surface forcing are diagnosed for 1193. c stable/ustable forcing conditions, and where hbl is relative 1194. c to grid points (caseA), so that conditional branches can be 1195. c avoided in later subroutines. 1196. 1197. c find bulk Richardson number at every grid level until > Ricr 1198. c 1199. c note: the reference depth is -epsilon/2.*zgrid(k), but the reference 1200. c u,v,t,s values are simply the surface layer values, 1201. c and not the averaged values from 0 to 2*ref.depth, 1202. c which is necessary for very fine grids(top layer < 2m thickness) 1203. c note: max values when Ricr never satisfied are 1204. c kbl=kmtj and hbl=-zgrid(kmtj) 1205. 1206. c indices for array Rib(k), the bulk Richardson number. 1207. ka = 1 1208. ku = 2 1209. 1210. c initialize hbl and kbl to bottomed out values 1211. Rib(ka) = 0.0 1212. kbl = kmtj 1213. hbl = -zgrid(kbl) 1214. 1215. kl = 1 1216. 10 kl = kl + 1 1217. 1218. c compute bfsfc = sw fraction at hbf * zgrid 1219. c call swfrac(zgrid(kl),kl,kmtj,kmax,bfsfc) ! inlined 1220. 1221. C**** Get actual k (since k is on staggered grid) 1222. kt = kl + NINT(0.5d0 + SIGN(0.5d0,-(ZE(kl)+zgrid(kl)))) 1223. 1224. C**** calculate fraction 1225. if (kt.gt.kmax) then 1226. bfsfc = 0. 1227. else 1228. if (kt.eq.kmax) then 1229. bfsfc = E0TO(kt) + (zgrid(kl)+ZE(kt-1))*EZTOI(kt) 1230. else 1231. bfsfc = E0TO(kt) + (zgrid(kl)+ZE(kt-1))*EZTO(kt) 1232. end if 1233. end if 1234. 1235. c use caseA as temporary array 1236. caseA = -zgrid(kl) 1237. 1238. c compute bfsfc= Bo + radiative contribution down to hbf * hbl 1239. bfsfc = Bo + Bosol * (1. - bfsfc) 1240. 1241. stable = 0.5 + SIGN( 5d-1, bfsfc ) 1242. sigma = stable * 1. + (1.-stable) * epsilon 1243. 1244. c compute velocity scales at sigma, for hbl= caseA = -zgrid(kl) 1245. c call wscale(sigma, caseA, ustar, bfsfc, wm, ws) 1246. c compute turbulent velocity scales. 1247. c use a 2D-lookup table for wm and ws as functions of ustar and 1248. c zetahat (=vonk*sigma*hbl*bfsfc). 1249. c use lookup table for zehat < zmax only; otherwise use 1250. c stable formulae 1251. 1252. zehat = vonk * sigma * caseA * bfsfc 1253. 1254. if (zehat.le.zmax) then 1255. zdiff = zehat-zmin 1256. iz = int( zdiff * rdeltaz ) 1257. iz = min( iz , nni ) 1258. iz = max( iz , 0 ) 1259. izp1=iz+1 1260. 1261. udiff = ustar-umin 1262. ju = int( udiff * rdeltau) 1263. ju = min( ju , nnj ) 1264. ju = max( ju , 0 ) 1265. jup1=ju+1 1266. 1267. zfrac = zdiff*rdeltaz - float(iz) 1268. ufrac = udiff*rdeltau - float(ju) 1269. 1270. fzfrac= 1.-zfrac 1271. wam = (fzfrac) * wmt(iz,jup1) + zfrac*wmt(izp1,jup1) 1272. wbm = (fzfrac) * wmt(iz,ju ) + zfrac*wmt(izp1,ju ) 1273. wm = (1.-ufrac)* wbm + ufrac*wam 1274. 1275. was = (fzfrac) * wst(iz,jup1) + zfrac*wst(izp1,jup1) 1276. wbs = (fzfrac) * wst(iz,ju ) + zfrac*wst(izp1,ju ) 1277. ws = (1.-ufrac)* wbs + ufrac*was 1278. else 1279. u3 = ustar*ustar*ustar 1280. wm = vonk * ustar * u3 / ( u3 + conc1*zehat ) 1281. ws = wm 1282. endif 1283. 1284. c compute the turbulent shear contribution to Rib 1285. bvsq =0.5* 1286. $ ( dbloc(kl-1) * byhwide(kl) + 1287. $ dbloc(kl ) * byhwide(kl+1) ) 1288. Vtsq = - zgrid(kl) * ws * sqrt(abs(bvsq)) * Vtc 1289. 1290. c compute bulk Richardson number at new level, dunder 1291. c note: Ritop needs to be zero on land and ocean bottom 1292. c points so that the following if statement gets triggered 1293. c correctly. otherwise, hbl might get set to (big) negative 1294. c values, that might exceed the limit for the "exp" function 1295. c in "swfrac" 1296. 1297. Rib(ku) = Ritop(kl) / (dVsq(kl)+Vtsq+epsl) 1298. 1299. c linearly interpolate to find hbl where Rib = Ricr 1300. if((kbl.eq.kmtj).and.(Rib(ku).gt.Ricr)) then 1301. hbl = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) * 1302. $ (Ricr - Rib(ka)) / (Rib(ku)-Rib(ka)) 1303. kbl = kl 1304. else 1305. ksave = ka 1306. ka = ku 1307. ku = ksave 1308. if (kl.lt.kmtj) go to 10 1309. end if 1310. 1311. c find stability and buoyancy forcing for boundary layer 1312. c call swfrac(-hbl,kbl-1,kmtj,bfsfc,kmax) ! inlined 1313. 1314. C**** Get actual k (since k is on staggered grid) 1315. kt = kbl-1+ NINT(0.5d0 + SIGN(0.5d0,-(ZE(kbl-1)-hbl))) 1316. 1317. C**** calculate fraction 1318. if (kt.gt.kmax) then 1319. bfsfc = 0. 1320. else 1321. if (kt.eq.kmax) then 1322. bfsfc = E0TO(kt) + (ZE(kt-1)-hbl)*EZTOI(kt) 1323. else 1324. bfsfc = E0TO(kt) + (ZE(kt-1)-hbl)*EZTO(kt) 1325. end if 1326. end if 1327. 1328. bfsfc = Bo + Bosol * (1. - bfsfc) 1329. stable = 0.5 + SIGN( 5d-1, bfsfc ) 1330. bfsfc = bfsfc + stable * epsl !ensures bfsfc never=0 1331. 1332. c check hbl limits for hekman or hmonob (NOT USED) 1333. 1334. c if(bfsfc.gt.0.0) then 1335. c hekman = cekman * ustar / (abs(Coriol)+epsl) 1336. c hmonob = cmonob * ustar*ustar*ustar 1337. c $ /(vonk * (bfsfc+epsl) ) 1338. c hlimit = stable * DMIN1(hekman,hmonob) + 1339. c $ (stable-1.) * zgrid(km) 1340. c hbl = DMIN1(hbl,hlimit) 1341. c hbl = DMAX1(hbl,-zgrid(1)) 1342. c endif 1343. c kbl = kmtj 1344. 1345. c find new kbl 1346. c do kl=2,kmtj 1347. c if((kbl.eq.kmtj).and.(-zgrid(kl).gt.hbl)) kbl = kl 1348. c end do 1349. 1350. c find stability and buoyancy forcing for final hbl values 1351. c call swfrac(-hbl,kbl-1,kmtj,kmax,bfsfc) 1352. 1353. c bfsfc = Bo + Bosol * (1. - bfsfc) 1354. c stable = 0.5 + SIGN( 5d-1, bfsfc ) 1355. c bfsfc = bfsfc + stable * epsl 1356. 1357. c determine caseA and caseB 1358. caseA = 0.5 + SIGN(5d-1,-zgrid(kbl) - 0.5*hwide(kbl) - hbl) 1359. 1360. byhbl = 1d0/hbl 1361. c compute boundary layer diffusivities 1362. 1363. c call blmix (km , km+1 , mdiff , zgrid, hwide , 1364. c $ ustar, bfsfc, hbl , stable, caseA, 1365. c $ visc , difs , dift , kbl , 1366. c $ gat1 , dat1 , dkm1 , blmc , ghats, 1367. c $ sigma, wm , ws ) 1368. c mixing coefficients within boundary layer depend on surface 1369. c forcing and the magnitude and gradient of interior mixing below 1370. c the boundary layer ("matching"). 1371. c Caution: if mixing bottoms out at hbl = -zgrid(km) then 1372. c fictitious layer at km+1 is needed with small but finite width 1373. c hwide(km+1) (eg. epsl = 1.e-20). 1374. 1375. c compute velocity scales at hbl 1376. 1377. sigma = stable * 1.0 + (1.-stable) * epsilon 1378. 1379. c call wscale(sigma, hbl, ustar, bfsfc, wm, ws) 1380. c compute turbulent velocity scales. 1381. c use a 2D-lookup table for wm and ws as functions of ustar and 1382. c zetahat (=vonk*sigma*hbl*bfsfc). 1383. c use lookup table for zehat < zmax only; otherwise use 1384. c stable formulae 1385. 1386. zehat = vonk * sigma * hbl * bfsfc 1387. 1388. if (zehat.le.zmax) then 1389. zdiff = zehat-zmin 1390. iz = int( zdiff * rdeltaz ) 1391. iz = min( iz , nni ) 1392. iz = max( iz , 0 ) 1393. izp1=iz+1 1394. 1395. udiff = ustar-umin 1396. ju = int( udiff * rdeltau) 1397. ju = min( ju , nnj ) 1398. ju = max( ju , 0 ) 1399. jup1=ju+1 1400. 1401. zfrac = zdiff*rdeltaz - float(iz) 1402. ufrac = udiff*rdeltau - float(ju) 1403. 1404. fzfrac= 1.-zfrac 1405. wam = (fzfrac) * wmt(iz,jup1) + zfrac*wmt(izp1,jup1) 1406. wbm = (fzfrac) * wmt(iz,ju ) + zfrac*wmt(izp1,ju ) 1407. wm = (1.-ufrac)* wbm + ufrac*wam 1408. 1409. was = (fzfrac) * wst(iz,jup1) + zfrac*wst(izp1,jup1) 1410. wbs = (fzfrac) * wst(iz,ju ) + zfrac*wst(izp1,ju ) 1411. ws = (1.-ufrac)* wbs + ufrac*was 1412. else 1413. u3 = ustar*ustar*ustar 1414. wm = vonk * ustar * u3 / ( u3 + conc1*zehat ) 1415. ws = wm 1416. endif 1417. 1418. kn = int(caseA+epsl) *(kbl -1) + 1419. $ (1-int(caseA+epsl)) * kbl 1420. 1421. c find the interior viscosities and derivatives at hbl 1422. delhat = 0.5*hwide(kn) - zgrid(kn) - hbl 1423. R = 1.0 - delhat * byhwide(kn) 1424. dvdzup = (visc(kn-1) - visc(kn)) * byhwide(kn) 1425. dvdzdn = (visc(kn) - visc(kn+1)) * byhwide(kn+1) 1426. viscp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup))+ 1427. $ R * (dvdzdn + abs(dvdzdn)) ) 1428. 1429. dvdzup = (difs(kn-1) - difs(kn)) * byhwide(kn) 1430. dvdzdn = (difs(kn) - difs(kn+1)) * byhwide(kn+1) 1431. difsp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup))+ 1432. $ R * (dvdzdn + abs(dvdzdn)) ) 1433. 1434. dvdzup = (dift(kn-1) - dift(kn)) * byhwide(kn) 1435. dvdzdn = (dift(kn) - dift(kn+1)) * byhwide(kn+1) 1436. diftp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup))+ 1437. $ R * (dvdzdn + abs(dvdzdn)) ) 1438. 1439. visch = visc(kn) + viscp * delhat 1440. difsh = difs(kn) + difsp * delhat 1441. difth = dift(kn) + diftp * delhat 1442. 1443. f1 = stable * conc1 * bfsfc / (ustar**4+epsl) 1444. 1445. bywm = 1./(wm+epsl) 1446. byws = 1./(ws+epsl) 1447. 1448. gat1(1) = visch * byhbl * bywm 1449. dat1(1) = -viscp * bywm + f1 * visch 1450. dat1(1) = min(dat1(1),0d0) 1451. 1452. gat1(2) = difsh * byhbl * byws 1453. dat1(2) = -difsp * byws + f1 * difsh 1454. dat1(2) = min(dat1(2),0d0) 1455. 1456. gat1(3) = difth * byhbl * byws 1457. dat1(3) = -diftp * byws + f1 * difth 1458. dat1(3) = min(dat1(3),0d0) 1459. 1460. do ki = 1,kbl-1 1461. 1462. c compute turbulent velocity scales on the interfaces 1463. sig = (-zgrid(ki) + 0.5 * hwide(ki)) * byhbl 1464. sigma= stable*sig + (1.-stable)*DMIN1(sig,epsilon) 1465. 1466. c call wscale(sigma, hbl, ustar, bfsfc, wm, ws) 1467. c compute turbulent velocity scales. 1468. c use a 2D-lookup table for wm and ws as functions of ustar and 1469. c zetahat (=vonk*sigma*hbl*bfsfc). 1470. c use lookup table for zehat < zmax only; otherwise use 1471. c stable formulae 1472. 1473. zehat = vonk * sigma * hbl * bfsfc 1474. 1475. if (zehat.le.zmax) then 1476. zdiff = zehat-zmin 1477. iz = int( zdiff * rdeltaz ) 1478. iz = min( iz , nni ) 1479. iz = max( iz , 0 ) 1480. izp1=iz+1 1481. 1482. udiff = ustar-umin 1483. ju = int( udiff * rdeltau) 1484. ju = min( ju , nnj ) 1485. ju = max( ju , 0 ) 1486. jup1=ju+1 1487. 1488. zfrac = zdiff*rdeltaz - float(iz) 1489. ufrac = udiff*rdeltau - float(ju) 1490. 1491. fzfrac= 1.-zfrac 1492. wam = (fzfrac) * wmt(iz,jup1) + zfrac*wmt(izp1,jup1) 1493. wbm = (fzfrac) * wmt(iz,ju ) + zfrac*wmt(izp1,ju ) 1494. wm = (1.-ufrac)* wbm + ufrac*wam 1495. 1496. was = (fzfrac) * wst(iz,jup1) + zfrac*wst(izp1,jup1) 1497. wbs = (fzfrac) * wst(iz,ju ) + zfrac*wst(izp1,ju ) 1498. ws = (1.-ufrac)* wbs + ufrac*was 1499. else 1500. u3 = ustar*ustar*ustar 1501. wm = vonk * ustar * u3 / ( u3 + conc1*zehat ) 1502. ws = wm 1503. endif 1504. 1505. c compute the dimensionless shape functions at the interfaces 1506. sig = (-zgrid(ki) + 0.5 * hwide(ki)) * byhbl 1507. a1 = sig - 2. 1508. a2 = 3.-2.*sig 1509. a3 = sig - 1. 1510. 1511. Gm = a1 + a2 * gat1(1) + a3 * dat1(1) 1512. Gs = a1 + a2 * gat1(2) + a3 * dat1(2) 1513. Gt = a1 + a2 * gat1(3) + a3 * dat1(3) 1514. 1515. c compute boundary layer diffusivities at the interfaces 1516. blmc(ki,1) = hbl * wm * sig * (1. + sig * Gm) 1517. blmc(ki,2) = hbl * ws * sig * (1. + sig * Gs) 1518. blmc(ki,3) = hbl * ws * sig * (1. + sig * Gt) 1519. 1520. c nonlocal transport term = ghats * o 1521. ghats(ki) = (1.-stable) * cg * byws * byhbl 1522. end do 1523. 1524. c find diffusivities at kbl-1 grid level 1525. sig = -zgrid(kbl-1) * byhbl 1526. sigma =stable * sig + (1.-stable) * MIN(sig,epsilon) 1527. 1528. c call wscale(sigma, hbl, ustar, bfsfc, wm, ws) 1529. c compute turbulent velocity scales. 1530. c use a 2D-lookup table for wm and ws as functions of ustar and 1531. c zetahat (=vonk*sigma*hbl*bfsfc). 1532. c use lookup table for zehat < zmax only; otherwise use 1533. c stable formulae 1534. 1535. zehat = vonk * sigma * hbl * bfsfc 1536. 1537. if (zehat.le.zmax) then 1538. zdiff = zehat-zmin 1539. iz = int( zdiff * rdeltaz ) 1540. iz = min( iz , nni ) 1541. iz = max( iz , 0 ) 1542. izp1=iz+1 1543. 1544. udiff = ustar-umin 1545. ju = int( udiff * rdeltau) 1546. ju = min( ju , nnj ) 1547. ju = max( ju , 0 ) 1548. jup1=ju+1 1549. 1550. zfrac = zdiff*rdeltaz - float(iz) 1551. ufrac = udiff*rdeltau - float(ju) 1552. 1553. fzfrac= 1.-zfrac 1554. wam = (fzfrac) * wmt(iz,jup1) + zfrac*wmt(izp1,jup1) 1555. wbm = (fzfrac) * wmt(iz,ju ) + zfrac*wmt(izp1,ju ) 1556. wm = (1.-ufrac)* wbm + ufrac*wam 1557. 1558. was = (fzfrac) * wst(iz,jup1) + zfrac*wst(izp1,jup1) 1559. wbs = (fzfrac) * wst(iz,ju ) + zfrac*wst(izp1,ju ) 1560. ws = (1.-ufrac)* wbs + ufrac*was 1561. else 1562. u3 = ustar*ustar*ustar 1563. wm = vonk * ustar * u3 / ( u3 + conc1*zehat ) 1564. ws = wm 1565. endif 1566. 1567. sig = -zgrid(kbl-1) * byhbl 1568. a1= sig - 2. 1569. a2 = 3.-2.*sig 1570. a3 = sig - 1. 1571. Gm = a1 + a2 * gat1(1) + a3 * dat1(1) 1572. Gs = a1 + a2 * gat1(2) + a3 * dat1(2) 1573. Gt = a1 + a2 * gat1(3) + a3 * dat1(3) 1574. dkm1(1) = hbl * wm * sig * (1. + sig * Gm) 1575. dkm1(2) = hbl * ws * sig * (1. + sig * Gs) 1576. dkm1(3) = hbl * ws * sig * (1. + sig * Gt) 1577. 1578. c call enhance (km , km+1 , mdiff , dkm1 , visc , 1579. c $ difs , dift , hbl , kbl , zgrid, caseA , 1580. c $ blmc ,ghats ) 1581. c enhance the diffusivity at the kbl-.5 interface 1582. 1583. if (kbl.le.kmtj) then 1584. ki = kbl-1 1585. delta = (hbl+zgrid(ki)) * byhwide(ki+1) 1586. 1587. dkmp5 = caseA*visc(ki) + (1.-caseA)*blmc(ki,1) 1588. dstar = (1.-delta)**2 * dkm1(1) + delta**2 * dkmp5 1589. blmc(ki,1) = (1.-delta) * visc(ki) + delta * dstar 1590. 1591. dkmp5 = caseA*difs(ki) + (1.-caseA)*blmc(ki,2) 1592. dstar = (1.-delta)**2 * dkm1(2) + delta**2 * dkmp5 1593. blmc(ki,2) = (1.-delta) * difs(ki) + delta * dstar 1594. 1595. dkmp5 = caseA*dift(ki) + (1.-caseA)*blmc(ki,3) 1596. dstar = (1.-delta)**2 * dkm1(3) + delta**2 * dkmp5 1597. blmc(ki,3) = (1.-delta) * dift(ki) + delta * dstar 1598. 1599. ghats(ki) = (1.-caseA) * ghats(ki) 1600. end if 1601. 1602. c combine interior and boundary layer coefficients and nonlocal term 1603. 1604. visc(1:kbl-1)=blmc(1:kbl-1,1) 1605. difs(1:kbl-1)=blmc(1:kbl-1,2) 1606. dift(1:kbl-1)=blmc(1:kbl-1,3) 1607. ghats(kbl:km)=0. 1608. 1609. return 1610. end 1611. 1612. subroutine wscale(sigma, hbl, ustar, bfsfc, wm , ws) 1613. 1614. c compute turbulent velocity scales. 1615. c use a 2D-lookup table for wm and ws as functions of ustar and 1616. c zetahat (=vonk*sigma*hbl*bfsfc). 1617. c 1618. c note: the lookup table is only used for unstable conditions 1619. c (zehat.le.0), in the stable domain wm (=ws) gets computed 1620. c directly. 1621. INCLUDE "KPP_1D.COM" 1622. INTENT (IN) sigma,hbl,ustar,bfsfc 1623. INTENT (OUT) wm,ws 1624. 1625. c input 1626. real*8 sigma ! normalized depth (d/hbl) 1627. real*8 hbl ! boundary layer depth (m) 1628. real*8 ustar ! surface friction velocity (m/s) 1629. real*8 bfsfc ! total surface buoyancy flux (m^2/s^3) 1630. c output 1631. real*8 wm,ws ! turbulent velocity scales at sigma 1632. c local 1633. real*8 zehat ! = zeta * ustar**3 1634. 1635. c use lookup table for zehat < zmax only; otherwise use 1636. c stable formulae 1637. 1638. zehat = vonk * sigma * hbl * bfsfc 1639. 1640. if (zehat.le.zmax) then 1641. zdiff = zehat-zmin 1642. iz = int( zdiff * rdeltaz ) 1643. iz = min( iz , nni ) 1644. iz = max( iz , 0 ) 1645. izp1=iz+1 1646. 1647. udiff = ustar-umin 1648. ju = int( udiff * rdeltau) 1649. ju = min( ju , nnj ) 1650. ju = max( ju , 0 ) 1651. jup1=ju+1 1652. 1653. zfrac = zdiff*rdeltaz - float(iz) 1654. ufrac = udiff*rdeltau - float(ju) 1655. 1656. fzfrac= 1.-zfrac 1657. wam = (fzfrac) * wmt(iz,jup1) + zfrac*wmt(izp1,jup1) 1658. wbm = (fzfrac) * wmt(iz,ju ) + zfrac*wmt(izp1,ju ) 1659. wm = (1.-ufrac)* wbm + ufrac*wam 1660. 1661. was = (fzfrac) * wst(iz,jup1) + zfrac*wst(izp1,jup1) 1662. wbs = (fzfrac) * wst(iz,ju ) + zfrac*wst(izp1,ju ) 1663. ws = (1.-ufrac)* wbs + ufrac*was 1664. else 1665. u3 = ustar*ustar*ustar 1666. wm = vonk * ustar * u3 / ( u3 + conc1*zehat ) 1667. ws = wm 1668. endif 1669. 1670. return 1671. end 1672. 1673. subroutine z121 (v,kmtj,km) 1674. c Apply 121 smoothing in k to 2-d array V(k=1,km) 1675. c top (0) value is used as a dummy 1676. c bottom (km+1) value is set to input value from above. 1677. 1678. PARAMETER (p5=5d-1, p25=2.5d-1) 1679. INTENT (IN) kmtj,km 1680. INTENT (INOUT) v 1681. c input 1682. real*8 V(0:km+1) ! 2-D array to be smoothed 1683. 1684. V(0) = p25 * V(1) 1685. V(kmtj+1) = V(kmtj) 1686. 1687. do k=1,kmtj 1688. tmp = V(k) 1689. V(k) = V(0) + p5 * V(k) + p25 * V(k+1) 1690. V(0) = p25 * tmp 1691. end do 1692. return 1693. end 1694. 1695. subroutine ddmix (km, alphaDT, betaDS, visc, difs, dift, kmtj) 1696. 1697. c Rrho dependent interior flux parameterization. 1698. c Add double-diffusion diffusivities to Ri-mix values at blending 1699. c interface and below. 1700. 1701. INCLUDE "KPP_1D.COM" 1702. INTENT (IN) km,alphaDT,betaDS,kmtj 1703. INTENT (INOUT) visc, difs, dift 1704. c input 1705. real*8 alphaDT(km) ! alpha * DT across interfaces 1706. real*8 betaDS(km) ! beta * DS across interfaces 1707. c output 1708. real*8 visc(0:km+1) ! interior viscosity (m^2/s) 1709. real*8 dift(0:km+1) ! interior thermal diffusivity (m^2/s) 1710. real*8 difs(0:km+1) ! interior scalar diffusivity (m^2/s) 1711. c local 1712. real*8 diffdd ! double diffusion diffusivity scale 1713. real*8 prandtl ! prandtl number 1714. 1715. do 100 ki= 1, kmtj 1716. 1717. c salt fingering case 1718. 1719. if((alphaDT(ki).gt.betaDS(ki)).and. 1720. $ (betaDS (ki).gt.0. )) then 1721. 1722. Rrho = MIN(alphaDT(ki) / betaDS(ki) , Rrho0) 1723. c diffdd = dsfmax*(1.0-((Rrho-1)/(Rrho0-1))**2)**pexp2 1724. diffdd = 1.0-((Rrho-1)/(Rrho0-1))**2 1725. diffdd = dsfmax*diffdd*diffdd*diffdd 1726. dift(ki) = dift(ki) + 0.7*diffdd 1727. difs(ki) = difs(ki) + diffdd 1728. 1729. c diffusive convection 1730. 1731. else if ((alphaDT(ki).lt.0.0).and.(betaDS(ki).lt.0.0) 1732. $ .and.(alphaDT(ki).lt.betaDS(ki)) ) then 1733. 1734. Rrho = alphaDT(ki) / betaDS(ki) 1735. diffdd=1.5d-6*9.0*0.101d0*exp(4.6d0*exp(-0.54d0*(1/Rrho-1))) 1736. prandtl = 0.15*Rrho 1737. if (Rrho.gt.0.5d0) prandtl = (1.85d0-0.85d0/Rrho)*Rrho 1738. dift(ki) = dift(ki) + diffdd 1739. difs(ki) = difs(ki) + prandtl*diffdd 1740. endif 1741. 100 continue 1742. return 1743. end 1744. 1745. subroutine kmixinit 1746. c initialize some constants for kmix subroutines, and initialize 1747. c for kmix subroutine "wscale" the 2D-lookup table for wm and ws 1748. c as functions of ustar and zetahat (=vonk*sigma*hbl*bfsfc). 1749. 1750. INCLUDE "KPP_1D.COM" 1751. c local 1752. real*8 zehat ! = zeta * ustar**3 1753. real*8 zeta ! = stability parameter d/L 1754. PARAMETER (RFRAC=.62d0,ZETA1=1.5d0, ZETA2=2d1, 1755. * ZE1=12d0,ZERAT=1.5d0) 1756. EF(Z) = RFRAC*EXP(-Z/ZETA1) + (1d0-RFRAC)*EXP(-Z/ZETA2) 1757. 1758. c define some non-dimensional constants and 1759. c the vertical mixing coefficients in m-k-s units 1760. 1761. c epsl = epsln 1762. Vtc = concv * sqrt(0.2/concs/epsilon) / vonk**2 / Ricr 1763. cg = cstar * vonk * (concs * vonk * epsilon)**(1./3.) 1764. 1765. c difm0 = vvcric * r10000 1766. c difs0 = vdcric * r10000 1767. c difmiw = fkpm * r10000 1768. c difsiw = fkph * r10000 1769. c difmcon = vvclim * r10000 1770. c difscon = vdclim * r10000 1771. 1772. c construct the wm and ws lookup tables 1773. 1774. c deltaz = (zmax-zmin)/(nni+1) 1775. c deltau = (umax-umin)/(nnj+1) 1776. c rdeltaz = 1./deltaz 1777. c rdeltau = 1./deltau 1778. 1779. do 100 i=0,nni+1 1780. zehat = deltaz*(i) + zmin 1781. do 90 j=0,nnj+1 1782. usta = deltau*(j) + umin 1783. zeta = zehat/(usta**3+epsl) 1784. 1785. if(zehat.ge.0.) then 1786. wmt(i,j) = vonk*usta/(1.+conc1*zeta) 1787. wst(i,j) = wmt(i,j) 1788. else 1789. if(zeta.gt.zetam) then 1790. wmt(i,j) = vonk* usta * sqrt(sqrt(1.-conc2*zeta)) 1791. else 1792. wmt(i,j) = vonk* (conam*usta**3-concm*zehat)**(1./3.) 1793. endif 1794. if(zeta.gt.zetas) then 1795. wst(i,j) = vonk* usta * sqrt(1.-conc3*zeta) 1796. else 1797. wst(i,j) = vonk* (conas*usta**3-concs*zehat)**(1./3.) 1798. endif 1799. endif 1800. 90 continue 1801. 100 continue 1802. 1803. c set up values used in swfrac 1804. do i=0,LMO 1805. ZE(i)=ZE1*(ZERAT**i-1d0)/(ZERAT-1d0) 1806. end do 1807. do i=1,ksrpd 1808. E0TO(i) = EF(ZE(i-1)) 1809. EZTOI(i)= E0TO(i) /(ZE(i)-ZE(i-1)) 1810. end do 1811. do i=1,ksrpd-1 1812. EZTO(i) = (E0TO(i)-E0TO(i+1))/(ZE(i)-ZE(i-1)) 1813. end do 1814. 1815. return 1816. end 1817. 1818. c block data kmixbdta 1819. c 1820. c include "KPP_1D.COM" 1821. c 1822. c data statements needed for "kmix" vertical mixing subroutines 1823. c in "kmixs.F" 1824. c 1825. c parameters for several subroutines 1826. c data epsilon / 0.1d0 / 1827. c data vonk / 0.4d0 / 1828. c data conc1 / 5d0 / 1829. c data conam,concm,conc2,zetam / 1.257d0,8.380d0,16d0,- 0.2d0 / 1830. c data conas,concs,conc3,zetas / -28.86d0,98.96d0,16d0, - 1d0 / 1831. c 1832. c parameters for subroutine "bldepth" 1833. c Ricr = critical bulk Richardson Number = 0.3 1834. c hbf = fraction of bounadry layer depth to 1835. c which absorbed solar radiation 1836. c contributes to surface buoyancy forcing = 1.0 1837. c 1838. c data Ricr / 0.3d0 / 1839. c data cekman / 0.7d0 / 1840. c data cmonob / 1d0 / 1841. c data concv / 1.8d0 / 1842. c data hbf / 1d0 / 1843. c 1844. c limits for lookup table of wm and ws in subroutine "kmixinit" 1845. c (zmin and zmax are in m3/s3, umin and umax are in m/s) 1846. c 1847. c data zmin,zmax /-4d-7, 0. / 1848. c data umin,umax / 0. , 4d-2/ 1849. c 1850. c parameters for subroutine "ri_iwmix" 1851. c to compute vertical mixing coefficients below boundary layer: 1852. c 1853. c Riinfty = local Richardson Number limit 1854. c for shear instability 1855. c rRiinfty = 1/Riinfty 1856. c BVSQcon = value of N^2 where the convection coefficients 1857. c first become their maxima (small negative number) 1858. c rBVSQcon = 1/BVSQcon 1859. c 1860. c data Riinfty / 0.7d0 / 1861. c data rRiinfty / 1.43d0 / 1862. c data BVSQcon / -1d-7 / 1863. c data rBVSQcon / -1d7 / 1864. c data num_v_smooth_Ri / 1 / 1865. c 1866. c Note: Data statements for vvcric, vdcric, fkpm, fkph, vvclim, 1867. c ^^^^^ and vdclim, formerly in this file, from blkdta.F 1868. c * The vertical mixing coefficients are defined 1869. c in (m2/s) units in "kmixinit". they are initialized 1870. c in "kmixbdta" and read via namelist "eddy" in (cm2/s) 1871. c units using their c-g-s names) 1872. c (m2/s) (cm2/s) 1873. c difm0 = viscosity max due to shear instability = vvcric 1874. c difs0 = diffusivity .. = vdcric 1875. c difmiw = viscosity background due to internal waves = fkpm 1876. c difsiw = diffusivity .. = fkph 1877. c difmcon = viscosity due to convective instability = vvclim 1878. c difscon = diffusivity .. = vdclim 1879. c 1880. c data vvcric / 50d0 / 1881. c data vdcric / 50d0 / 1882. c data fkpm / 10d0 / 1883. c data fkph / 0.3d0 / 1884. c data vvclim / 1000d0 / 1885. c data vdclim / 1000d0 / 1886. c 1887. c parameters for subroutine "ddmix" 1888. c to compute additional diffusivity due to double diffusion: 1889. c Rrho0 = limit for double diffusive density ratio 1890. c dsfmax = maximum diffusivity in case of salt fingering (m2/s) 1891. c 1892. c data Rrho0 / 1.9d0 / 1893. c data dsfmax / 0.001d0 / 1894. c 1895. c parameters for subroutine "blmix" 1896. c to compute mixing within boundary layer: 1897. c cstar = proportionality coefficient 1898. c for nonlocal transport 1899. c 1900. c data cstar / 10d0 / 1901. c 1902. c end 1903. 1904. subroutine swfrac(z,k,kmtj,kmax,bfsfc) 1905. 1906. C**** Calculate fraction of solar energy penetrating to depth z 1907. C**** using linear interpolation for depths between grid levels 1908. C**** There is a slight error since 'z' is scaled by free surface height 1909. C**** and ZE is not. 1910. C**** 1911. C**** k is grid box containing point 1912. C**** kmtj is number of grid points in column 1913. C**** kmax is the number of grid points in column that recieve 1914. C**** solar radiation 1915. C**** 1916. INCLUDE "KPP_1D.COM" 1917. INTENT (IN) z,k,kmtj,kmax 1918. INTENT (OUT) bfsfc 1919. real*8 z,bfsfc 1920. 1921. C**** Get actual k (since k is on staggered grid) 1922. kt = k + NINT(0.5d0 + SIGN(0.5d0,-(ZE(k)+z))) 1923. 1924. C**** calculate fraction 1925. if (kt.gt.kmax) then 1926. bfsfc = 0. 1927. else 1928. if (kt.eq.kmax) then 1929. bfsfc = E0TO(kt) + (z+ZE(kt-1))*EZTOI(kt) 1930. else 1931. bfsfc = E0TO(kt) + (z+ZE(kt-1))*EZTO(kt) 1932. end if 1933. end if 1934. 1935. return 1936. end 4000. 4001. SUBROUTINE OCONV 4002. C**** 4003. C**** OCONV uses vertical diffusion coefficients from KVMIX 4004. C**** using KPP boundary layer code. If LCONV = .TRUE. also 4005. C**** OCONV allows oceanic convective plumes to sink from below boundary 4006. C**** layer to reach another level at which the plume's 4007. C**** potential density is equal to that of the surrounding water. 4008. C**** Both processes are checked and applied on every 4009. C**** horizontal quarter box. 4010. C**** 4011. INCLUDE 'C070.COM' 4012. LOGICAL*4 QPOLE, LCONV,LDD 4013. COMMON /WORK01/UT(IM,JM,LMO),VT(IM,JM,LMO),UL(LMO,IM+2),MML(LMO), 4014. * G0ML(LMO,2,2),S0ML(LMO,2,2),GZML(LMO,2,2),SZML(LMO,2,2), 4015. * BYMML(LMO),DTBYDZ(LMO),BYDZ2(LMO), 4016. * RAVM(IM+1),RAMV(IM+1),ID(IM+1),LMUV(IM+1), 4017. * UL0(LMO,IM+2),G0ML0(LMO,2,2),S0ML0(LMO,2,2),MML0(LMO), 4018. * BYMML0(LMO),MMLT(LMO),BYMMLT(LMO) 4019. COMMON /DIFVCB/ G0M1(IM,JM,3),MO1(IM,JM),GXM1(IM,JM),GYM1(IM,JM) 4020. * ,SXM1(IM,JM),SYM1(IM,JM),UO1(IM,JM),VO1(IM,JM),AKVM(0:LMO+1) 4021. * ,AKVG(0:LMO+1),AKVS(0:LMO+1),GHATM(LMO),GHATG(LMO),GHATS(LMO) 4022. COMMON /FLUXCB/ SROX(IM,JM),DMUA(IM,JM,2),DMVA(IM,JM,2), 4023. * W0(IM,JM,4),E0(IM,JM,4),E1(IM,JM,4),DMUI(IM,JM),DMVI(IM,JM) 4024. C**** CONV parameters: BETA controls degree of convection (default 0.5). 4025. PARAMETER (BETA=5d-1, BYBETA=1d0/BETA) 4026. C**** KPP variables 4027. PARAMETER(epsln=1d-20,OMEGA=TWOPI*366d0/(365d0*SDAY)) 4028. REAL*8 zgrid(0:LMO+1),hwide(0:LMO+1),Shsq(LMO),dVsq(LMO) 4029. * ,talpha(LMO),sbeta(LMO),dbloc(LMO),dbsfc(LMO),Ritop(LMO), 4030. * alphaDT(LMO),betaDS(LMO),ghat(LMO),byhwide(0:LMO+1) 4031. REAL*8 G(LMO),S(LMO),T(LMO),BYRHO(LMO),RHO(LMO),P(LMO) 4032. DATA LCONV/.FALSE./,LDD/.FALSE./ 4033. DATA IFIRST/1/,ILAST/0/,JLAST/0/ 4034. C**** initiallise kpp routines 4035. IF (IFIRST.eq.1) THEN 4036. call kmixinit 4037. BYDTS = 1d0/DTS 4038. IFIRST = 0. 4039. END IF 4040. C**** Load UO,VO into UT,VT. UO,VO will be updated, while UT,VT 4041. C**** will be fixed during convection. 4042. DO I=1,IM*JM*LMO 4043. UT(I,1,1) = UO(I,1,1) 4044. VT(I,1,1) = VO(I,1,1) 4045. END DO 4046. C**** 4047. C**** Outside loop over J 4048. C**** 4049. DO 790 J=2,JM 4050. C**** coriolis parameter, defined at tracer point 4051. Coriol = 2d0*OMEGA*SINP(J) 4052. 4053. QPOLE = J.EQ.JM 4054. IF (QPOLE) THEN 4055. C**** 4056. C**** Define parameters and currents at the North Pole 4057. C**** 4058. LMIJ=LMM(1,JM) 4059. IQ=1 4060. JQ=1 4061. KMUV=IM+1 4062. DO 110 I=1,IM 4063. ID(I) = I + IM*(JM-2) + IM*JM*LMO 4064. LMUV(I) = LMV(I,JM-1) 4065. RAVM(I) = 1d0/IM 4066. RAMV(I) = RAMVS(JM) 4067. UL0(1,I) = VO(I,JM-1,1) 4068. UL(1,I) = VO1(I,JM-1) 4069. DO 110 L=2,LMIJ 4070. UL0(L,I) = VO(I,JM-1,L) 4071. UL(L,I) = UL0(L,I) 4072. 110 CONTINUE 4073. ID(IM+1) = 1 + IM*(JM-1) 4074. LMUV(IM+1) = LMU(1,JM) 4075. RAVM(IM+1) = 1d0 4076. RAMV(IM+1) = 1d0 4077. UL0(1,IM+1) = UO(1,JM,1) 4078. MML0(1) = MO(1,JM,1)*DXYP(JM) 4079. BYMML0(1) = 1d0/MML0(1) 4080. DTBYDZ(1) = DTS/MO(1,JM,1) 4081. G0ML0(1,1,1)= G0M(1,JM,1) 4082. S0ML0(1,1,1)= S0M(1,JM,1) 4083. GZML(1,1,1) = GZM(1,JM,1) 4084. SZML(1,1,1) = SZM(1,JM,1) 4085. UL(1,IM+1) = UO1(1,JM) 4086. MMLT(1) = MO1(1,JM)*DXYP(JM) 4087. BYMMLT(1) = 1d0/MMLT(1) 4088. S0ML(1,1,1) = S0M(1,JM,1) 4089. DO L=2,LMIJ 4090. UL0(L,IM+1) = UO(1,JM,L) 4091. MML0(L) = MO(1,JM,L)*DXYP(JM) 4092. BYMML0(L) = 1d0/MML0(L) 4093. DTBYDZ(L) = DTS/MO(1,JM,L) 4094. BYDZ2(L-1) = 2d0/(MO(1,JM,L-1)+MO(1,JM,L)) 4095. G0ML0(L,1,1)= G0M(1,JM,L) 4096. S0ML0(L,1,1)= S0M(1,JM,L) 4097. GZML(L,1,1) = GZM(1,JM,L) 4098. SZML(L,1,1) = SZM(1,JM,L) 4099. UL(L,IM+1) = UL0(L,IM+1) 4100. MMLT(L) = MML0(L) 4101. BYMMLT(L) = BYMML0(L) 4102. S0ML(L,1,1) = S0ML0(L,1,1) 4103. END DO 4104. DO L=1,MIN(3,LMIJ) 4105. G0ML(L,1,1) = G0M1(1,JM,L) 4106. END DO 4107. DO L=3,LMIJ 4108. G0ML(L,1,1) = G0ML0(L,1,1) 4109. END DO 4110. C**** Calculate whole box turbulent stresses sqrt(|tau_0|/rho_0) 4111. C**** except for RHO factor and sqrt: U2rho = Ustar**2 * rho 4112. C**** DM[UV]A are defined on t-grid. DM[UV]I defined on u,v grid 4113. C**** Note that rotational ice stresses are not calculated at the pole 4114. UISTR = 0 4115. VISTR = 0 4116. DO I=1,IM 4117. UISTR = UISTR + DMVI(I,JM-1)*COSI(I) 4118. VISTR = VISTR - DMVI(I,JM-1)*SINI(I) 4119. END DO 4120. UISTR = UISTR/IM 4121. VISTR = VISTR/IM 4122. U2rho = SQRT( 4123. * (DMUA(1,JM,1)*(1d0-RSI(1,JM)) + UISTR*RSI(1,JM))**2+ 4124. * (DMVA(1,JM,1)*(1d0-RSI(1,JM)) + VISTR*RSI(1,JM))**2)*BYDTS 4125. C**** Calculate surface freshwater and heat fluxes 4126. DELTAFW = (MO(1,JM,1) - MO1(1,JM))*BYDTS 4127. DELTAE = (G0ML0(1,1,1) - G0ML(1,1,1))*BYDXYP(JM)*BYDTS 4128. DELTASR = SROX(1,JM)*(1d0-RSI(1,JM))*BYDTS ! W/m^2 4129. I=1 4130. GOTO 500 4131. END IF 4132. C**** 4133. C**** Define parameters and currents away from the poles 4134. C**** 4135. IM1=IM 4136. I=1 4137. 210 IF(LMM(I,J).LE.1) GO TO 730 4138. LMIJ=LMM(I,J) 4139. KMUV=4 4140. ID(1) = IM1 + IM*(J-1) 4141. ID(2) = I + IM*(J-1) 4142. ID(3) = I + IM*(J-2) + IM*JM*LMO 4143. ID(4) = I + IM*(J-1) + IM*JM*LMO 4144. LMUV(1) = LMU(IM1,J) 4145. LMUV(2) = LMU(I ,J) 4146. LMUV(3) = LMV(I,J-1) 4147. LMUV(4) = LMV(I,J ) 4148. MML0(1) = 2.5d-1*MO(I,J,1)*DXYP(J) 4149. BYMML0(1) = 1d0/MML0(1) 4150. DTBYDZ(1) = DTS/MO(I,J,1) 4151. MMLT(1) = 2.5d-1*MO1(I,J)*DXYP(J) 4152. BYMMLT(1) = 1d0/MMLT(1) 4153. DO 220 L=2,LMIJ 4154. MML0(L) = 2.5d-1*MO(I,J,L)*DXYP(J) 4155. BYMML0(L) = 1d0/MML0(L) 4156. DTBYDZ(L) = DTS/MO(I,J,L) 4157. BYDZ2(L-1)= 2d0/(MO(I,J,L)+MO(I,J,L-1)) 4158. MMLT(L) = MML0(L) 4159. BYMMLT(L) = BYMML0(L) 4160. 220 CONTINUE 4161. 4162. C**** Calculate whole box turbulent stresses sqrt(|tau_0|/rho_0) 4163. C**** except for RHO factor and sqrt: U2rho = Ustar**2 * rho 4164. C**** DM[UV]A are defined on t-grid. DM[UV]I defined on u,v grid 4165. UISTR = 0 4166. VISTR = 0 4167. IF (RSI(I,J).gt.0) THEN 4168. IF (LMU(I,J).gt.0) UISTR = UISTR + DMUI(I,J) 4169. IF (LMU(IM1,J).gt.0) UISTR = UISTR + DMUI(IM1,J) 4170. ANSTR=1.-SIGN(0.25,LMU(I,J)-0.5)-SIGN(0.25,LMU(IM1,J)-0.5) 4171. UISTR = UISTR*ANSTR*RSI(I,J) 4172. IF (LMV(I,J).gt.0) VISTR = VISTR + DMVI(I,J) 4173. IF (LMV(I,J-1).gt.0) VISTR = VISTR + DMVI(I,J-1) 4174. ANSTR=1.-SIGN(0.25,LMV(I,J)-0.5)-SIGN(0.25,LMV(I,J-1)-0.5) 4175. VISTR = VISTR*ANSTR*RSI(I,J) 4176. END IF 4177. U2rho = SQRT((DMUA(I,J,1)*(1d0-RSI(I,J)) + UISTR)**2 + 4178. * (DMVA(I,J,1)*(1d0-RSI(I,J)) + VISTR)**2)*BYDTS 4179. C**** Calculate surface freshwater flux and Solar forcing 4180. DELTAFW = (MO(I,J,1) - MO1(I,J))*BYDTS ! kg/m^2 s 4181. DELTASR = SROX(I,J)*(1d0-RSI(I,J))*BYDTS ! W/m^2 4182. 4183. C**** Loop over quarter boxes 4184. JQ=1 4185. 230 RJ = BETA*(2d0*JQ-3d0) 4186. RAVM(3) = (1.25d0-5d-1*JQ) 4187. RAMV(3) = (1.25d0-5d-1*JQ)*RAMVS(J)*5d-1 4188. RAVM(4) = (5d-1*JQ-2.5d-1) 4189. RAMV(4) = (5d-1*JQ-2.5d-1)*RAMVN(J)*5d-1 4190. IQ=1 4191. 240 RI = BETA*(2d0*IQ-3d0) 4192. RAVM(1) = (1.25d0 - 5d-1*IQ) 4193. RAMV(1) = (1.25d0 - 5d-1*IQ)*2.5d-1 4194. RAVM(2) = (5d-1*IQ-2.5d-1) 4195. RAMV(2) = (5d-1*IQ-2.5d-1)*2.5d-1 4196. UL0(1,1) = UT(IM1,J,1) 4197. UL0(1,2) = UT(I ,J,1) 4198. UL0(1,3) = VT(I,J-1,1) 4199. UL0(1,4) = VT(I,J ,1) 4200. G0ML0(1,IQ,JQ)=2.5d-1*(G0M(I,J,1) + RI*GXM(I,J,1) + RJ*GYM(I,J,1)) 4201. S0ML0(1,IQ,JQ)=2.5d-1*(S0M(I,J,1) + RI*SXM(I,J,1) + RJ*SYM(I,J,1)) 4202. GZML(1,IQ,JQ)=2.5d-1* GZM(I,J,1) 4203. SZML(1,IQ,JQ)=2.5d-1* SZM(I,J,1) 4204. UL(1,1) = UO1(IM1,J) 4205. UL(1,2) = UO1(I ,J) 4206. UL(1,3) = VO1(I,J-1) 4207. UL(1,4) = VO1(I,J ) 4208. S0ML(1,IQ,JQ)=2.5d-1*(S0M(I,J,1) + RI*SXM1(I,J) + RJ*SYM1(I,J)) 4209. DO 250 L=2,LMIJ 4210. UL0(L,1) = UT(IM1,J,L) 4211. UL0(L,2) = UT(I ,J,L) 4212. UL0(L,3) = VT(I,J-1,L) 4213. UL0(L,4) = VT(I,J ,L) 4214. G0ML0(L,IQ,JQ)=2.5d-1*(G0M(I,J,L) + RI*GXM(I,J,L) + RJ*GYM(I,J,L)) 4215. S0ML0(L,IQ,JQ)=2.5d-1*(S0M(I,J,L) + RI*SXM(I,J,L) + RJ*SYM(I,J,L)) 4216. GZML(L,IQ,JQ)=2.5d-1* GZM(I,J,L) 4217. SZML(L,IQ,JQ)=2.5d-1* SZM(I,J,L) 4218. UL(L,1) = UL0(L,1) 4219. UL(L,2) = UL0(L,2) 4220. UL(L,3) = UL0(L,3) 4221. UL(L,4) = UL0(L,4) 4222. S0ML(L,IQ,JQ) = S0ML0(L,IQ,JQ) 4223. 250 CONTINUE 4224. G0ML(1,IQ,JQ)=2.5d-1*(G0M1(I,J,1)+ RI*GXM1(I,J) + RJ*GYM1(I,J)) 4225. DO L=2,MIN(3,LMIJ) 4226. G0ML(L,IQ,JQ)=2.5d-1*(G0M1(I,J,L)+RI*GXM(I,J,L)+RJ*GYM(I,J,L)) 4227. END DO 4228. DO L=3,LMIJ 4229. G0ML(L,IQ,JQ) = G0ML0(L,IQ,JQ) 4230. END DO 4231. C**** Calculate surface heat fluxe 4232. DELTAE=4d0*(G0ML0(1,IQ,JQ)-G0ML(1,IQ,JQ))*BYDXYP(J)*BYDTS !W/m^2 4233. C**** 4234. C**** Vertical mixing dependent on KPP boundary layer scheme 4235. C**** 4236. 500 IF(LMIJ.LE.1) GO TO 725 4237. C**** Do not recalculate quantities good for all quarter boxes 4238. IF (I.ne.ILAST .or. J.ne.JLAST) THEN 4239. 4240. C**** Z0 = BLDATA(8)/GRAV-ZSOLID Ocean depth (m) 4241. C**** scale depths using Z0/ZE(LMIJ) 4242. 4243. ZSCALE = (BLDATA(I,J,8) - GRAV*ZSOLID(I,J))/(ZE(LMIJ)*GRAV) 4244. zgrid(0) = epsln 4245. hwide(0) = epsln 4246. byhwide(0) = 0. 4247. P(1) = 5d-1*MO(I,J,1)*GRAV 4248. DO L=1,LMO-1 4249. zgrid(L) = -0.5*ZSCALE*(ZE(L-1) + ZE(L)) ! tracer level 4250. hwide(L) = zgrid(L-1) - zgrid(L) ! tracer level 4251. byhwide(L) = 1d0/hwide(L) 4252. P(L+1) = P(L) + 5d-1*GRAV*(MO(I,J,L)+MO(I,J,L+1)) 4253. END DO 4254. zgrid(LMO) = -0.5*ZSCALE*(ZE(LMO-1) + ZE(LMO)) ! tracer level 4255. hwide(LMO) = zgrid(LMO-1) - zgrid(LMO) ! tracer level 4256. byhwide(LMO) = 1d0/hwide(LMO) 4257. zgrid(LMO+1) = -ZE(LMO)*ZSCALE 4258. hwide(LMO+1) = epsln 4259. byhwide(LMO+1) = 0. 4260. 4261. END IF 4262. ILAST = I 4263. JLAST = J 4264. 4265. C**** If swfrac is not used add solar from next two levels to Bo 4266. C**** E0TO1 = 0.20876 (see OSOURC) CAUTION: Resolution dependent! 4267. C DELTAE = DELTAE + 0.20876d0*DELTASR 4268. C**** else remove solar from Bo 4269. DELTAE = DELTAE - 0.79124d0*DELTASR 4270. 4271. C**** Start iteration for diffusivities 4272. MML = MMLT ! start with pre-source mass 4273. BYMML = BYMMLT 4274. ITER = 0 4275. HBL = 0 4276. 510 ITER =ITER + 1 4277. HBLP = HBL 4278. IF (ITER.eq.2) THEN 4279. MML = MML0 ! continue with post-source mass 4280. BYMML = BYMML0 4281. END IF 4282. C**** Initiallise fields 4283. Shsq=0. 4284. dVsq=0. 4285. dbloc=0. 4286. dbsfc=0. 4287. Ritop=0. 4288. 4289. C**** velocity shears 4290. DO L=1,LMIJ-1 4291. DO K=1,KMUV 4292. Shsq(L) = Shsq(L) + RAVM(K)*(UL(L,K)-UL(L+1,K))**2 !interface 4293. dVsq(L) = dVsq(L) + RAVM(K)*(UL(1,K)-UL(L ,K))**2 !tracer pts 4294. END DO 4295. END DO 4296. DO K=1,KMUV 4297. dVsq(LMIJ)=dVsq(LMIJ)+RAVM(K)*(UL(1,K)-UL(LMIJ,K))**2 4298. END DO 4299. 4300. C**** density related quantities 4301. C**** 4302. C**** density of surface layer (kg/m3) 4303. C**** rho = rho{G(1),S(1),P(1)} 4304. C**** local buoyancy gradient at km interfaces: (m/s2) 4305. C**** dbloc = g/rho{k+1,k+1} * [ rho{k,k+1}-rho{k+1,k+1} ] 4306. C**** buoyancy difference with respect to "zref", the surface (m/s2) 4307. C**** dbsfc = g/rho{k,k} * [ rho{k,k}-rho{1,k} ] (CORRECT) 4308. C**** thermal expansion coefficient without 1/rho factor (kg/m3/C) 4309. C**** talpha= d(rho{k,k})/d(T(k)) 4310. C**** salt expansion coefficient without 1/rho factor (kg/m3/PSU) 4311. C**** sbeta = d(rho{k,k})/d(S(k)) 4312. 4313. G(1) = G0ML(1,IQ,JQ)*BYMML(1) 4314. S(1) = S0ML(1,IQ,JQ)*BYMML(1) 4315. BYRHO(1) = VOLGSP(G(1),S(1),P(1)) 4316. RHO(1) = 1d0/BYRHO(1) 4317. talpha(1) = ALPHAGSP(G(1),S(1),P(1)) ! >0 4318. sbeta(1) = BETAGSP(G(1),S(1),P(1)) 4319. DO L=2,LMIJ 4320. G(L) = G0ML(L,IQ,JQ)*BYMML(L) 4321. S(L) = S0ML(L,IQ,JQ)*BYMML(L) 4322. BYRHO(L) = VOLGSP(G(L),S(L),P(L)) 4323. RHO(L) = 1d0/BYRHO(L) 4324. RHOM = 1d0/VOLGSP(G(L-1),S(L-1),P(L)) ! L-1 density wrt P_L 4325. RHO1 = 1d0/VOLGSP(G(1),S(1),P(L)) ! surf density wrt P_L 4326. dbsfc(L) = GRAV*(1d0-RHO1*BYRHO(L)) 4327. dbloc(L-1) = GRAV*(1d0-RHOM*BYRHO(L)) 4328. C**** numerator of bulk richardson number on grid levels 4329. Ritop(L) = (zgrid(1)-zgrid(L)) * dbsfc(L) 4330. END DO 4331. 4332. C**** surface wind turbulent friction speed (m/s) = sqrt( |tau_0|/rho ) 4333. Ustar = SQRT(U2rho*BYRHO(1)) 4334. 4335. C**** surface buoyancy forcing turbulent and solar (m^2/s^3) 4336. C**** Bo = -g*(talpha * wtsfc - sbeta * wssfc)/rho 4337. C**** Bosol = -g*(talpha * wtsol )/rho 4338. C**** Bo includes all buoyancy and heat forcing 4339. 4340. BYSHC = 1d0/SHCGS(G1,S(1)) 4341. 4342. Bo = GRAV*(-talpha(1)*DELTAE *BYSHC + sbeta(1)*DELTAFW*S(1)*1d3) 4343. * *BYRHO(1)**2 4344. Bosol = - GRAV* talpha(1)*DELTASR*BYSHC 4345. * *BYRHO(1)**2 4346. 4347. C**** Double diffusive option (ddmix) needs alphaDT,betaDS 4348. C**** alphaDT = mean -talpha * delta(temp.) at interfaces (kg/m3) 4349. C**** betaDS = mean sbeta * delta(salt) at interfaces (kg/m3) 4350. 4351. if (LDD) then 4352. T(1) = TEMGSP(G(1),S(1),P(1)) 4353. do L=2,LMIJ 4354. T(L) = TEMGSP(G(L),S(L),P(L)) 4355. talpha(L)= ALPHAGSP(G(L),S(L),P(L)) ! >0 4356. sbeta(L) = BETAGSP(G(L),S(L),P(L)) 4357. alphaDT(L-1) = - 5d-1 * (talpha(L-1) + talpha(L)) 4358. $ * (T(L-1) - T(L)) 4359. betaDS (L-1) = 5d-1 * (sbeta (L-1) + sbeta(L)) 4360. $ * (S(L-1) - S(L))*1d3 4361. end do 4362. end if 4363. 4364. CALL KPPMIX(LDD,LMO,zgrid,hwide,LMIJ,Shsq,dVsq,Ustar,Bo,Bosol 4365. * ,alphaDT,betaDS,dbloc,Ritop,Coriol,byhwide, 4366. * AKVM,AKVS,AKVG,GHAT,HBL,KBL) 4367. 4368. C**** Calculate non-local transport terms for each scalar 4369. C**** ghat[sg] = kv * ghat * (J,kg) 4370. C**** Correct units for diffusivities (m^2/s) => (kg^2/m^4 s) 4371. C**** ghat (s/m^2) => (s m^4/kg^2) 4372. DO L=1,LMIJ-1 4373. R = 5d-1*(RHO(L)+RHO(L+1)) 4374. R2 = R**2 4375. GHATM(L) = 0. ! no non-local momentum transport 4376. C**** GHAT terms must be zero for consistency with OSOURC 4377. GHATG(L) = 0. !AKVG(L)*GHAT(L)*DELTAE*DXYP(J) 4378. GHATS(L) = 0. !AKVS(L)*GHAT(L)*DELTAFW*S(1)*DXYP(J) 4379. AKVM(L) = AKVM(L)*R2 4380. AKVG(L) = AKVG(L)*R2 4381. AKVS(L) = AKVS(L)*R2 4382. END DO 4383. 4384. C**** For each field (U,G,S + TRACERS) call VDIFF 4385. C**** Momentum 4386. DO K=1,KMUV 4387. CALL VDIFF(UL(1,K),AKVM(1),GHATM,DTBYDZ,BYDZ2,DTS 4388. * ,LMIJ,UL0(1,K)) 4389. END DO 4390. C**** Enthalpy 4391. CALL VDIFFS(G0ML(1,IQ,JQ),AKVG(1),GHATG,DTBYDZ,BYDZ2,DTS 4392. * ,LMIJ,G0ML0(1,IQ,JQ)) 4393. C**** Salinity 4394. CALL VDIFFS(S0ML(1,IQ,JQ),AKVS(1),GHATS,DTBYDZ,BYDZ2,DTS 4395. * ,LMIJ,S0ML0(1,IQ,JQ)) 4396. IF ((ITER.eq.1 .or. ABS(HBLP-HBL).gt.(ZE(KBL)-ZE(KBL-1))*0.25) 4397. * .and. ITER.lt.4) GO TO 510 4398. C IF (ITER.eq.4) PRINT*,'HBL ',I,J,KBL,HBLP,HBL 4399. c IF (ITER.eq.1) GO TO 510 4400. C**** Gradients of scalars? 4401. C**** Choice: 1) Apply interpolated KV to linear profile 4402. C**** 2) Double up resolution and rederive gradients 4403. C**** 3) Rederive gradients 4404. C**** Choice 1B: Same as above but implicitly 4405. GZML(1,IQ,JQ) = 0. 4406. SZML(1,IQ,JQ) = 0. 4407. DO L=2,LMIJ-1 4408. DTBYDZ2 = 2d0*DTBYDZ(L)**2*BYDTS 4409. GZML(L,IQ,JQ)=GZML(L,IQ,JQ)/(1d0+DTBYDZ2*(AKVG(L-1)+AKVG(L))) 4410. SZML(L,IQ,JQ)=SZML(L,IQ,JQ)/(1d0+DTBYDZ2*(AKVS(L-1)+AKVS(L))) 4411. END DO 4412. DTBYDZ2 = 4d0*DTBYDZ(LMIJ)**2*BYDTS 4413. GZML(LMIJ,IQ,JQ)=GZML(LMIJ,IQ,JQ)/(1d0+DTBYDZ2*AKVG(LMIJ-1)) 4414. SZML(LMIJ,IQ,JQ)=SZML(LMIJ,IQ,JQ)/(1d0+DTBYDZ2*AKVS(LMIJ-1)) 4415. C**** Choice 2: too difficult 4416. C**** Choice 3 4417. C GZML(1,IQ,JQ) = 0. 4418. C SZML(1,IQ,JQ) = 0. 4419. C DO L=2,LMIJ-1 4420. C GZML(L,IQ,JQ) = 2.5d-1*(G0ML(L-1,IQ,JQ) - G0ML(L+1,IQ,JQ)) 4421. C SZML(L,IQ,JQ) = 2.5d-1*(S0ML(L-1,IQ,JQ) - S0ML(L+1,IQ,JQ)) 4422. C END DO 4423. C GZML(LMIJ,IQ,JQ)= 5d-1*(G0ML(LMIJ-1,IQ,JQ) - G0ML(LMIJ,IQ,JQ)) 4424. C SZML(LMIJ,IQ,JQ)= 5d-1*(S0ML(LMIJ-1,IQ,JQ) - S0ML(LMIJ,IQ,JQ)) 4425. C**** Diagnostics for non-local transport and vertical diffusion 4426. DO 530 L=1,LMIJ 4427. C OIJL(I,J,L,13) = OIJL(I,J,L,13) + AKVG(L)*GHATG(L) 4428. C OIJL(I,J,L,14) = OIJL(I,J,L,14) + AKVS(L)*GHATS(L) 4429. OIJL(I,J,L,13) = OIJL(I,J,L,13) + AKVM(L) ! need new common 4430. c OIJL(I,J,L,16) = OIJL(I,J,L,16) + AKVS(L) ! need new common 4431. 530 OIJL(I,J,L,14) = OIJL(I,J,L,14) + AKVG(L) 4432. C**** Set diagnostics 4433. C AIJ(I,J,77) = AIJ(I,J,77) + HBL ! boundary layer depth 4434. C AIJ(I,J,78) = AIJ(I,J,78) + Bo ! surface buoyancy forcing 4435. c AIJ(I,J,79) = AIJ(I,J,79) + Bosol ! solar buoyancy forcing 4436. c AIJ(I,J,80) = AIJ(I,J,80) + Ustar ! turbulent friction speed 4437. C**** 4438. C**** Update current prognostic variables 4439. C**** 4440. 700 DO 710 K=1,KMUV 4441. DO 710 L=1,LMUV(K) 4442. 710 UO(ID(K),1,L) = UO(ID(K),1,L) + RAMV(K)*(UL(L,K)-UT(ID(K),1,L)) 4443. 725 IF(QPOLE) GO TO 750 4444. C**** End of quarter box loops 4445. IQ=IQ+1 4446. IF(IQ.LE.2) GO TO 240 4447. JQ=JQ+1 4448. IF(JQ.LE.2) GO TO 230 4449. C**** 4450. C**** Recreate main prognostic variables of potential enthalpy and 4451. C**** salt from those on the quarter boxes 4452. C**** 4453. DO 720 L=1,LMIJ 4454. G0M(I,J,L)= G0ML(L,2,2)+G0ML(L,2,1)+G0ML(L,1,2)+G0ML(L,1,1) 4455. NSIGG = EXPONENT(G0M(I,J,L)) -2 - 42 4456. GXM(I,J,L)=(G0ML(L,2,2)+G0ML(L,2,1)-G0ML(L,1,2)-G0ML(L,1,1)) 4457. * *BYBETA 4458. GYM(I,J,L)=(G0ML(L,2,2)-G0ML(L,2,1)+G0ML(L,1,2)-G0ML(L,1,1)) 4459. * *BYBETA 4460. IF (NSIGG+23.gt.EXPONENT(GXM(I,J,L))) GXM(I,J,L) = 4461. * SCALE(REAL(NINT(SCALE(GXM(I,J,L),-NSIGG))),NSIGG) 4462. IF (NSIGG+23.gt.EXPONENT(GYM(I,J,L))) GYM(I,J,L) = 4463. * SCALE(REAL(NINT(SCALE(GYM(I,J,L),-NSIGG))),NSIGG) 4464. GZM(I,J,L)= GZML(L,2,2)+GZML(L,2,1)+GZML(L,1,2)+GZML(L,1,1) 4465. S0M(I,J,L)= S0ML(L,2,2)+S0ML(L,2,1)+S0ML(L,1,2)+S0ML(L,1,1) 4466. NSIGS = EXPONENT(S0M(I,J,L)) - 2 - 42 4467. SXM(I,J,L)=(S0ML(L,2,2)+S0ML(L,2,1)-S0ML(L,1,2)-S0ML(L,1,1)) 4468. * *BYBETA 4469. SYM(I,J,L)=(S0ML(L,2,2)-S0ML(L,2,1)+S0ML(L,1,2)-S0ML(L,1,1)) 4470. * *BYBETA 4471. IF (NSIGS+23.gt.EXPONENT(SXM(I,J,L))) SXM(I,J,L) = 4472. * SCALE(REAL(NINT(SCALE(SXM(I,J,L),-NSIGS))),NSIGS) 4473. IF (NSIGS+23.gt.EXPONENT(SYM(I,J,L))) SYM(I,J,L) = 4474. * SCALE(REAL(NINT(SCALE(SYM(I,J,L),-NSIGS))),NSIGS) 4475. 720 SZM(I,J,L)= SZML(L,2,2)+SZML(L,2,1)+SZML(L,1,2)+SZML(L,1,1) 4476. C**** End of I loop 4477. 730 IM1=I 4478. I=I+1 4479. IF(I.LE.IM) GO TO 210 4480. GO TO 790 4481. C**** 4482. C**** Load the prognostic variables of potential enthalpy and 4483. C**** salt from the column arrays at the poles 4484. C**** 4485. 750 DO 760 L=1,LMIJ 4486. G0M(1,JM,L) = G0ML(L,1,1) 4487. GZM(1,JM,L) = GZML(L,1,1) 4488. S0M(1,JM,L) = S0ML(L,1,1) 4489. 760 SZM(1,JM,L) = SZML(L,1,1) 4490. C**** End of outside J loop 4491. 790 CONTINUE 4492. RETURN 4493. END 4494. 4495. SUBROUTINE STCONV 4496. C**** 4497. C**** STCONV uses vertical diffusion coeeficients from STMIX 4498. C**** using KPP boundary layer code and (LCONV = TRUE) also applies 4499. C**** penetrating convection to each end of the straits. 4500. C**** 4501. INCLUDE 'C070.COM' 4502. COMMON /WORK01/ UL(LMO,2),MML(LMO),G0ML(LMO,2),S0ML(LMO,2) 4503. * ,GZML(LMO,2),SZML(LMO,2),BYMML(LMO),DTBYDZ(LMO) 4504. * ,BYDZ2(LMO),UL0(LMO),G0ML0(LMO),S0ML0(LMO) 4505. COMMON /DIFVCB/ AKVM(0:LMO+1),AKVG(0:LMO+1),AKVS(0:LMO+1) 4506. * ,GHAT(LMO) 4507. C**** CONV parameters: BETA controls degree of convection (default 0.5). 4508. PARAMETER (BETA=5d-1,BYBETA=1d0/BETA) 4509. LOGICAL*4 LCONV,LDD 4510. REAL*8 zgrid(0:LMO+1),hwide(0:LMO+1),Shsq(LMO),dVsq(LMO) 4511. * ,talpha(LMO),sbeta(LMO),dbloc(LMO),dbsfc(LMO),Ritop(LMO), 4512. * alphaDT(LMO),betaDS(LMO),byhwide(0:LMO+1) 4513. REAL*8 G(LMO),S(LMO),T(LMO),BYRHO(LMO),RHO(LMO),P(LMO) 4514. PARAMETER(epsln=1d-20,OMEGA=TWOPI*366d0/(365d0*SDAY)) 4514.5 DATA LCONV/.FALSE./,LDD/.FALSE./ 4515. DATA IFIRST/1/ 4516. 4517. IF (IFIRST.eq.1) THEN 4518. IFIRST=0 4519. zgrid(0) = epsln 4520. hwide(0) = epsln 4521. byhwide(0) = 0. 4522. DO L=1,LMO 4523. zgrid(L) = -0.5*(ZE(L-1) + ZE(L)) ! tracer level 4524. hwide(L) = zgrid(L-1) - zgrid(L) ! tracer level 4525. byhwide(L) = 1d0/hwide(L) 4526. END DO 4527. zgrid(LMO+1) = -ZE(LMO) 4528. hwide(LMO+1) = epsln 4529. byhwide(LMO+1) = 0. 4530. 4531. Ustar = 0. ! No wind/ice stress in straits 4532. Bo = 0 ! No buoyancy forcing 4533. Bosol = 0 4534. BYDTS = 1d0/DTS 4535. END IF 4536. C**** 4537. C**** Outside loop over straits 4538. C**** 4539. DO 790 N=1,NMST 4540. C**** 4541. C**** Define parameters and currents 4542. C**** 4543. LMIJ=LMST(N) 4544. P(1) = 5d-1*GRAV*MMST(1,N)/(DIST(N)*WIST(N)) 4545. DO 220 L=1,LMIJ-1 4546. P(L+1) = P(L) + 5d-1*GRAV*(MMST(L,N)+MMST(L+1,N)) 4547. * /(DIST(N)*WIST(N)) 4548. DTBYDZ(L) = DTS*DIST(N)*WIST(N)/MMST(L,N) 4549. BYDZ2(L) = 2d0*DIST(N)*WIST(N)/(MMST(L+1,N)+MMST(L,N)) 4550. MML(L) = 5d-1*MMST(L,N) 4551. 220 BYMML(L) = 1d0/MML(L) 4552. MML(LMIJ) = 5d-1*MMST(LMIJ,N) 4553. BYMML(LMIJ) = 1d0/MML(LMIJ) 4554. DTBYDZ(LMIJ) = DTS*DIST(N)*WIST(N)/MMST(LMIJ,N) 4555. C**** Loop over half boxes 4556. IQ=1 4557. 240 RI = BETA*(2d0*IQ-3d0) 4558. DO 250 L=1,LMIJ 4559. UL(L,IQ) = 5d-1* MUST(L,N)*DIST(N)*BYMML(L) 4560. G0ML(L,IQ) = 5d-1*(G0MST(L,N) + RI*GXMST(L,N)) 4561. S0ML(L,IQ) = 5d-1*(S0MST(L,N) + RI*SXMST(L,N)) 4562. GZML(L,IQ) = 5d-1* GZMST(L,N) 4563. 250 SZML(L,IQ) = 5d-1* SZMST(L,N) 4564. C**** 4565. C**** Vertical mixing derived from KPP scheme 4566. C**** 4567. 500 IF(LMIJ.LE.1) GO TO 700 4568. 4569. C**** Coriolis parameter, defined at strait end 4570. Coriol = 2d0*OMEGA*SINP(JST(N,IQ)) 4571. C**** Save original values 4572. UL0 = UL(1:LMO,IQ) 4573. G0ML0 = G0ML(1:LMO,IQ) 4574. S0ML0 = S0ML(1:LMO,IQ) 4575. C**** Start iteration for diffusivities 4576. ITER = 0 4577. HBL = 0 4578. 510 ITER =ITER + 1 4579. HBLP = HBL 4580. C**** Initiallise fields 4581. Shsq=0. 4582. dVsq=0. 4583. dbloc=0. 4584. dbsfc=0. 4585. Ritop=0. 4586. 4587. C**** velocity shears 4588. 4589. DO L=1,LMIJ-1 4590. Shsq(L) = (UL(L,IQ)-UL(L+1,IQ))**2 !interface 4591. dvsq(L) = (UL(1,IQ)-UL(L ,IQ))**2 !tracer pts 4592. END DO 4593. dVsq(LMIJ) = (UL(1,IQ)-UL(LMIJ,IQ))**2 ! tracer pts 4594. 4595. C**** density related quantities 4596. C**** 4597. C**** density of surface layer (kg/m3) 4598. C**** rho = rho{G1),S(1),zt(1)} 4599. C**** local buoyancy gradient at km interfaces: (m/s2) 4600. C**** dbloc = g/rho{k+1,k+1} * [ rho{k,k+1}-rho{k+1,k+1} ] 4601. C**** buoyancy difference with respect to "zref", the surface (m/s2) 4602. C**** dbsfc = g/rho{k,k} * [ rho{k,k}-rho{1,k} ] CORRECT 4603. C**** thermal expansion coefficient without 1/rho factor (kg/m3/C) 4604. C**** talpha= d(rho{k,k})/d(T(k)) 4605. C**** salt expansion coefficient without 1/rho factor (kg/m3/PSU) 4606. C**** sbeta = d(rho{k,k})/d(S(k)) 4607. 4608. G(1) = G0ML(1,IQ)*BYMML(1) 4609. S(1) = S0ML(1,IQ)*BYMML(1) 4610. BYRHO(1) = VOLGSP(G(1),S(1),P(1)) 4611. RHO(1) = 1d0/BYRHO(1) 4612. DO L=2,LMIJ 4613. G(L) = G0ML(L,IQ)*BYMML(L) 4614. S(L) = S0ML(L,IQ)*BYMML(L) 4615. BYRHO(L) = VOLGSP(G(L),S(L),P(L)) 4616. RHO(L) = 1d0/BYRHO(L) 4617. RHOM = 1d0/VOLGSP(G(L-1),S(L-1),P(L)) ! L-1 density wrt P(L) 4618. RHO1 = 1d0/VOLGSP(G(1),S(1),P(L)) ! surf density wrt P(L) 4619. dbsfc(L) = GRAV*(1d0-RHO1*BYRHO(L)) 4620. dbloc(L-1) = GRAV*(1d0-RHOM*BYRHO(L)) 4621. C**** numerator of bulk richardson number on grid levels 4622. Ritop(L) = (zgrid(1)-zgrid(L)) * dbsfc(L) 4623. END DO 4624. 4625. C**** Double diffusive option (ddmix) needs alphaDT,betaDS 4626. C**** alphaDT = mean -talpha * delta(temp.) at interfaces (kg/m3) 4627. C**** betaDS = mean sbeta * delta(salt) at interfaces (kg/m3) 4628. 4629. if (LDD) then 4630. T(1) = TEMGSP(G(1),S(1),P(1)) 4631. talpha(1) = ALPHAGSP(G(1),S(1),P(1)) ! >0 4632. sbeta(1) = BETAGSP(G(1),S(1),P(1)) 4633. do L=2,LMIJ 4634. T(L) = TEMGSP(G(L),S(L),P(L)) 4635. talpha(L) = ALPHAGSP(G(L),S(L),P(L)) ! >0 4636. sbeta(L) = BETAGSP(G(L),S(L),P(L)) 4637. alphaDT(L-1) = - 5d-1 * (talpha(L-1) + talpha(L)) 4638. $ * (T(L-1) - T(L)) 4639. betaDS (L-1) = 5d-1 * (sbeta (L-1) + sbeta(L)) 4640. $ * (S(L-1) - S(L))*1d3 4641. end do 4642. end if 4643. 4644. C**** Get diffusivities for the whole column 4645. 4646. CALL KPPMIX(LDD,LMO,zgrid,hwide,LMIJ,Shsq,dVsq,Ustar,Bo,Bosol 4647. * ,alphaDT,betaDS,dbloc,Ritop,Coriol,byhwide, 4648. * AKVM,AKVS,AKVG,GHAT,HBL,KBL) 4649. 4650. C**** Correct units for diffusivities (m^2/s) => (kg^2/m^4 s) 4651. DO L=1,LMIJ-1 4652. R = 5d-1*(RHO(L)+RHO(L+1)) 4653. R2 = R**2 4654. AKVM(L) = AKVM(L)*R2 4655. AKVG(L) = AKVG(L)*R2 4656. AKVS(L) = AKVS(L)*R2 4657. GHAT(L) = 0. ! no non-local transports since no surface fluxes 4658. END DO 4659. 4660. C**** For each field (U,G,S + TRACERS) call VDIFF 4661. C**** Momentum 4662. CALL VDIFF(UL(1,IQ),AKVM(1),GHAT,DTBYDZ,BYDZ2,DTS 4663. * ,LMIJ,UL0) 4664. C**** Enthalpy 4665. CALL VDIFFS(G0ML(1,IQ),AKVG(1),GHAT,DTBYDZ,BYDZ2,DTS 4666. * ,LMIJ,G0ML0) 4667. C**** Salinity 4668. CALL VDIFFS(S0ML(1,IQ),AKVS(1),GHAT,DTBYDZ,BYDZ2,DTS 4669. * ,LMIJ,S0ML0) 4670. IF ((ITER.eq.1 .or. ABS(HBLP-HBL).gt.(ZE(KBL)-ZE(KBL-1))*0.02 ) 4671. * .and. ITER.lt.4) GO TO 510 4672. C**** Gradients of scalars? 4673. C**** Choice 1B: Same as above but implicitly 4674. GZML(1,IQ) = 0. 4675. SZML(1,IQ) = 0. 4676. DO L=2,LMIJ-1 4677. DTBYDZ2 = 2d0*DTBYDZ(L)**2*BYDTS 4678. GZML(L,IQ)=GZML(L,IQ)/(1d0+DTBYDZ2*(AKVG(L-1)+AKVG(L))) 4679. SZML(L,IQ)=SZML(L,IQ)/(1d0+DTBYDZ2*(AKVS(L-1)+AKVS(L))) 4680. END DO 4681. DTBYDZ2 = 4d0*DTBYDZ(LMIJ)**2*BYDTS 4682. GZML(LMIJ,IQ)=GZML(LMIJ,IQ)/(1d0+DTBYDZ2*AKVG(LMIJ-1)) 4683. SZML(LMIJ,IQ)=SZML(LMIJ,IQ)/(1d0+DTBYDZ2*AKVS(LMIJ-1)) 4684. DO L=1,LMIJ 4685. C OLNST(L,N,4) = OLNST(L,N,4) + AKVS(L) ! new common 4686. C OLNST(L,N,6) = OLNST(L,N,6) + AKVG(L) ! needed 4687. OLNST(L,N,5) = OLNST(L,N,5) + AKVM(L) 4688. END DO 4689. C**** Set diagnostics (no such things!) 4690. C AST(N,1) = AST(N,1) + HBL ! boundary layer depth 4691. 4692. C**** End of half box loops 4693. 700 IQ=IQ+1 4694. IF(IQ.LE.2) GO TO 240 4695. C**** 4696. C**** Recreate main prognostic variables of potential enthalpy and 4697. C**** salt from those on the half boxes 4698. C**** 4699. DO 720 L=1,LMIJ 4700. MUST(L,N) = ( UL(L,2) + UL(L,1))*MML(L)/DIST(N) 4701. G0MST(L,N) = G0ML(L,2) + G0ML(L,1) 4702. NSIGG = EXPONENT(G0MST(L,N)) -1 - 42 4703. GXMST(L,N) = (G0ML(L,2) - G0ML(L,1))*BYBETA 4704. IF (NSIGG+30.gt.EXPONENT(GXMST(L,N))) GXMST(L,N) = 4705. * SCALE(dble(NINT(SCALE(GXMST(L,N),-NSIGG))),NSIGG) 4706. GZMST(L,N) = GZML(L,2) + GZML(L,1) 4707. S0MST(L,N) = S0ML(L,2) + S0ML(L,1) 4708. NSIGS = EXPONENT(S0MST(L,N)) -1 - 42 4709. SXMST(L,N) = (S0ML(L,2) - S0ML(L,1))*BYBETA 4710. IF (NSIGS+30.gt.EXPONENT(SXMST(L,N))) SXMST(L,N) = 4711. * SCALE(dble(NINT(SCALE(SXMST(L,N),-NSIGS))),NSIGS) 4712. 720 SZMST(L,N) = SZML(L,2) + SZML(L,1) 4713. C**** End of outside loop over straits 4714. 790 CONTINUE 4715. RETURN 4716. END 9300. 9301. SUBROUTINE VDIFF(U,K,GHAT,DTBYDZ,BYDZ2,DT,LMIJ,U0) 9302. C**** Implicit vertical diffusion + non local transport for velocity 9302.1 PARAMETER (LMO=13) 9303. REAL*8 U(LMO),K(LMO),GHAT(LMO),DTBYDZ(LMO),BYDZ2(LMO) 9304. REAL*8 A(LMO),B(LMO),C(LMO),R(LMO),DT,U0(LMO),gam(LMO),bet 9304.1 INTENT (IN) K,GHAT,DTBYDZ,BYDZ2,DT,LMIJ,U0 9304.2 INTENT (OUT) U 9305. C**** U0,U input and output field (velocity or concentration) 9306. C**** K vertical diffusivity ((z)^2/ s) 9307. C**** GHAT non local transport of scalar 9308. C**** kv * ghats * surface flux 9309. C**** DTBYDZ DT/DZ_L 9310. C**** BYDZ2 1d0/DZ_L+1/2 9311. C**** DT timestep (s) 9312. C**** Boundary conditions assumed to be no-flux at Z=0, Z=Z(LMIJ) 9313. C**** Calculate operators for tridiagonal solver 9314. A(1) = 0 9315. B(1) = 1d0 + DTBYDZ(1)*BYDZ2(1)*K(1) 9316. C(1) = - DTBYDZ(1)*BYDZ2(1)*K(1) 9317. R(1) = U0(1) - DTBYDZ(1)*GHAT(1) 9318. DO L=2,LMIJ-1 9319. A(L) = - DTBYDZ(L)* BYDZ2(L-1)*K(L-1) 9320. B(L) = 1d0 + DTBYDZ(L)*(BYDZ2(L-1)*K(L-1)+BYDZ2(L)*K(L)) 9321. C(L) = - DTBYDZ(L)* BYDZ2(L)*K(L) 9322. R(L) = U0(L) + DTBYDZ(L)*(GHAT(L-1) - GHAT(L)) 9323. END DO 9324. A(LMIJ) = - DTBYDZ(LMIJ)*BYDZ2(LMIJ-1)*K(LMIJ-1) 9325. B(LMIJ) = 1d0 + DTBYDZ(LMIJ)*BYDZ2(LMIJ-1)*K(LMIJ-1) 9326. C(LMIJ) = 0 9327. R(LMIJ) = U0(LMIJ) + DTBYDZ(LMIJ)*GHAT(LMIJ-1) 9328. 9329. C CALL TRIDIAG(A,B,C,R,U,LMIJ) ! INLINED 9408. bet=b(1) 9409. u(1)=r(1)/bet 9410. do j=2,LMIJ 9411. gam(j)=c(j-1)/bet 9412. bet=b(j)-a(j)*gam(j) 9414. u(j)=(r(j)-a(j)*u(j-1))/bet 9415. end do 9416. do j=LMIJ-1,1,-1 9417. u(j)=u(j)-gam(j+1)*u(j+1) 9418. end do 9330. 9331. RETURN 9332. END 9350. 9351. SUBROUTINE VDIFFS(U,K,GHAT,DTBYDZ,BYDZ2,DT,LMIJ,U0) 9352. C**** Implicit vertical diffusion + non local transport for tracers 9352.1 PARAMETER (LMO=13) 9353. REAL*8 U(LMO),K(LMO),GHAT(LMO),DTBYDZ(LMO),BYDZ2(LMO) 9354. REAL*8 A(LMO),B(LMO),C(LMO),R(LMO),DT,U0(LMO),gam(LMO),bet 9354.1 INTENT (IN) K,GHAT,DTBYDZ,BYDZ2,DT,LMIJ,U0 9354.2 INTENT (OUT) U 9355. C**** U0,U input and output fields (total tracer) 9356. C**** K vertical diffusivity ((z)^2/ s) 9357. C**** GHAT non local transport of scalar 9358. C**** kv * ghats * surface flux 9359. C**** DTBYDZ DT/DZ_L 9360. C**** BYDZ2 1d0/DZ_L+1/2 9361. C**** DT timestep (s) 9362. C**** Boundary conditions assumed to be no-flux at Z=0, Z=Z(LMIJ) 9363. C**** Calculate operators for tridiagonal solver 9364. A(1) = 0 9365. B(1) = 1d0 + DTBYDZ(1)*BYDZ2(1)*K(1) 9366. C(1) = - DTBYDZ(2)*BYDZ2(1)*K(1) 9367. R(1) = U0(1) - DT * GHAT(1) 9368. DO L=2,LMIJ-1 9369. A(L) = - DTBYDZ(L-1)* BYDZ2(L-1)*K(L-1) 9370. B(L) = 1d0 + DTBYDZ(L )*(BYDZ2(L-1)*K(L-1)+BYDZ2(L)*K(L)) 9371. C(L) = - DTBYDZ(L+1)* BYDZ2(L)*K(L) 9372. R(L) = U0(L) + DT * (GHAT(L-1) - GHAT(L)) 9373. END DO 9374. A(LMIJ) = - DTBYDZ(LMIJ-1)*BYDZ2(LMIJ-1)*K(LMIJ-1) 9375. B(LMIJ) = 1d0 + DTBYDZ(LMIJ )*BYDZ2(LMIJ-1)*K(LMIJ-1) 9376. C(LMIJ) = 0 9377. R(LMIJ) = U0(LMIJ) + DT * GHAT(LMIJ-1) 9378. 9379. C CALL TRIDIAG(A,B,C,R,U,LMIJ) ! INLINED 9408. bet=b(1) 9409. u(1)=r(1)/bet 9410. do j=2,LMIJ 9411. gam(j)=c(j-1)/bet 9412. bet=b(j)-a(j)*gam(j) 9414. u(j)=(r(j)-a(j)*u(j-1))/bet 9415. end do 9416. do j=LMIJ-1,1,-1 9417. u(j)=u(j)-gam(j+1)*u(j+1) 9418. end do 9380. 9381. RETURN 9382. END 9400. 9401. SUBROUTINE TRIDIAG(A,B,C,R,U,N) 9402. C**** Solve a tridiagonal system (A,B,C)U = R 9403. C**** Routine taken from Numerical Recipes (Press et al, 1990) 9403.1 INTENT (IN) A,B,C,R,N 9403.2 INTENT (OUT) U 9404. INTEGER n,NMAX,j 9405. PARAMETER (NMAX=5000) 9406. REAL*8 A(N),B(N),C(N),R(N),U(N),bet,gam(NMAX) 9407. IF (B(1).eq.0) STOP "Rewrite equations in TRIDIAG" 9408. bet=b(1) 9409. u(1)=r(1)/bet 9410. do j=2,n 9411. gam(j)=c(j-1)/bet 9412. bet=b(j)-a(j)*gam(j) 9413. C if (bet.eq.0) STOP "TRIDIAG failed" 9414. u(j)=(r(j)-a(j)*u(j-1))/bet 9415. end do 9416. do j=n-1,1,-1 9417. u(j)=u(j)-gam(j+1)*u(j+1) 9418. end do 9419. return 9420. end 9421.