| 1233 | == [=#step10 Sanity Check] |
| 1234 | |
| 1235 | Introducing a form of the proposed `domain_substitute.h90` file to the final version and running through a preprocssor should recover the equivalent of the original file (barring white space changes and line concatenations). For example: |
| 1236 | |
| 1237 | {{{ |
| 1238 | cat domain_substitute.h90 |
| 1239 | #define kJs 2 |
| 1240 | #define kIs 2 |
| 1241 | #define kJe jpjm1 |
| 1242 | #define kIe jpim1 |
| 1243 | #define DO_2D_00_00 DO jj = kJs ,kJe ; DO ji = kIs ,kIe |
| 1244 | #define DO_2D_10_10 DO jj = kJs-1, kJe ; DO ji = kIs-1, kIe |
| 1245 | #define DO_2D_01_01 DO jj = kJs , kJe+1 ; DO ji = kIs , kIe+1 |
| 1246 | #define DO_2D_11_11 DO jj = kJs-1, kJe+1 ; DO ji = kIs-1, kIe+1 |
| 1247 | |
| 1248 | #define DO_3D_00_00(ks,ke) DO jk = ks, ke ; DO_2D_00_00 |
| 1249 | #define DO_3D_10_10(ks,ke) DO jk = ks, ke ; DO_2D_10_10 |
| 1250 | #define DO_3D_01_01(ks,ke) DO jk = ks, ke ; DO_2D_01_01 |
| 1251 | #define DO_3D_11_11(ks,ke) DO jk = ks, ke ; DO_2D_11_11 |
| 1252 | |
| 1253 | #define END_2D END DO ; END DO |
| 1254 | #define END_3D END DO ; END DO ; END DO |
| 1255 | |
| 1256 | ed - TESTDO_FILES_3D/traldf_iso.F90 << EOF |
| 1257 | > 0a |
| 1258 | > #include "domain_substute.h90" |
| 1259 | > . |
| 1260 | > w |
| 1261 | > q |
| 1262 | > EOF |
| 1263 | |
| 1264 | gfortran -E -P TESTDO_FILES_3D/traldf_iso.F90 > SANITY/traldf_iso.f90 |
| 1265 | diff -u TESTDO_FILES/traldf_iso.F90 SANITY/traldf_iso.f90 > sanity.patch |
| 1266 | |
| 1267 | }}} |
| 1268 | |
| 1269 | And this does appear to be the case: |
| 1270 | |
| 1271 | {{{#!diff |
| 1272 | --- TESTDO_FILES/traldf_iso.F90 2019-03-04 12:59:44.000000000 +0000 |
| 1273 | +++ SANITY/traldf_iso.f90 2019-03-04 17:17:32.000000000 +0000 |
| 1274 | @@ -1,3 +1,6 @@ |
| 1275 | + |
| 1276 | + |
| 1277 | + |
| 1278 | MODULE traldf_iso |
| 1279 | !!====================================================================== |
| 1280 | !! *** MODULE traldf_iso *** |
| 1281 | @@ -39,7 +42,17 @@ |
| 1282 | LOGICAL :: l_hst ! flag to compute heat transport |
| 1283 | |
| 1284 | !! * Substitutions |
| 1285 | -# include "vectopt_loop_substitute.h90" |
| 1286 | + !!---------------------------------------------------------------------- |
| 1287 | + !! *** vectopt_loop_substitute *** |
| 1288 | + !!---------------------------------------------------------------------- |
| 1289 | + !! ** purpose : substitute the inner loop start/end indices with CPP macro |
| 1290 | + !! allow unrolling of do-loop (useful with vector processors) |
| 1291 | + !!---------------------------------------------------------------------- |
| 1292 | + !!---------------------------------------------------------------------- |
| 1293 | + !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
| 1294 | + !! $Id: vectopt_loop_substitute.h90 10068 2018-08-28 14:09:04Z nicolasmartin $ |
| 1295 | + !! Software governed by the CeCILL license (see ./LICENSE) |
| 1296 | + !!---------------------------------------------------------------------- |
| 1297 | !!---------------------------------------------------------------------- |
| 1298 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
| 1299 | !! $Id: traldf_iso.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ |
| 1300 | @@ -143,58 +156,42 @@ |
| 1301 | ! |
| 1302 | IF( kpass == 1 ) THEN !== first pass only ==! |
| 1303 | ! |
| 1304 | - DO jk = 2, jpkm1 |
| 1305 | - DO jj = 2, jpjm1 |
| 1306 | - DO ji = fs_2, fs_jpim1 ! vector opt. |
| 1307 | - ! |
| 1308 | - zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & |
| 1309 | - & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) |
| 1310 | - zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & |
| 1311 | - & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) |
| 1312 | - ! |
| 1313 | - zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & |
| 1314 | - & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku |
| 1315 | - zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & |
| 1316 | - & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv |
| 1317 | - ! |
| 1318 | - ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & |
| 1319 | - & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) |
| 1320 | - END DO |
| 1321 | - END DO |
| 1322 | - END DO |
| 1323 | + DO jk = 2, jpkm1 ; DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 |
| 1324 | + ! |
| 1325 | + zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & |
| 1326 | + & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) |
| 1327 | + zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & |
| 1328 | + & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) |
| 1329 | + ! |
| 1330 | + zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & |
| 1331 | + & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku |
| 1332 | + zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & |
| 1333 | + & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv |
| 1334 | + ! |
| 1335 | + ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & |
| 1336 | + & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) |
| 1337 | + END DO ; END DO ; END DO |
| 1338 | ! |
| 1339 | IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient |
| 1340 | - DO jk = 2, jpkm1 |
| 1341 | - DO jj = 2, jpjm1 |
| 1342 | - DO ji = fs_2, fs_jpim1 |
| 1343 | - akz(ji,jj,jk) = 0.25_wp * ( & |
| 1344 | - & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & |
| 1345 | - & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & |
| 1346 | - & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & |
| 1347 | - & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) |
| 1348 | - END DO |
| 1349 | - END DO |
| 1350 | - END DO |
| 1351 | + DO jk = 2, jpkm1 ; DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 |
| 1352 | + akz(ji,jj,jk) = 0.25_wp * ( & |
| 1353 | + & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & |
| 1354 | + & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & |
| 1355 | + & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & |
| 1356 | + & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) |
| 1357 | + END DO ; END DO ; END DO |
| 1358 | ! |
| 1359 | IF( ln_traldf_blp ) THEN ! bilaplacian operator |
| 1360 | - DO jk = 2, jpkm1 |
| 1361 | - DO jj = 1, jpjm1 |
| 1362 | - DO ji = 1, fs_jpim1 |
| 1363 | - akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & |
| 1364 | - & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) |
| 1365 | - END DO |
| 1366 | - END DO |
| 1367 | - END DO |
| 1368 | + DO jk = 2, jpkm1 ; DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 |
| 1369 | + akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & |
| 1370 | + & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) |
| 1371 | + END DO ; END DO ; END DO |
| 1372 | ELSEIF( ln_traldf_lap ) THEN ! laplacian operator |
| 1373 | - DO jk = 2, jpkm1 |
| 1374 | - DO jj = 1, jpjm1 |
| 1375 | - DO ji = 1, fs_jpim1 |
| 1376 | - ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) |
| 1377 | - zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) |
| 1378 | - akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt |
| 1379 | - END DO |
| 1380 | - END DO |
| 1381 | - END DO |
| 1382 | + DO jk = 2, jpkm1 ; DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 |
| 1383 | + ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) |
| 1384 | + zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) |
| 1385 | + akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt |
| 1386 | + END DO ; END DO ; END DO |
| 1387 | ENDIF |
| 1388 | ! |
| 1389 | ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit |
| 1390 | @@ -215,28 +212,20 @@ |
| 1391 | !!end |
| 1392 | |
| 1393 | ! Horizontal tracer gradient |
| 1394 | - DO jk = 1, jpkm1 |
| 1395 | - DO jj = 1, jpjm1 |
| 1396 | - DO ji = 1, fs_jpim1 ! vector opt. |
| 1397 | - zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) |
| 1398 | - zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) |
| 1399 | - END DO |
| 1400 | - END DO |
| 1401 | - END DO |
| 1402 | + DO jk = 1, jpkm1 ; DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 |
| 1403 | + zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) |
| 1404 | + zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) |
| 1405 | + END DO ; END DO ; END DO |
| 1406 | IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient |
| 1407 | - DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) |
| 1408 | - DO ji = 1, fs_jpim1 ! vector opt. |
| 1409 | - zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) |
| 1410 | - zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) |
| 1411 | - END DO |
| 1412 | - END DO |
| 1413 | + DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 |
| 1414 | + zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) |
| 1415 | + zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) |
| 1416 | + END DO ; END DO |
| 1417 | IF( ln_isfcav ) THEN ! first wet level beneath a cavity |
| 1418 | - DO jj = 1, jpjm1 |
| 1419 | - DO ji = 1, fs_jpim1 ! vector opt. |
| 1420 | - IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) |
| 1421 | - IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) |
| 1422 | - END DO |
| 1423 | - END DO |
| 1424 | + DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 |
| 1425 | + IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) |
| 1426 | + IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) |
| 1427 | + END DO ; END DO |
| 1428 | ENDIF |
| 1429 | ENDIF |
| 1430 | ! |
| 1431 | @@ -252,36 +241,32 @@ |
| 1432 | IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) |
| 1433 | ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) |
| 1434 | ENDIF |
| 1435 | - DO jj = 1 , jpjm1 !== Horizontal fluxes |
| 1436 | - DO ji = 1, fs_jpim1 ! vector opt. |
| 1437 | - zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) |
| 1438 | - zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) |
| 1439 | - ! |
| 1440 | - zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & |
| 1441 | - & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) |
| 1442 | - ! |
| 1443 | - zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & |
| 1444 | - & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) |
| 1445 | - ! |
| 1446 | - zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku |
| 1447 | - zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv |
| 1448 | - ! |
| 1449 | - zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & |
| 1450 | - & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & |
| 1451 | - & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) |
| 1452 | - zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & |
| 1453 | - & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & |
| 1454 | - & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) |
| 1455 | - END DO |
| 1456 | - END DO |
| 1457 | + DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 |
| 1458 | + zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) |
| 1459 | + zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) |
| 1460 | + ! |
| 1461 | + zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & |
| 1462 | + & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) |
| 1463 | + ! |
| 1464 | + zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & |
| 1465 | + & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) |
| 1466 | + ! |
| 1467 | + zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku |
| 1468 | + zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv |
| 1469 | + ! |
| 1470 | + zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & |
| 1471 | + & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & |
| 1472 | + & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) |
| 1473 | + zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & |
| 1474 | + & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & |
| 1475 | + & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) |
| 1476 | + END DO ; END DO |
| 1477 | ! |
| 1478 | - DO jj = 2 , jpjm1 !== horizontal divergence and add to pta |
| 1479 | - DO ji = fs_2, fs_jpim1 ! vector opt. |
| 1480 | - pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & |
| 1481 | - & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & |
| 1482 | - & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) |
| 1483 | - END DO |
| 1484 | - END DO |
| 1485 | + DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 |
| 1486 | + pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & |
| 1487 | + & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & |
| 1488 | + & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) |
| 1489 | + END DO ; END DO |
| 1490 | END DO ! End of slab |
| 1491 | |
| 1492 | !!---------------------------------------------------------------------- |
| 1493 | @@ -295,75 +280,55 @@ |
| 1494 | ! ! Surface and bottom vertical fluxes set to zero |
| 1495 | ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp |
| 1496 | |
| 1497 | - DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) |
| 1498 | - DO jj = 2, jpjm1 |
| 1499 | - DO ji = fs_2, fs_jpim1 ! vector opt. |
| 1500 | - ! |
| 1501 | - zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & |
| 1502 | - & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) |
| 1503 | - zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & |
| 1504 | - & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) |
| 1505 | - ! |
| 1506 | - zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & |
| 1507 | - & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku |
| 1508 | - zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & |
| 1509 | - & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv |
| 1510 | - ! |
| 1511 | - zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked |
| 1512 | - zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) |
| 1513 | - ! |
| 1514 | - ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & |
| 1515 | - & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & |
| 1516 | - & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & |
| 1517 | - & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) |
| 1518 | - END DO |
| 1519 | - END DO |
| 1520 | - END DO |
| 1521 | + DO jk = 2, jpkm1 ; DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 |
| 1522 | + ! |
| 1523 | + zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & |
| 1524 | + & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) |
| 1525 | + zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & |
| 1526 | + & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) |
| 1527 | + ! |
| 1528 | + zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & |
| 1529 | + & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku |
| 1530 | + zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & |
| 1531 | + & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv |
| 1532 | + ! |
| 1533 | + zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked |
| 1534 | + zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) |
| 1535 | + ! |
| 1536 | + ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & |
| 1537 | + & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & |
| 1538 | + & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & |
| 1539 | + & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) |
| 1540 | + END DO ; END DO ; END DO |
| 1541 | ! !== add the vertical 33 flux ==! |
| 1542 | IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz |
| 1543 | - DO jk = 2, jpkm1 |
| 1544 | - DO jj = 1, jpjm1 |
| 1545 | - DO ji = fs_2, fs_jpim1 |
| 1546 | - ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & |
| 1547 | - & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & |
| 1548 | - & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) |
| 1549 | - END DO |
| 1550 | - END DO |
| 1551 | - END DO |
| 1552 | + DO_3D_10_00( 2, jpkm1 ) |
| 1553 | + ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & |
| 1554 | + & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & |
| 1555 | + & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) |
| 1556 | + END DO ; END DO ; END DO |
| 1557 | ! |
| 1558 | ELSE ! bilaplacian |
| 1559 | SELECT CASE( kpass ) |
| 1560 | CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 |
| 1561 | - DO jk = 2, jpkm1 |
| 1562 | - DO jj = 1, jpjm1 |
| 1563 | - DO ji = fs_2, fs_jpim1 |
| 1564 | - ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & |
| 1565 | - & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & |
| 1566 | - & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) |
| 1567 | - END DO |
| 1568 | - END DO |
| 1569 | - END DO |
| 1570 | + DO_3D_10_00( 2, jpkm1 ) |
| 1571 | + ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & |
| 1572 | + & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & |
| 1573 | + & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) |
| 1574 | + END DO ; END DO ; END DO |
| 1575 | CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. |
| 1576 | - DO jk = 2, jpkm1 |
| 1577 | - DO jj = 1, jpjm1 |
| 1578 | - DO ji = fs_2, fs_jpim1 |
| 1579 | - ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & |
| 1580 | - & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & |
| 1581 | - & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) |
| 1582 | - END DO |
| 1583 | - END DO |
| 1584 | - END DO |
| 1585 | + DO_3D_10_00( 2, jpkm1 ) |
| 1586 | + ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & |
| 1587 | + & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & |
| 1588 | + & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) |
| 1589 | + END DO ; END DO ; END DO |
| 1590 | END SELECT |
| 1591 | ENDIF |
| 1592 | ! |
| 1593 | - DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! |
| 1594 | - DO jj = 2, jpjm1 |
| 1595 | - DO ji = fs_2, fs_jpim1 ! vector opt. |
| 1596 | - pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & |
| 1597 | - & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) |
| 1598 | - END DO |
| 1599 | - END DO |
| 1600 | - END DO |
| 1601 | + DO jk = 1, jpkm1 ; DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 |
| 1602 | + pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & |
| 1603 | + & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) |
| 1604 | + END DO ; END DO ; END DO |
| 1605 | ! |
| 1606 | IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==! |
| 1607 | ( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass (bilaplacian) ==! |
| 1608 | }}} |