New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps (diff) – NEMO

Changes between Version 16 and Version 17 of 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps


Ignore:
Timestamp:
2019-03-04T18:50:27+01:00 (5 years ago)
Author:
acc
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps

    v16 v17  
    8181{{{#!Fold title=Preview tag=h2 
    8282[=#preview] 
     83[#step1 The refactoring script] \\ 
     84[#step2 The refactoring script explained]\\ 
     85[#step3 Notes on testing regular expressions]\\ 
     86[#step4 Some contrived tests]\\ 
     87[#step5 Results of real files]\\ 
     88\\ 
     89[#step6 Automating the tiling changes]\\ 
     90[#step7 Results on traldf_iso.F90]\\ 
     91[#step8 do2dfinder.pl]\\ 
     92[#step9 do3dfinder.pl]\\ 
     93[#step10 Sanity checks and domain_substitute.h90]\\ 
     94 
     95 
    8396 
    8497Part of the reorganisation for RK3 requires the refactoring of arrays such as un, ub into a single, 4 dimensional array with a time-level dimension. It is expected that much of the work required here can be automated to the extent that it is feasible to re-apply these changes after the annual merge. Below is a working example of how this might be achieved. Perl is used to carry out the pattern matching and substitution because of its ability to match patterns extending over several lines. A random subset of source files are used in this example and serve to illustrate the successes and caveats for the method.  
    … …  
    96109}}} 
    97110 
    98 '''The refactoring script''' 
     111[=#step1 '''The refactoring script'''] 
    99112 
    100113{{{#!bash 
    … …  
    121134}}} 
    122135 
    123 '''The refactoring script explained''' 
     136[=#step2 '''The refactoring script explained'''] 
    124137 
    125138{{{#!bash 
    … …  
    175188}}} 
    176189 
    177 '''Notes on testing regular expressions''' 
     190[=#step3 '''Notes on testing regular expressions'''] 
    178191 
    179192Testing and deciphering the regular expression used in the LHS of the perl substitute command is made easier by the availability of on-line testers. below is a screenshot from regex101.com which helps illustrate and explain the regular expression used here: 
    … …  
    181194[[Image(regex101_example_sm.png)]] 
    182195 
    183 ''' Some contrived tests:''' 
     196[=#step4 ''' Some contrived tests:'''] 
    184197 
    185198{{{#!f 
    … …  
    256269So all changes were made correctly and even those entries which were potential pitfalls (pun and sbc_fwb) were correctly ignored. Time to try a real set: 
    257270 
    258 '''The results on the sample set of files (patch.list):''' 
     271[=#step5 '''The results on the sample set of files (patch.list):'''] 
    259272 
    260273{{{#!diff 
    … …  
    563576So far so good.... 
    564577 
    565 == Automating the tiling changes 
     578== [=#step6 Automating the tiling changes] 
    566579 
    567580Here is a almost complete attempt at automating the loop changes. Earlier versions (now superceded) maintained the DO loop ranges as arguments to the macros. These arguments are now interptreted and converted to the binary representative form suggested by Gurvan. The logic for this is basic at present and possibly easily fooled (but works on the examples used so far). I've persisted with a two-stage conversion with a script to convert 2D loops and then a second script to convert 3D loops. This makes the scripts readable and allows easier verification. The two scripts are named `do2dfinder.pl` and `do3dfinder.pl` and are included below. Firstly here is an example of the scripts in action on the following test file: 
    … …  
    665678+     END_3D 
    666679}}} 
    667 traldf_iso.F90 provides a more stringent test: 
     680[=#step7 traldf_iso.F90 provides a more stringent test:] 
    668681{{{ 
    669682   perl do2dfinder.pl TESTDO_FILES/traldf_iso.F90 > TESTDO_FILES_2D/traldf_iso.F90 
    … …  
    984997              ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==! 
    985998}}} 
    986 And finally the two scripts that achieve this. Note the logic that may need tightening at sections 5 and 6 in the first script: 
     999[=#step8 And finally the two scripts that achieve this. Note the logic that may need tightening at sections 5 and 6 in the first script:] 
    9871000{{{#!perl 
    9881001#cat do2dfinder.pl 
    … …  
    11081121} 
    11091122}}} 
    1110 and 
     1123[=#step9 and] 
    11111124{{{#!perl 
    11121125#cat do3dfinder.pl 
    … …  
    12181231} 
    12191232}}} 
     1233== [=#step10 Sanity Check] 
     1234 
     1235Introducing 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{{{ 
     1238cat 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 
     1256ed - TESTDO_FILES_3D/traldf_iso.F90 << EOF 
     1257> 0a 
     1258> #include "domain_substute.h90" 
     1259> . 
     1260> w 
     1261> q 
     1262> EOF 
     1263 
     1264gfortran -E -P TESTDO_FILES_3D/traldf_iso.F90 > SANITY/traldf_iso.f90 
     1265diff -u TESTDO_FILES/traldf_iso.F90 SANITY/traldf_iso.f90 > sanity.patch 
     1266 
     1267}}} 
     1268 
     1269And 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}}} 
    12201609 
    12211610