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


Ignore:
Timestamp:
2019-03-01T14:59:39+01:00 (5 years ago)
Author:
acc
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps

    v10 v11  
    553553       END DO 
    554554       ! 
    555 }} 
     555}}} 
    556556 
    557557 
    … …  
    562562 
    563563So far so good.... 
     564 
     565'''Automating the tiling changes''' 
     566 
     567 
     568Here 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: 
     569{{{ 
     570# 
     571open(F,$ARGV[0]) || die "Cannot open $ARGV[0]: $!"; 
     572while(<F>) { 
     573   if ( $_ =~ /^\s*DO\s* jj/i) { 
     574      # Start processing loop if line contains DO jj (case and whitespace independent) 
     575      # 
     576      # 1. Store the current line 
     577      # 
     578      $jline = $_; 
     579      # 
     580      # 2. Read the next line and check if it contains DO ji 
     581      # 
     582      my $iline = <F> || die "DO jj line at end of file?"; 
     583      if ( $iline =~ /^\s*do\s*ji/i) { 
     584         # 
     585         # 3. Initialise a count to track any nested do loops 
     586         # 
     587         my $docount = 0; 
     588         # 
     589         # 4. Store the loop limits from the two lines stored and remove spaces and new-lines 
     590         # 
     591         ($jargs = $jline) =~ s/(^.*)=([^\!\n]*)(\!*.*)/\2/; 
     592         ($iargs = $iline) =~ s/(^.*)=([^\!\n]*)(\!*.*)/\2/; 
     593         chomp($iargs); 
     594         chomp($jargs); 
     595         $iargs =~ s/^\s+//; $iargs =~ s/\s+$//; 
     596         $jargs =~ s/^\s+//; $jargs =~ s/\s+$//; 
     597         # 
     598         # 5. Store the leading indentation for the outer loop 
     599         # 
     600         ($jspac = $jline) =~ s/(^[\s]*)([^\s]*).*/\1/; 
     601         chomp($jspac); 
     602         # 
     603         # 6. Construct a DO_2D line to replace the original statements 
     604         # 
     605         print $jspac,"DO_2D( ",$iargs," , ",$jargs," )\n"; 
     606         # 
     607         # 7. Now process the loop contents until the matching pair of END DO statements 
     608         # 
     609         while ( $docount >= 0 || ! ( $iline =~ /^\s*end\s*do/i ) ) { 
     610            $iline = <F> || die "reached end of file within do loop?" ; 
     611            # 
     612            # 8. Increment a counter if another DO statement is found 
     613            # 
     614            if ( $iline =~ /^\s*do\s*/i )  { $docount++ }; 
     615            # 
     616            # 9. Decrement a counter if a END DO statement is found 
     617            # 
     618            if ( $iline =~ /^\s*end\s*do/i )  { $docount-- }; 
     619            # 
     620            # 10. A negative counter means the matching END DO for the ji loop has been reached 
     621            # 
     622            if ( $docount < 0 ) { 
     623               # 
     624               # 11. Check the next line is the expected 2nd END DO. 
     625               #     Output END_2D statement if it is 
     626               # 
     627               $jline = <F> || die "reached end of file looking for end do?" ; 
     628               if ( ! ($jline =~ /^\s*end\s*do/i) )  { 
     629                  die "END DOs are not consecutive"; 
     630               } else { 
     631                  print $jspac,"END_2D\n"; 
     632               } 
     633 
     634            } else { 
     635               # 
     636               # 12. This is a line inside the loop. Remove three leading spaces (if any) and output. 
     637               # 
     638               $iline =~ s/^\s\s\s//; 
     639               print $iline; 
     640            } 
     641         } 
     642      } else { 
     643         # 
     644         # 13. Consecutive DO statements were not found. Do not process these loops. 
     645         # 
     646         print $jline; 
     647         print $iline; 
     648      } 
     649   } else { 
     650      # 
     651      # 14. Code outside of a DO construct. Leave unchanged. 
     652      # 
     653      print $_; 
     654   } 
     655} 
     656}}} 
     657Secondly, here is a simple test file: 
     658{{{ 
     659!   some random text 
     660!   followed by a valid loop pair 
     661     DO jj  = 2,jpjm1     ! with comments 
     662        DO ji = 2,jpim1 
     663           some loop content 
     664           more loop content 
     665        END DO 
     666     END DO 
     667!   followed by an invalid loop pair 
     668     DO jj  = 2,jpjm1 
     669        j = jj-1 
     670 
     671        DO ji = 2,jpim1 
     672           some loop content 
     673           more loop content 
     674        END DO 
     675     END DO 
     676!   followed by an valid loop pair with a nested do 
     677     DO jj  = 2,jpjm1 
     678        DO ji = 2,jpim1 
     679           some loop content 
     680           do jn = 1, jptrc 
     681              yet more loop content 
     682           end do 
     683           more loop content 
     684        END DO 
     685     END DO 
     686}}} 
     687Followed by the result of running: `perl do2dfinder.pl testdo.F90 > testdo_after.F90` and the difference: 
     688{{{ 
     689!   some random text 
     690!   followed by a valid loop pair 
     691     DO_2D( 2,jpim1 , 2,jpjm1 ) 
     692        some loop content 
     693        more loop content 
     694     END_2D 
     695!   followed by an invalid loop pair 
     696     DO jj  = 2,jpjm1 
     697        j = jj-1 
     698 
     699        DO ji = 2,jpim1 
     700           some loop content 
     701           more loop content 
     702        END DO 
     703     END DO 
     704!   followed by an valid loop pair with a nested do 
     705     DO_2D( 2,jpim1 , 2,jpjm1 ) 
     706        some loop content 
     707        do jn = 1, jptrc 
     708           yet more loop content 
     709        end do 
     710        more loop content 
     711     END_2D 
     712}}} 
     713{{{#!diff 
     714--- testdo.F90  2019-02-28 15:52:33.000000000 +0000 
     715+++ testdo_after.F90    2019-02-28 15:53:44.000000000 +0000 
     716@@ -1,11 +1,9 @@ 
     717 !   some random text 
     718 !   followed by a valid loop pair 
     719-     DO jj  = 2,jpjm1     ! with comments 
     720-        DO ji = 2,jpim1 
     721-           some loop content 
     722-           more loop content 
     723-        END DO 
     724-     END DO 
     725+     DO_2D( 2,jpim1 , 2,jpjm1 ) 
     726+        some loop content 
     727+        more loop content 
     728+     END_2D 
     729 !   followed by an invalid loop pair 
     730      DO jj  = 2,jpjm1 
     731         j = jj-1 
     732@@ -16,12 +14,10 @@ 
     733         END DO 
     734      END DO 
     735 !   followed by an valid loop pair with a nested do 
     736-     DO jj  = 2,jpjm1 
     737-        DO ji = 2,jpim1 
     738-           some loop content 
     739-           do jn = 1, jptrc 
     740-              yet more loop content 
     741-           end do 
     742-           more loop content 
     743-        END DO 
     744-     END DO 
     745+     DO_2D( 2,jpim1 , 2,jpjm1 ) 
     746+        some loop content 
     747+        do jn = 1, jptrc 
     748+           yet more loop content 
     749+        end do 
     750+        more loop content 
     751+     END_2D 
     752}}} 
     753 
     754And finally, the results on a real case dynvor.F90 (just a sample of the differences): 
     755{{{#!diff 
     756--- dynvor.F90  2019-02-28 16:01:28.000000000 +0000 
     757+++ dynvor_after.F90    2019-02-28 16:01:55.000000000 +0000 
     758@@ -225,18 +225,14 @@ 
     759       SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
     760       CASE ( np_RVO )                           !* relative vorticity 
     761          DO jk = 1, jpkm1                                 ! Horizontal slab 
     762-            DO jj = 1, jpjm1 
     763-               DO ji = 1, jpim1 
     764-                  zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
     765-                     &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     766-               END DO 
     767-            END DO 
     768+            DO_2D( 1, jpim1 , 1, jpjm1 ) 
     769+               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
     770+                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     771+            END_2D 
     772             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
     773-               DO jj = 1, jpjm1 
     774-                  DO ji = 1, jpim1 
     775-                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     776-                  END DO 
     777-               END DO 
     778+               DO_2D( 1, jpim1 , 1, jpjm1 ) 
     779+                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     780+               END_2D 
     781             ENDIF 
     782          END DO 
     783 
     784@@ -244,18 +240,14 @@ 
     785 
     786       CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     787          DO jk = 1, jpkm1                                 ! Horizontal slab 
     788-            DO jj = 1, jpjm1 
     789-               DO ji = 1, jpim1                          ! relative vorticity 
     790-                  zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
     791-                     &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
     792-               END DO 
     793-            END DO 
     794+            DO_2D( 1, jpim1 , 1, jpjm1 ) 
     795+               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
     796+                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
     797+            END_2D 
     798             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
     799-               DO jj = 1, jpjm1 
     800-                  DO ji = 1, jpim1 
     801-                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     802-                  END DO 
     803-               END DO 
     804+               DO_2D( 1, jpim1 , 1, jpjm1 ) 
     805+                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     806+               END_2D 
     807             ENDIF 
     808          END DO 
     809 
     810@@ -271,50 +263,40 @@ 
     811          CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     812             zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 
     813          CASE ( np_RVO )                           !* relative vorticity 
     814-            DO jj = 2, jpj 
     815-               DO ji = 2, jpi   ! vector opt. 
     816-                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
     817-                     &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     818-               END DO 
     819-            END DO 
     820+            DO_2D( 2, jpi , 2, jpj ) 
     821+               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
     822+                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     823+            END_2D 
     824          CASE ( np_MET )                           !* metric term 
     825-            DO jj = 2, jpj 
     826-               DO ji = 2, jpi 
     827-                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
     828-                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
     829-               END DO 
     830-            END DO 
     831+            DO_2D( 2, jpi , 2, jpj ) 
     832+               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
     833+                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
     834+            END_2D 
     835          CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     836-            DO jj = 2, jpj 
     837-               DO ji = 2, jpi   ! vector opt. 
     838-                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
     839-                     &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     840-               END DO 
     841-            END DO 
     842+            DO_2D( 2, jpi , 2, jpj ) 
     843+               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
     844+                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     845+            END_2D 
     846          CASE ( np_CME )                           !* Coriolis + metric 
     847-            DO jj = 2, jpj 
     848-               DO ji = 2, jpi   ! vector opt. 
     849-                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
     850-                       &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
     851-                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
     852-               END DO 
     853-            END DO 
     854+            DO_2D( 2, jpi , 2, jpj ) 
     855+               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
     856+                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
     857+                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
     858+            END_2D 
     859          CASE DEFAULT                                             ! error 
     860             CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
     861          END SELECT 
     862          ! 
     863          !                                   !==  compute and add the vorticity term trend  =! 
     864-         DO jj = 2, jpjm1 
     865-            DO ji = 2, jpim1   ! vector opt. 
     866-               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
     867-                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
     868-                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
     869-                  ! 
     870-               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
     871-                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
     872-                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
     873-            END DO 
     874-         END DO 
     875+         DO_2D( 2, jpim1 , 2, jpjm1 ) 
     876+            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
     877+               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
     878+               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
     879+               ! 
     880+            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
     881+               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
     882+               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
     883+         END_2D 
     884          !                                             ! =============== 
     885       END DO                                           !   End of slab 
     886       !                                                ! =============== 
     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 
     896-            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 
     907-            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  =! 
     954-         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 ) 
     962-            END DO 
     963-         END DO 
     964+         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 ) 
     971+         END_2D 
     972          !                                             ! =============== 
     973       END DO                                           !   End of slab 
     974       !                                                ! =============== 
     975}}} 
     976 
    564977 
    565978