Changes between Version 12 and Version 13 of 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps


Ignore:
Timestamp:
2019-03-04T15:19:10+01:00 (17 months ago)
Author:
acc
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • 2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps

    v12 v13  
    565565== Automating the tiling changes 
    566566 
    567  
    568 Here 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: 
     567Here 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: 
     568{{{ 
     569cat TESTDO_FILES/testdo.F90 
     570!   some random text 
     571!   followed by a valid loop pair 
     572     DO jj  = 2,jpjm1     ! with comments 
     573        DO ji = 2,jpim1 
     574           some loop content 
     575           more loop content 
     576        END DO 
     577     END DO 
     578!   followed by an invalid loop pair 
     579     DO jj  = 2,jpjm1 
     580        j = jj-1 
     581 
     582        DO ji = 2,jpim1 
     583           some loop content 
     584           more loop content 
     585        END DO 
     586     END DO 
     587!   followed by a valid loop pair with a nested do 
     588     DO jj  = 1,jpjm1 
     589        DO ji = 1,jpi 
     590           some loop content 
     591           do jn = 1, jptrc 
     592              yet more loop content 
     593           end do 
     594           more loop content 
     595        END DO 
     596     END DO 
     597!   followed by a valid 3D loop 
     598     DO jk  = 2,jpkm1 
     599        DO jj  = 1,jpjm1 
     600           DO ji = 1,jpi 
     601              some loop content 
     602              more loop content 
     603           END DO 
     604        END DO 
     605     END DO 
     606}}} 
     607This file is processed as follows: 
     608{{{ 
     609  perl do2dfinder.pl TESTDO_FILES/testdo.F90 > TESTDO_FILES_2D/testdo.F90 
     610  perl do3dfinder.pl TESTDO_FILES_2D/testdo.F90 > TESTDO_FILES_3D/testdo.F90 
     611  diff -u TESTDO_FILES/testdo.F90 TESTDO_FILES_3D/testdo.F90 > testdo.patch 
     612}}} 
     613and the resulting changes are: 
     614{{{#!diff 
     615--- TESTDO_FILES/testdo.F90     2019-03-04 13:43:28.000000000 +0000 
     616+++ TESTDO_FILES_3D/testdo.F90  2019-03-04 13:43:38.000000000 +0000 
     617@@ -1,11 +1,9 @@ 
     618 !   some random text 
     619 !   followed by a valid loop pair 
     620-     DO jj  = 2,jpjm1     ! with comments 
     621-        DO ji = 2,jpim1 
     622-           some loop content 
     623-           more loop content 
     624-        END DO 
     625-     END DO 
     626+     DO_2D_00_00 
     627+        some loop content 
     628+        more loop content 
     629+     END_2D 
     630 !   followed by an invalid loop pair 
     631      DO jj  = 2,jpjm1 
     632         j = jj-1 
     633@@ -16,21 +14,15 @@ 
     634         END DO 
     635      END DO 
     636 !   followed by a valid loop pair with a nested do 
     637-     DO jj  = 1,jpjm1 
     638-        DO ji = 1,jpi 
     639-           some loop content 
     640-           do jn = 1, jptrc 
     641-              yet more loop content 
     642-           end do 
     643-           more loop content 
     644-        END DO 
     645-     END DO 
     646+     DO_2D_10_11 
     647+        some loop content 
     648+        do jn = 1, jptrc 
     649+           yet more loop content 
     650+        end do 
     651+        more loop content 
     652+     END_2D 
     653 !   followed by a valid 3D loop 
     654-     DO jk  = 2,jpkm1 
     655-        DO jj  = 1,jpjm1 
     656-           DO ji = 1,jpi 
     657-              some loop content 
     658-              more loop content 
     659-           END DO 
     660-        END DO 
     661-     END DO 
     662+     DO_3D_10_11( 2,jpkm1 ) 
     663+        some loop content 
     664+        more loop content 
     665+     END_3D 
     666}}} 
     667traldf_iso.F90 provides a more strigent test: 
     668{{{ 
     669   perl do2dfinder.pl TESTDO_FILES/traldf_iso.F90 > TESTDO_FILES_2D/traldf_iso.F90 
     670   perl do3dfinder.pl TESTDO_FILES_2D/traldf_iso.F90 > TESTDO_FILES_3D/traldf_iso.F90 
     671   diff -u TESTDO_FILES/traldf_iso.F90 TESTDO_FILES_3D/traldf_iso.F90 > traldf_iso.patch 
     672}}} 
     673{{{#!diff 
     674--- TESTDO_FILES/traldf_iso.F90 2019-03-04 12:59:44.000000000 +0000 
     675+++ TESTDO_FILES_3D/traldf_iso.F90      2019-03-04 13:26:08.000000000 +0000 
     676@@ -143,58 +143,42 @@ 
     677       ! 
     678       IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
     679          ! 
     680-         DO jk = 2, jpkm1 
     681-            DO jj = 2, jpjm1 
     682-               DO ji = fs_2, fs_jpim1   ! vector opt. 
     683-                  ! 
     684-                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     685-                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     686-                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     687-                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     688-                     ! 
     689-                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     690-                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     691-                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     692-                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     693-                     ! 
     694-                  ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     695-                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     696-               END DO 
     697-            END DO 
     698-         END DO 
     699+         DO_3D_00_00( 2, jpkm1 ) 
     700+            ! 
     701+            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     702+               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     703+            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     704+               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     705+               ! 
     706+            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     707+               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     708+            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     709+               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     710+               ! 
     711+            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     712+               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     713+         END_3D 
     714          ! 
     715          IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
     716-            DO jk = 2, jpkm1 
     717-               DO jj = 2, jpjm1 
     718-                  DO ji = fs_2, fs_jpim1 
     719-                     akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     720-                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     721-                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     722-                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     723-                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     724-                  END DO 
     725-               END DO 
     726-            END DO 
     727+            DO_3D_00_00( 2, jpkm1 ) 
     728+               akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     729+                  &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     730+                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     731+                  &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     732+                  &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     733+            END_3D 
     734             ! 
     735             IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
     736-               DO jk = 2, jpkm1 
     737-                  DO jj = 1, jpjm1 
     738-                     DO ji = 1, fs_jpim1 
     739-                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     740-                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     741-                     END DO 
     742-                  END DO 
     743-               END DO 
     744+               DO_3D_10_10( 2, jpkm1 ) 
     745+                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     746+                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     747+               END_3D 
     748             ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
     749-               DO jk = 2, jpkm1 
     750-                  DO jj = 1, jpjm1 
     751-                     DO ji = 1, fs_jpim1 
     752-                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     753-                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     754-                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     755-                     END DO 
     756-                  END DO 
     757-               END DO 
     758+               DO_3D_10_10( 2, jpkm1 ) 
     759+                  ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     760+                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     761+                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     762+               END_3D 
     763            ENDIF 
     764            ! 
     765          ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
     766@@ -215,28 +199,20 @@ 
     767          !!end 
     768 
     769          ! Horizontal tracer gradient 
     770-         DO jk = 1, jpkm1 
     771-            DO jj = 1, jpjm1 
     772-               DO ji = 1, fs_jpim1   ! vector opt. 
     773-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     774-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     775-               END DO 
     776-            END DO 
     777-         END DO 
     778+         DO_3D_10_10( 1, jpkm1 ) 
     779+            zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     780+            zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     781+         END_3D 
     782          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
     783-            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell) 
     784-               DO ji = 1, fs_jpim1   ! vector opt. 
     785-                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
     786-                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     787-               END DO 
     788-            END DO 
     789+            DO_2D_10_10 
     790+               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
     791+               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     792+            END_2D 
     793             IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
     794-               DO jj = 1, jpjm1 
     795-                  DO ji = 1, fs_jpim1   ! vector opt. 
     796-                     IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
     797-                     IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
     798-                  END DO 
     799-               END DO 
     800+               DO_2D_10_10 
     801+                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
     802+                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
     803+               END_2D 
     804             ENDIF 
     805          ENDIF 
     806          ! 
     807@@ -252,36 +228,32 @@ 
     808             IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
     809             ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
     810             ENDIF 
     811-            DO jj = 1 , jpjm1            !==  Horizontal fluxes 
     812-               DO ji = 1, fs_jpim1   ! vector opt. 
     813-                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
     814-                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     815-                  ! 
     816-                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     817-                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     818-                  ! 
     819-                  zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     820-                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     821-                  ! 
     822-                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     823-                  zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     824-                  ! 
     825-                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     826-                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     827-                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     828-                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     829-                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     830-                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
     831-               END DO 
     832-            END DO 
     833-            ! 
     834-            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
     835-               DO ji = fs_2, fs_jpim1   ! vector opt. 
     836-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     837-                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     838-                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     839-               END DO 
     840-            END DO 
     841+            DO_2D_10_10 
     842+               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
     843+               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     844+               ! 
     845+               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     846+                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     847+               ! 
     848+               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     849+                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     850+               ! 
     851+               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     852+               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     853+               ! 
     854+               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     855+                  &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     856+                  &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     857+               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     858+                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     859+                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
     860+            END_2D 
     861+            ! 
     862+            DO_2D_00_00 
     863+               pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     864+                  &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     865+                  &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     866+            END_2D 
     867          END DO                                        !   End of slab 
     868 
     869          !!---------------------------------------------------------------------- 
     870@@ -295,75 +267,55 @@ 
     871          !                          ! Surface and bottom vertical fluxes set to zero 
     872          ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
     873 
     874-         DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1) 
     875-            DO jj = 2, jpjm1 
     876-               DO ji = fs_2, fs_jpim1   ! vector opt. 
     877-                  ! 
     878-                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     879-                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     880-                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     881-                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     882-                     ! 
     883-                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     884-                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     885-                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     886-                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     887-                     ! 
     888-                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     889-                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     890-                  ! 
     891-                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     892-                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     893-                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     894-                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     895-               END DO 
     896-            END DO 
     897-         END DO 
     898+         DO_3D_00_00( 2, jpkm1 ) 
     899+            ! 
     900+            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     901+               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     902+            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     903+               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     904+               ! 
     905+            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     906+               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     907+            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     908+               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     909+               ! 
     910+            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     911+            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     912+            ! 
     913+            ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     914+               &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     915+               &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     916+               &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     917+         END_3D 
     918          !                                !==  add the vertical 33 flux  ==! 
     919          IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
     920-            DO jk = 2, jpkm1 
     921-               DO jj = 1, jpjm1 
     922-                  DO ji = fs_2, fs_jpim1 
     923-                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     924-                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
     925-                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     926-                  END DO 
     927-               END DO 
     928-            END DO 
     929+            DO_3D_10_00( 2, jpkm1 ) 
     930+               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     931+                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
     932+                  &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     933+            END_3D 
     934             ! 
     935          ELSE                                   ! bilaplacian 
     936             SELECT CASE( kpass ) 
     937             CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
     938-               DO jk = 2, jpkm1 
     939-                  DO jj = 1, jpjm1 
     940-                     DO ji = fs_2, fs_jpim1 
     941-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
     942-                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
     943-                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     944-                     END DO 
     945-                  END DO 
     946-               END DO 
     947+               DO_3D_10_00( 2, jpkm1 ) 
     948+                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
     949+                     &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
     950+                     &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     951+               END_3D 
     952             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     953-               DO jk = 2, jpkm1 
     954-                  DO jj = 1, jpjm1 
     955-                     DO ji = fs_2, fs_jpim1 
     956-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
     957-                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
     958-                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     959-                     END DO 
     960-                  END DO 
     961-               END DO 
     962+               DO_3D_10_00( 2, jpkm1 ) 
     963+                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
     964+                     &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
     965+                     &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     966+               END_3D 
     967             END SELECT 
     968          ENDIF 
     969          ! 
     970-         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
     971-            DO jj = 2, jpjm1 
     972-               DO ji = fs_2, fs_jpim1   ! vector opt. 
     973-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     974-                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     975-               END DO 
     976-            END DO 
     977-         END DO 
     978+         DO_3D_00_00( 1, jpkm1 ) 
     979+            pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     980+               &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     981+         END_3D 
     982          ! 
     983          IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
     984              ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==! 
     985}}} 
     986And finally the two scripts that achieve this: 
    569987{{{ 
    570988cat do2dfinder.pl 
     
    5901008         # 4. Store the loop limits from the two lines stored and remove spaces and new-lines 
    5911009         # 
    592          ($jargs = $jline) =~ s/(^.*)=([^\!\n]*)(\!*.*)/\2/; 
    593          ($iargs = $iline) =~ s/(^.*)=([^\!\n]*)(\!*.*)/\2/; 
     1010         ($jargs = $jline) =~ s/(^.*?)=([^\!\n]*)(\!*.*)/\2/; 
     1011         ($iargs = $iline) =~ s/(^.*?)=([^\!\n]*)(\!*.*)/\2/; 
    5941012         chomp($iargs); 
    5951013         chomp($jargs); 
    5961014         $iargs =~ s/^\s+//; $iargs =~ s/\s+$//; 
    5971015         $jargs =~ s/^\s+//; $jargs =~ s/\s+$//; 
    598          # 
    599          # 5. Store the leading indentation for the outer loop 
     1016         @ilims = split(/,/,$iargs); 
     1017         @jlims = split(/,/,$jargs); 
     1018         # 
     1019         # 5. Convert i-limits into Gurvan's indexing 
     1020         # i.e. lower limit: [2,fs_2] --> 0 ; 1--> 1 
     1021         # i.e. upper limit: *m* --> 0 else --> 1 
     1022         # 
     1023         if( @ilims[0] =~ /2/ ) { 
     1024           @ilims[0] = 0; 
     1025         } elsif ( @ilims[0] == 1 ) { 
     1026           @ilims[0] = 1; 
     1027         } 
     1028         if( @ilims[1] =~ /m/ ) { 
     1029           @ilims[1] = 0; 
     1030         } else { 
     1031           @ilims[1] = 1; 
     1032         } 
     1033         $iargs = join('',@ilims[0], @ilims[1]); 
     1034         # 
     1035         # 6. Convert j-limits into Gurvan's indexing 
     1036         # i.e. lower limit: [2,fs_2] --> 0 ; 1--> 1 
     1037         # i.e. upper limit: *m* --> 0 else --> 1 
     1038         # 
     1039         if( @jlims[0] =~ /2/ ) { 
     1040           @jlims[0] = 0; 
     1041         } elsif ( @jlims[0] == 1 ) { 
     1042           @jlims[0] = 1; 
     1043         } 
     1044         if( @jlims[1] =~ /m/ ) { 
     1045           @jlims[1] = 0; 
     1046         } else { 
     1047           @jlims[1] = 1; 
     1048         } 
     1049         $jargs = join('',@jlims[0], @jlims[1]); 
     1050         # 
     1051         # 7. Store the leading indentation for the outer loop 
    6001052         # 
    6011053         ($jspac = $jline) =~ s/(^[\s]*)([^\s]*).*/\1/; 
    6021054         chomp($jspac); 
    6031055         # 
    604          # 6. Construct a DO_2D line to replace the original statements 
    605          # 
    606          print $jspac,"DO_2D( ",$iargs," , ",$jargs," )\n"; 
    607          # 
    608          # 7. Now process the loop contents until the matching pair of END DO statements 
     1056         # 8. Construct a DO_2D line to replace the original statements 
     1057         # 
     1058         print $jspac,join('_',"DO_2D",$jargs,$iargs),"\n"; 
     1059         # 
     1060         # 9. Now process the loop contents until the matching pair of END DO statements 
    6091061         # 
    6101062         while ( $docount >= 0 || ! ( $iline =~ /^\s*end\s*do/i ) ) { 
    6111063            $iline = <F> || die "reached end of file within do loop?" ; 
    6121064            # 
    613             # 8. Increment a counter if another DO statement is found 
     1065            # 10. Increment a counter if another DO statement is found 
    6141066            # 
    6151067            if ( $iline =~ /^\s*do\s*/i )  { $docount++ }; 
    6161068            # 
    617             # 9. Decrement a counter if a END DO statement is found 
     1069            # 11. Decrement a counter if a END DO statement is found 
    6181070            # 
    6191071            if ( $iline =~ /^\s*end\s*do/i )  { $docount-- }; 
    6201072            # 
    621             # 10. A negative counter means the matching END DO for the ji loop has been reached 
     1073            # 12. A negative counter means the matching END DO for the ji loop has been reached 
    6221074            # 
    6231075            if ( $docount < 0 ) { 
    6241076               # 
    625                # 11. Check the next line is the expected 2nd END DO. 
     1077               # 13. Check the next line is the expected 2nd END DO. 
    6261078               #     Output END_2D statement if it is 
    6271079               # 
     
    6351087            } else { 
    6361088               # 
    637                # 12. This is a line inside the loop. Remove three leading spaces (if any) and output. 
     1089               # 14. This is a line inside the loop. Remove three leading spaces (if any) and output. 
    6381090               # 
    6391091               $iline =~ s/^\s\s\s//; 
     
    6431095      } else { 
    6441096         # 
    645          # 13. Consecutive DO statements were not found. Do not process these loops. 
     1097         # 15. Consecutive DO statements were not found. Do not process these loops. 
    6461098         # 
    6471099         print $jline; 
     
    6501102   } else { 
    6511103      # 
    652       # 14. Code outside of a DO construct. Leave unchanged. 
     1104      # 16. Code outside of a DO construct. Leave unchanged. 
    6531105      # 
    6541106      print $_; 
     
    6561108} 
    6571109}}} 
    658 Secondly, here is a simple test file: 
    659 {{{ 
    660 !   some random text 
    661 !   followed by a valid loop pair 
    662      DO jj  = 2,jpjm1     ! with comments 
    663         DO ji = 2,jpim1 
    664            some loop content 
    665            more loop content 
    666         END DO 
    667      END DO 
    668 !   followed by an invalid loop pair 
    669      DO jj  = 2,jpjm1 
    670         j = jj-1 
    671  
    672         DO ji = 2,jpim1 
    673            some loop content 
    674            more loop content 
    675         END DO 
    676      END DO 
    677 !   followed by an valid loop pair with a nested do 
    678      DO jj  = 2,jpjm1 
    679         DO ji = 2,jpim1 
    680            some loop content 
    681            do jn = 1, jptrc 
    682               yet more loop content 
    683            end do 
    684            more loop content 
    685         END DO 
    686      END DO 
    687 }}} 
    688 Followed by the result of running: `perl do2dfinder.pl testdo.F90 > testdo_after.F90` and the difference: 
    689 {{{ 
    690 !   some random text 
    691 !   followed by a valid loop pair 
    692      DO_2D( 2,jpim1 , 2,jpjm1 ) 
    693         some loop content 
    694         more loop content 
    695      END_2D 
    696 !   followed by an invalid loop pair 
    697      DO jj  = 2,jpjm1 
    698         j = jj-1 
    699  
    700         DO ji = 2,jpim1 
    701            some loop content 
    702            more loop content 
    703         END DO 
    704      END DO 
    705 !   followed by an valid loop pair with a nested do 
    706      DO_2D( 2,jpim1 , 2,jpjm1 ) 
    707         some loop content 
    708         do jn = 1, jptrc 
    709            yet more loop content 
    710         end do 
    711         more loop content 
    712      END_2D 
    713 }}} 
    714 {{{#!diff 
    715 --- testdo.F90  2019-02-28 15:52:33.000000000 +0000 
    716 +++ testdo_after.F90    2019-02-28 15:53:44.000000000 +0000 
    717 @@ -1,11 +1,9 @@ 
    718  !   some random text 
    719  !   followed by a valid loop pair 
    720 -     DO jj  = 2,jpjm1     ! with comments 
    721 -        DO ji = 2,jpim1 
    722 -           some loop content 
    723 -           more loop content 
    724 -        END DO 
    725 -     END DO 
    726 +     DO_2D( 2,jpim1 , 2,jpjm1 ) 
    727 +        some loop content 
    728 +        more loop content 
    729 +     END_2D 
    730  !   followed by an invalid loop pair 
    731       DO jj  = 2,jpjm1 
    732          j = jj-1 
    733 @@ -16,12 +14,10 @@ 
    734          END DO 
    735       END DO 
    736  !   followed by an valid loop pair with a nested do 
    737 -     DO jj  = 2,jpjm1 
    738 -        DO ji = 2,jpim1 
    739 -           some loop content 
    740 -           do jn = 1, jptrc 
    741 -              yet more loop content 
    742 -           end do 
    743 -           more loop content 
    744 -        END DO 
    745 -     END DO 
    746 +     DO_2D( 2,jpim1 , 2,jpjm1 ) 
    747 +        some loop content 
    748 +        do jn = 1, jptrc 
    749 +           yet more loop content 
    750 +        end do 
    751 +        more loop content 
    752 +     END_2D 
    753 }}} 
    754  
    755 And finally, the results on a real case dynvor.F90 (just a sample of the differences): 
    756 {{{#!diff 
    757 --- dynvor.F90  2019-02-28 16:01:28.000000000 +0000 
    758 +++ dynvor_after.F90    2019-02-28 16:01:55.000000000 +0000 
    759 @@ -225,18 +225,14 @@ 
    760        SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
    761        CASE ( np_RVO )                           !* relative vorticity 
    762           DO jk = 1, jpkm1                                 ! Horizontal slab 
    763 -            DO jj = 1, jpjm1 
    764 -               DO ji = 1, jpim1 
    765 -                  zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
    766 -                     &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    767 -               END DO 
    768 -            END DO 
    769 +            DO_2D( 1, jpim1 , 1, jpjm1 ) 
    770 +               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
    771 +                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    772 +            END_2D 
    773              IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    774 -               DO jj = 1, jpjm1 
    775 -                  DO ji = 1, jpim1 
    776 -                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    777 -                  END DO 
    778 -               END DO 
    779 +               DO_2D( 1, jpim1 , 1, jpjm1 ) 
    780 +                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    781 +               END_2D 
    782              ENDIF 
    783           END DO 
    784  
    785 @@ -244,18 +240,14 @@ 
    786  
    787        CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    788           DO jk = 1, jpkm1                                 ! Horizontal slab 
    789 -            DO jj = 1, jpjm1 
    790 -               DO ji = 1, jpim1                          ! relative vorticity 
    791 -                  zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
    792 -                     &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    793 -               END DO 
    794 -            END DO 
    795 +            DO_2D( 1, jpim1 , 1, jpjm1 ) 
    796 +               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
    797 +                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    798 +            END_2D 
    799              IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    800 -               DO jj = 1, jpjm1 
    801 -                  DO ji = 1, jpim1 
    802 -                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    803 -                  END DO 
    804 -               END DO 
    805 +               DO_2D( 1, jpim1 , 1, jpjm1 ) 
    806 +                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    807 +               END_2D 
    808              ENDIF 
    809           END DO 
    810  
    811 @@ -271,50 +263,40 @@ 
    812           CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    813              zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 
    814           CASE ( np_RVO )                           !* relative vorticity 
    815 -            DO jj = 2, jpj 
    816 -               DO ji = 2, jpi   ! vector opt. 
    817 -                  zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    818 -                     &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
    819 -               END DO 
    820 -            END DO 
    821 +            DO_2D( 2, jpi , 2, jpj ) 
    822 +               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    823 +                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
    824 +            END_2D 
    825           CASE ( np_MET )                           !* metric term 
    826 -            DO jj = 2, jpj 
    827 -               DO ji = 2, jpi 
    828 -                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    829 -                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
    830 -               END DO 
    831 -            END DO 
    832 +            DO_2D( 2, jpi , 2, jpj ) 
    833 +               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    834 +                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
    835 +            END_2D 
    836           CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    837 -            DO jj = 2, jpj 
    838 -               DO ji = 2, jpi   ! vector opt. 
    839 -                  zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    840 -                     &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
    841 -               END DO 
    842 -            END DO 
    843 +            DO_2D( 2, jpi , 2, jpj ) 
    844 +               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    845 +                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
    846 +            END_2D 
    847           CASE ( np_CME )                           !* Coriolis + metric 
    848 -            DO jj = 2, jpj 
    849 -               DO ji = 2, jpi   ! vector opt. 
    850 -                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    851 -                       &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    852 -                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
    853 -               END DO 
    854 -            END DO 
    855 +            DO_2D( 2, jpi , 2, jpj ) 
    856 +               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    857 +                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    858 +                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
    859 +            END_2D 
    860           CASE DEFAULT                                             ! error 
    861              CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    862           END SELECT 
    863           ! 
    864           !                                   !==  compute and add the vorticity term trend  =! 
    865 -         DO jj = 2, jpjm1 
    866 -            DO ji = 2, jpim1   ! vector opt. 
    867 -               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
    868 -                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
    869 -                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
    870 -                  ! 
    871 -               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
    872 -                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
    873 -                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
    874 -            END DO 
    875 -         END DO 
    876 +         DO_2D( 2, jpim1 , 2, jpjm1 ) 
    877 +            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
    878 +               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
    879 +               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
    880 +               ! 
    881 +            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
    882 +               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
    883 +               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
    884 +         END_2D 
    885           !                                             ! =============== 
    886        END DO                                           !   End of slab 
    887        !                                                ! =============== 
    888  
    889 }}} 
    890  
    891 ''' Extending for the 3D loops ''' 
    892  
    893 Dealing 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  
    913 Only 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  
     1110and 
    9151111{{{ 
    9161112cat do3dfinder.pl 
     
    9291125      # 
    9301126      my $iline = <F> || die "DO jk line at end of file?"; 
    931       if ( $iline =~ /^\s*DO_2D\s*\(/i) { 
     1127      if ( $iline =~ /^\s*DO_2D_.*/i) { 
    9321128         my $isinavlid = 0; 
    9331129         # 
     
    9361132         my $docount = 0; 
    9371133         # 
    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); 
     1134         # 4. Store the loop limits from the k-loop line and remove spaces and new-lines 
     1135         # 
     1136         ($jargs = $jline) =~ s/(^.*?)=([^\!\n]*)(\!*.*)/\2/; 
    9431137         chomp($jargs); 
    944          $iargs =~ s/^\s+//; $iargs =~ s/\s+$//; 
    9451138         $jargs =~ s/^\s+//; $jargs =~ s/\s+$//; 
    9461139         # 
     
    9601153         print $orig $iline; 
    9611154         my $altr = new IO::Scalar \$astr; 
    962          print $altr $jspac,"DO_3D( ",$iargs," , ",$jargs," )\n"; 
     1155         $iline =~ s/DO_2D/DO_3D/; 
     1156         $iline =~ s/^\s+//; $iline =~ s/\s+$//; 
     1157         $iline =~ s/DO_2D/DO_3D/; chomp($iline); 
     1158         print $altr $jspac,$iline,"( ",$jargs," )\n"; 
    9631159         # 
    9641160         # 7. Now process the loop contents until the matching  END_2D statement 
     
    10221218} 
    10231219}}} 
    1024 These 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 }}} 
    1030 And 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) ) 
    1044 -            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) 
    1052 -            END DO 
    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 
    1084 -         DO jj = 2, jpjm1 
    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) 
    1088 -            END DO 
    1089 -         END DO 
    1090 +         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
    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) 
    1093 +         END_2D 
    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 }}} 
    11311220 
    11321221