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 11 and Version 12 of 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps


Ignore:
Timestamp:
2019-03-01T15:45:10+01:00 (5 years ago)
Author:
acc
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps

    v11 v12  
    563563So far so good.... 
    564564 
    565 '''Automating the tiling changes''' 
     565== Automating the tiling changes 
    566566 
    567567 
    568568Here is a first attempt at automating the loop changes. Just for the 2D loops so far but it shows the possibilities. For now, I've assumed we don't need all the explicit DO_2D_00_01 type macros but just go for a generic version with arguments. This was proposed for the 3D version so why not for the 2D cases?. TBD. Firstly here is the annotated perl script: 
    569569{{{ 
     570cat do2dfinder.pl 
    570571# 
    571572open(F,$ARGV[0]) || die "Cannot open $ARGV[0]: $!"; 
    … …  
    885886       END DO                                           !   End of slab 
    886887       !                                                ! =============== 
    887 @@ -365,43 +347,33 @@ 
    888           CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    889              zwz(:,:) = ff_f(:,:) 
    890           CASE ( np_RVO )                           !* relative vorticity 
    891 -            DO jj = 1, jpjm1 
    892 -               DO ji = 1, fs_jpim1   ! vector opt. 
    893 -                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    894 -                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    895 -               END DO 
     888 
     889}}} 
     890 
     891''' Extending for the 3D loops ''' 
     892 
     893Dealing with the 2D loops can be considered as the first stage of converting suitable 3D loops. In the simplest cases a small variation of the do2dfinder.pl perl script which subsequently looks for consecutive `DO jk =` `DO_2D` statements instead of `DO jj`-`DO ji` pairs should work. However, finding such consecutive statements, in the 3D case, is a much less certain indication of a valid loop to process. Take, for example, this snippet of `dynadv_cen2.F90` which has been processed for 2D loops: 
     894 
     895{{{ 
     896      DO jk = 2, jpkm1                    ! interior advective fluxes 
     897         DO_2D( 2, jpi , 2, jpj ) 
     898            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     899         END_2D 
     900         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     901            zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
     902            zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     903         END_2D 
     904      END DO 
     905      DO jk = 1, jpkm1                    ! divergence of vertical momentum flux divergence 
     906         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     907            ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     908            va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     909         END_2D 
     910      END DO 
     911}}} 
     912 
     913Only the second 3D loop can be collapsed to a single construct. Unfortunately, this decision can't be made until the end of the loop is reached so a modification is required that allows two versions of a loop to be mantained in memory and the appropriate set written out on decision. The IO:stringy perl package provides the ideal tools by allowing string variables to be treated like files. Here is the modified script that will apply the 2nd stage conversion of 3D loops were possible: 
     914 
     915{{{ 
     916cat do3dfinder.pl 
     917# 
     918use IO::Scalar; 
     919open(F,$ARGV[0]) || die "Cannot open $ARGV[0]: $!"; 
     920while(<F>) { 
     921   if ( $_ =~ /^\s*DO\s* jk/i) { 
     922      # Start processing loop if line contains DO jk (case and whitespace independent) 
     923      # 
     924      # 1. Store the current line 
     925      # 
     926      $jline = $_; 
     927      # 
     928      # 2. Read the next line and check if it contains DO_2D 
     929      # 
     930      my $iline = <F> || die "DO jk line at end of file?"; 
     931      if ( $iline =~ /^\s*DO_2D\s*\(/i) { 
     932         my $isinavlid = 0; 
     933         # 
     934         # 3. Initialise a count to track any nested do loops 
     935         # 
     936         my $docount = 0; 
     937         # 
     938         # 4. Store the loop limits from the two lines stored and remove spaces and new-lines 
     939         # 
     940         ($jargs = $jline) =~ s/(^.*)=([^\!\n]*)(\!*.*)/\2/; 
     941         ($iargs = $iline) =~ s/(^.*)\(([^\!\n]*)\)(\!*.*)/\2/; 
     942         chomp($iargs); 
     943         chomp($jargs); 
     944         $iargs =~ s/^\s+//; $iargs =~ s/\s+$//; 
     945         $jargs =~ s/^\s+//; $jargs =~ s/\s+$//; 
     946         # 
     947         # 5. Store the leading indentation for the outer loop 
     948         # 
     949         ($jspac = $jline) =~ s/(^[\s]*)([^\s]*).*/\1/; 
     950         chomp($jspac); 
     951         # 
     952         # 6. Construct a DO_3D line to replace the original statements 
     953         # 
     954         # Keep two versions of output until we know if it is transformable 
     955         # 
     956         my $ostr = ""; 
     957         my $astr = ""; 
     958         my $orig = new IO::Scalar \$ostr; 
     959         print $orig $jline; 
     960         print $orig $iline; 
     961         my $altr = new IO::Scalar \$astr; 
     962         print $altr $jspac,"DO_3D( ",$iargs," , ",$jargs," )\n"; 
     963         # 
     964         # 7. Now process the loop contents until the matching  END_2D statement 
     965         # 
     966         while ( $docount >= 0 || ! ( $iline =~ /^\s*END_2D/i ) ) { 
     967            $iline = <F> || eval{ $isinvalid = 1 }; 
     968            #print $orig $iline; 
     969            # 
     970            # 8. Increment a counter if another DO_2D statement is found 
     971            # 
     972            if ( $iline =~ /^\s*do_2d/i )  { $docount++ }; 
     973            # 
     974            # 9. Decrement a counter if a END DO statement is found 
     975            # 
     976            if ( $iline =~ /^\s*end_2d/i )  { $docount-- }; 
     977            # 
     978            # 10. A negative counter means the matching END_2D for the ji loop has been reached 
     979            # 
     980            if ( $docount < 0 ) { 
     981               # 
     982               # 11. Check the next line is the expected END DO for the jk loop. 
     983               #     Output END_3D statement if it is 
     984               # 
     985               $jline = <F> || eval {$isinvalid = 1} ; 
     986               if ( ! ($jline =~ /^\s*end\s*do/i) )  { 
     987                  $isinvalid = 1 ; 
     988                  print $orig $iline; 
     989                  print $orig $jline; 
     990               } else { 
     991                  print $altr $jspac,"END_3D\n"; 
     992               } 
     993               if ( $isinvalid == 0 ) { 
     994                  print $altr; 
     995               } else { 
     996                  print $orig; 
     997                  $isinvalid = 0; 
     998               } 
     999 
     1000            } else { 
     1001               # 
     1002               # 12. This is a line inside the loop. Remove three leading spaces (if any) and output. 
     1003               # 
     1004               print $orig $iline; 
     1005               $iline =~ s/^\s\s\s//; 
     1006               print $altr $iline; 
     1007            } 
     1008         } 
     1009      } else { 
     1010         # 
     1011         # 13. Consecutive DO statements were not found. Do not process these loops. 
     1012         # 
     1013         print $jline; 
     1014         print $iline; 
     1015      } 
     1016   } else { 
     1017      # 
     1018      # 14. Code outside of a DO construct. Leave unchanged. 
     1019      # 
     1020      print $_; 
     1021   } 
     1022} 
     1023}}} 
     1024These scripts can be run sequentially; for example: 
     1025{{{ 
     1026  cp TEST_FILES_ORG/dynadv_cen2.F90 TESTDO_FILES/ 
     1027  perl do2dfinder.pl TESTDO_FILES/dynadv_cen2.F90 > TESTDO_FILES_2D/dynadv_cen2.F90 
     1028  perl do3dfinder.pl TESTDO_FILES_2D/dynadv_cen2.F90 > TESTDO_FILES_3D/dynadv_cen2.F90 
     1029}}} 
     1030And the difference between the original and final version is: 
     1031{{{#!diff 
     1032--- TESTDO_FILES/dynadv_cen2.F90        2019-03-01 13:35:29.000000000 +0000 
     1033+++ TESTDO_FILES_3D/dynadv_cen2.F90     2019-03-01 13:40:02.000000000 +0000 
     1034@@ -68,22 +68,18 @@ 
     1035       DO jk = 1, jpkm1                    ! horizontal transport 
     1036          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
     1037          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     1038-         DO jj = 1, jpjm1                 ! horizontal momentum fluxes (at T- and F-point) 
     1039-            DO ji = 1, fs_jpim1   ! vector opt. 
     1040-               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
     1041-               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
     1042-               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
     1043-               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
    8961044-            END DO 
    897 +            DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
    898 +               zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
    899 +                  &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    900 +            END_2D 
    901           CASE ( np_MET )                           !* metric term 
    902 -            DO jj = 1, jpjm1 
    903 -               DO ji = 1, fs_jpim1   ! vector opt. 
    904 -                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    905 -                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    906 -               END DO 
     1045-         END DO 
     1046-         DO jj = 2, jpjm1                 ! divergence of horizontal momentum fluxes 
     1047-            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1048-               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
     1049-                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     1050-               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
     1051-                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
    9071052-            END DO 
    908 +            DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
    909 +               zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    910 +                  &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    911 +            END_2D 
    912           CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    913 -            DO jj = 1, jpjm1 
    914 -               DO ji = 1, fs_jpim1   ! vector opt. 
    915 -                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
    916 -                     &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    917 -               END DO 
    918 -            END DO 
    919 +            DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
    920 +               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
    921 +                  &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    922 +            END_2D 
    923           CASE ( np_CME )                           !* Coriolis + metric 
    924 -            DO jj = 1, jpjm1 
    925 -               DO ji = 1, fs_jpim1   ! vector opt. 
    926 -                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    927 -                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    928 -               END DO 
    929 -            END DO 
    930 +            DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
    931 +               zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
    932 +                  &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 
    933 +            END_2D 
    934           CASE DEFAULT                                             ! error 
    935              CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    936           END SELECT 
    937           ! 
    938           IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    939 -            DO jj = 1, jpjm1 
    940 -               DO ji = 1, fs_jpim1   ! vector opt. 
    941 -                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    942 -               END DO 
    943 -            END DO 
    944 +            DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
    945 +               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    946 +            END_2D 
    947           ENDIF 
    948  
    949           IF( ln_sco ) THEN 
    950 @@ -413,16 +385,14 @@ 
    951              zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
    952           ENDIF 
    953           !                                   !==  compute and add the vorticity term trend  =! 
     1053-         END DO 
     1054+         DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
     1055+            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
     1056+            zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
     1057+            zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
     1058+            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     1059+         END_2D 
     1060+         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     1061+            ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
     1062+               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     1063+            va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
     1064+               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     1065+         END_2D 
     1066       END DO 
     1067       ! 
     1068       IF( l_trddyn ) THEN           ! trends: send trend to trddyn for diagnostic 
     1069@@ -96,41 +92,29 @@ 
     1070       ! 
     1071       !                             !==  Vertical advection  ==! 
     1072       ! 
     1073-      DO jj = 2, jpjm1                    ! surface/bottom advective fluxes set to zero 
     1074-         DO ji = fs_2, fs_jpim1 
     1075-            zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp 
     1076-            zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp 
     1077-         END DO 
     1078-      END DO 
     1079+      DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     1080+         zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp 
     1081+         zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp 
     1082+      END_2D 
     1083       IF( ln_linssh ) THEN                ! linear free surface: advection through the surface 
    9541084-         DO jj = 2, jpjm1 
    955 -            DO ji = fs_2, fs_jpim1   ! vector opt. 
    956 -               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
    957 -               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
    958 -               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    959 -               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    960 -               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    961 -               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
     1085-            DO ji = fs_2, fs_jpim1 
     1086-               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 
     1087-               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 
    9621088-            END DO 
    9631089-         END DO 
    9641090+         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
    965 +            zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
    966 +            zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
    967 +            zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    968 +            zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    969 +            pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    970 +            pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
     1091+            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 
     1092+            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 
    9711093+         END_2D 
    972           !                                             ! =============== 
    973        END DO                                           !   End of slab 
    974        !                                                ! =============== 
    975 }}} 
    976  
     1094       ENDIF 
     1095       DO jk = 2, jpkm1                    ! interior advective fluxes 
     1096-         DO jj = 2, jpj                       ! 1/4 * Vertical transport 
     1097-            DO ji = 2, jpi 
     1098-               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     1099-            END DO 
     1100-         END DO 
     1101-         DO jj = 2, jpjm1 
     1102-            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1103-               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
     1104-               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     1105-            END DO 
     1106-         END DO 
     1107-      END DO 
     1108-      DO jk = 1, jpkm1                    ! divergence of vertical momentum flux divergence 
     1109-         DO jj = 2, jpjm1 
     1110-            DO ji = fs_2, fs_jpim1   ! vector opt. 
     1111-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     1112-               va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     1113-            END DO 
     1114-         END DO 
     1115+         DO_2D( 2, jpi , 2, jpj ) 
     1116+            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     1117+         END_2D 
     1118+         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     1119+            zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
     1120+            zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     1121+         END_2D 
     1122       END DO 
     1123+      DO_3D( fs_2, fs_jpim1 , 2, jpjm1 , 1, jpkm1 ) 
     1124+         ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     1125+         va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     1126+      END_3D 
     1127       ! 
     1128       IF( l_trddyn ) THEN                 ! trends: send trend to trddyn for diagnostic 
     1129          zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
     1130}}} 
    9771131 
    9781132