Version 24 (modified by acc, 8 months ago) (diff)

KERNEL-02_Storkey_Coward_IMMERSE_first_steps

The PI is responsible to closely follow the progress of the action, and especially to contact NEMO project manager if the delay on preview (or review) are longer than the 2 weeks expected.

Help

The action has to be detailed briefly in the 'Summary' for later inclusion in other pages. To do so, edit 'Summary' section as a common wiki page and set links for the development ticket and branch.

Out of this, the rest of the page ('Abstract'|'Preview'|'Tests'|'Review') can be edited on-line inside the form fields considering the following color code given hereafter: PI(S), Previewer(s) and Reviewer(s).
Record your modifications for the section you have edited by clicking on the corresponding button at the end of the section with 'Save …' button. Just above, the log record will be updated.

The informations inside the form fields and this wiki page itself are stored in 2 separate databases. As a consequence, there is absolutely no risk to make any modification in the page itself as long as you don't rename the page or modify the source code of {{{#!TracForm ... }}} processors.

Summary

Action KERNEL-02_Storkey_Coward_IMMERSE_first_steps
PI(S) Dave Storkey, Andrew Coward

Digest

Reorganisation of code to prepare for implementation of RK3 timestepping and tiling (IMMERSE WP3 and WP4)

  1. State variables to be renamed with an extra time dimension and time level swapping achieved by swapping time level indices.
  2. Prognostic fields to be passed to tendency routines from stp as arguments.
  3. Introduction of DO loop macros in preparation for the implementation of tiling.

Report on 2019 progress (items 1. and 2.) now available: https://forge.ipsl.jussieu.fr/nemo/wiki/2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps/Report

Dependencies
Target
Trac Ticket #2258
SVN branch branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps
Previewer(s) Gurvan Madec
Reviewer(s) Gurvan Madec
Link

'.' => '/nemo/wiki/2019WP/KERNEL-02_Storkey_Coward_IMMERSE_first_steps'

Abstract

This section should be completed before starting to develop the code, in order to find agreement on the method beforehand.

Description

Implementation

Reference manual and web pages updates

Updated on 12/03/2019 16:49:00 by davestorkey

Once the PI has completed this section, he should send a mail to the previewer(s) asking them to preview the work within two weeks.

Preview

1. Variable names and extra dimensions
The refactoring script
The refactoring script explained
Notes on testing regular expressions
Some contrived tests
Results on real files

2. Do loop macro substitution
Automating the tiling changes
Results on traldf_iso.F90
do2dfinder.pl
do3dfinder.pl
Sanity checks and domain_substitute.h90

Part of the reorganisation for RK3 requires the refactoring of arrays such as un, ub into a single, 4 dimensional array with a time-level dimension. It is expected that much of the work required here can be automated to the extent that it is feasible to re-apply these changes after the annual merge. Below is a working example of how this might be achieved. Perl is used to carry out the pattern matching and substitution because of its ability to match patterns extending over several lines. A random subset of source files are used in this example and serve to illustrate the successes and caveats for the method.

Version 2 -Tweaked refactoring script to handle indirect addressing (i.e. brackets within array indices)

Step 1 Perl is used in a 'edit in place' mode so the original files will be overwritten. Step 1 is therefore to create copies of the test files:

#!/bin/bash
mkdir TEST_FILES
cp FLO/flo_oce.F90 FLO/floats.F90 SBC/sbcfwb.F90 DYN/dynadv_ubs.F90 DYN/dynkeg.F90 DYN/dynvor.F90 DYN/dynadv_cen2.F90 TRA/trabbl.F90 TEST_FILES
cp FLO/flo_oce.F90 FLO/floats.F90 SBC/sbcfwb.F90 DYN/dynadv_ubs.F90 DYN/dynkeg.F90 DYN/dynvor.F90 DYN/dynadv_cen2.F90 TRA/trabbl.F90 TEST_FILES_ORG

The refactoring script

#!/bin/bash
#
 INVARS=( ub  vb  wb  un  vn  wn )
OUTVARS=( uu  vv  ww  uu  vv  ww )
  TLEVS=( Nnn Nnn Nnn Nii Nii Nii )
#
rm patch.list
for f in TEST_FILES/*.F90
do
 echo "{{{#!diff" >> patch.list
 echo "Index: "$f >> patch.list
 echo "==============================" >> patch.list
 n=0
 for n in `seq 0 1 $(( ${#INVARS[*]} - 1 ))`
 do
  perl -0777 -pi -e 's@([+.(,\s\-\/\*\%])'${INVARS[$n]}'(\s*)(\(((?:[^()]++|(?3))*)\))@\1'${OUTVARS[$n]}'\2\(\4,'${TLEVS[$n]}'\)@g'  $f
 done
 diff -u TEST_FILES_ORG/`basename $f` $f >> patch.list
 echo "}}}" >> patch.list
done

The refactoring script explained

# Some bash arrays to list the old names, new names and associated time-level index. 
# The choice of names here is not meant to reflect a desired choice. Note all three arrays 
# must have the same number of entries

   INVARS=( ub  vb  wb  un  vn  wn )
  OUTVARS=( uu  vv  ww  uu  vv  ww )
    TLEVS=( Nnn Nnn Nnn Nii Nii Nii )

#
# Lines referring to patch.list are just generating some output for this Wiki page. 
# They can be ignored for the purposes of explaining the script

# Loop over all files in test directory

  for f in TEST_FILES/*.F90
  do

# Loop over each input variable name

   for n in `seq 0 1 $(( ${#INVARS[*]} - 1 ))`
   do

# The perl command:

# -0777 -pi  are the options that force replace in place operation. 
# Could elect to redirect to new files

# The substitute command matches the INVAR string preceded by any of: +.(,whitespace-/*%
# Any match here goes into the first pattern space (\1)

# The INVAR string can have any amount of whitespace (including none) between it and an 
# opening bracket. This whitespace is preserved in pattern space 2

# Everything following the opening bracket to matching closing bracket 
# is stored in pattern space 3. This requires a recursion and pattern space 4 ends up with
# the interior of the outer brackets

# The RHS of the substitute command rebuilds the line replacing INVARS with OUTVARS and
# adding a ,TLEVS before the closing bracket

   perl -0777 -pi -e 's@([+.(,\s\-\/\*\%])'${INVARS[$n]}'(\s*)(\(((?:[^()]++|(?3))*)\))@\1'${OUTVARS[$n]}'\2\(\4,'${TLEVS[$n]}'\)@g'  $f

# End of variable loop

   done
 
# End of file loop

  done

Notes on testing regular expressions

Testing and deciphering the regular expression used in the LHS of the perl substitute command is made easier by the availability of on-line testers. below is a screenshot from regex101.com which helps illustrate and explain the regular expression used here:

regex tester

Some contrived tests:

!cat TEST_FILES_ORG/contrived_tests.F90
  un(:,:,:)                    ! The simplest test. Should ==> uu(:,:,:,jtn)
   un ( ji   , jj,   : )        ! Check alternative simple case ==> uu ( ji   , jj,   :,jtn)
   a+ub(:,:,jk)                 ! Check preceeding operators are correctly recognised
   a-vb(:,:,jk)                 ! Check preceeding operators are correctly recognised
   a*vn(:,:,jk)                 ! Check preceeding operators are correctly recognised
   a/wb(:,:,jk)                 ! Check preceeding operators are correctly recognised
   a%wn(:,:,jk)                 ! Check preceeding operators are correctly recognised
 .OR.wn(:,:,jk)                 ! Check preceeding operators are correctly recognised
    (wn(:,:,jk) + wn(:,:,jk-1)) ! Check preceeding brackets are correctly recognised
   un ( ji+1 ,     &            ! Check that entries over
  &     jj, jk - 1 )            !                    multiple lines are handled
   wn( mi0(ii) , mj0(jj ))      ! Brackets within arguments may break [ does this occur?]
  pun(:,:,:)                    ! target preceded by non-whitespace or operator. Should be unchanged
  sbc_fwb(:,:,:)                ! target preceded by non-whitespace or operator. Should be unchanged

Result of the script on the contrived tests:

   uu(:,:,:,Nii)                    ! The simplest test. Should ==> uu(:,:,:,jtn)
   uu ( ji   , jj,   : ,Nii)        ! Check alternative simple case ==> uu ( ji   , jj,   :,jtn)
   a+uu(:,:,jk,Nnn)                 ! Check preceeding operators are correctly recognised
   a-vv(:,:,jk,Nnn)                 ! Check preceeding operators are correctly recognised
   a*vv(:,:,jk,Nii)                 ! Check preceeding operators are correctly recognised
   a/ww(:,:,jk,Nnn)                 ! Check preceeding operators are correctly recognised
   a%ww(:,:,jk,Nii)                 ! Check preceeding operators are correctly recognised
 .OR.ww(:,:,jk,Nii)                 ! Check preceeding operators are correctly recognised
    (ww(:,:,jk,Nii) + ww(:,:,jk-1,Nii)) ! Check preceeding brackets are correctly recognised
   uu ( ji+1 ,     &            ! Check that entries over
  &     jj, jk - 1 ,Nii)            !                    multiple lines are handled
   ww( mi0(ii) , mj0(jj ),Nii)      ! Brackets within arguments may break [ does this occur?]
  pun(:,:,:)                    ! target preceded by non-whitespace or operator. Should be unchanged
  sbc_fwb(:,:,:)                ! target preceded by non-whitespace or operator. Should be unchanged
  • contrived_tests.F90

    ==============================
    old new  
    1    un(:,:,:)                    ! The simplest test. Should ==> uu(:,:,:,jtn) 
    2    un ( ji   , jj,   : )        ! Check alternative simple case ==> uu ( ji   , jj,   :,jtn) 
    3    a+ub(:,:,jk)                 ! Check preceeding operators are correctly recognised 
    4    a-vb(:,:,jk)                 ! Check preceeding operators are correctly recognised 
    5    a*vn(:,:,jk)                 ! Check preceeding operators are correctly recognised 
    6    a/wb(:,:,jk)                 ! Check preceeding operators are correctly recognised 
    7    a%wn(:,:,jk)                 ! Check preceeding operators are correctly recognised 
    8  .OR.wn(:,:,jk)                 ! Check preceeding operators are correctly recognised 
    9     (wn(:,:,jk) + wn(:,:,jk-1)) ! Check preceeding brackets are correctly recognised 
    10    un ( ji+1 ,     &            ! Check that entries over 
    11   &     jj, jk - 1 )            !                    multiple lines are handled 
    12    wn( mi0(ii) , mj0(jj ))      ! Brackets within arguments may break [ does this occur?] 
     1   uu(:,:,:,Nii)                    ! The simplest test. Should ==> uu(:,:,:,jtn) 
     2   uu ( ji   , jj,   : ,Nii)        ! Check alternative simple case ==> uu ( ji   , jj,   :,jtn) 
     3   a+uu(:,:,jk,Nnn)                 ! Check preceeding operators are correctly recognised 
     4   a-vv(:,:,jk,Nnn)                 ! Check preceeding operators are correctly recognised 
     5   a*vv(:,:,jk,Nii)                 ! Check preceeding operators are correctly recognised 
     6   a/ww(:,:,jk,Nnn)                 ! Check preceeding operators are correctly recognised 
     7   a%ww(:,:,jk,Nii)                 ! Check preceeding operators are correctly recognised 
     8 .OR.ww(:,:,jk,Nii)                 ! Check preceeding operators are correctly recognised 
     9    (ww(:,:,jk,Nii) + ww(:,:,jk-1,Nii)) ! Check preceeding brackets are correctly recognised 
     10   uu ( ji+1 ,     &            ! Check that entries over 
     11  &     jj, jk - 1 ,Nii)            !                    multiple lines are handled 
     12   ww( mi0(ii) , mj0(jj ),Nii)      ! Brackets within arguments may break [ does this occur?] 
    1313  pun(:,:,:)                    ! target preceded by non-whitespace or operator. Should be unchanged 
    1414  sbc_fwb(:,:,:)                ! target preceded by non-whitespace or operator. Should be unchanged 

So all changes were made correctly and even those entries which were potential pitfalls (pun and sbc_fwb) were correctly ignored. Time to try a real set:

The results on the sample set of files (patch.list):

  • dynadv_cen2.F90

    ==============================
    old new  
    6666      !                             !==  Horizontal advection  ==! 
    6767      ! 
    6868      DO jk = 1, jpkm1                    ! horizontal transport 
    69          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    70          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     69         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) 
     70         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) 
    7171         DO jj = 1, jpjm1                 ! horizontal momentum fluxes (at T- and F-point) 
    7272            DO ji = 1, fs_jpim1   ! vector opt. 
    73                zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
    74                zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
    75                zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
    76                zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     73               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( uu(ji,jj,jk,Nii) + uu(ji+1,jj  ,jk,Nii) ) 
     74               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( uu(ji,jj,jk,Nii) + uu(ji  ,jj+1,jk,Nii) ) 
     75               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vv(ji,jj,jk,Nii) + vv(ji+1,jj  ,jk,Nii) ) 
     76               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vv(ji,jj,jk,Nii) + vv(ji  ,jj+1,jk,Nii) ) 
    7777            END DO 
    7878         END DO 
    7979         DO jj = 2, jpjm1                 ! divergence of horizontal momentum fluxes 
     
    105105      IF( ln_linssh ) THEN                ! linear free surface: advection through the surface 
    106106         DO jj = 2, jpjm1 
    107107            DO ji = fs_2, fs_jpim1 
    108                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) 
    109                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) 
     108               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1,Nii) + e1e2t(ji+1,jj) * ww(ji+1,jj,1,Nii) ) * uu(ji,jj,1,Nii) 
     109               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1,Nii) + e1e2t(ji,jj+1) * ww(ji,jj+1,1,Nii) ) * vv(ji,jj,1,Nii) 
    110110            END DO 
    111111         END DO 
    112112      ENDIF 
    113113      DO jk = 2, jpkm1                    ! interior advective fluxes 
    114114         DO jj = 2, jpj                       ! 1/4 * Vertical transport 
    115115            DO ji = 2, jpi 
    116                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     116               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk,Nii) 
    117117            END DO 
    118118         END DO 
    119119         DO jj = 2, jpjm1 
    120120            DO ji = fs_2, fs_jpim1   ! vector opt. 
    121                zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
    122                zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     121               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( uu(ji,jj,jk,Nii) + uu(ji,jj,jk-1,Nii) ) 
     122               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vv(ji,jj,jk,Nii) + vv(ji,jj,jk-1,Nii) ) 
    123123            END DO 
    124124         END DO 
    125125      END DO 
  • dynadv_ubs.F90

    ==============================
    old new  
    101101      DO jk = 1, jpkm1                       !  Laplacian of the velocity  ! 
    102102         !                                   ! =========================== ! 
    103103         !                                         ! horizontal volume fluxes 
    104          zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    105          zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     104         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) 
     105         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) 
    106106         ! 
    107107         DO jj = 2, jpjm1                          ! laplacian 
    108108            DO ji = fs_2, fs_jpim1   ! vector opt. 
    109                zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj  ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
    110                zlv_vv(ji,jj,jk,1) = ( vb (ji  ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji  ,jj-1,jk) ) * vmask(ji,jj,jk) 
    111                zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
    112                   &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk) 
    113                zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
    114                   &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk) 
     109               zlu_uu(ji,jj,jk,1) = ( uu (ji+1,jj  ,jk,Nnn) - 2.*uu (ji,jj,jk,Nnn) + uu (ji-1,jj  ,jk,Nnn) ) * umask(ji,jj,jk) 
     110               zlv_vv(ji,jj,jk,1) = ( vv (ji  ,jj+1,jk,Nnn) - 2.*vv (ji,jj,jk,Nnn) + vv (ji  ,jj-1,jk,Nnn) ) * vmask(ji,jj,jk) 
     111               zlu_uv(ji,jj,jk,1) = ( uu (ji  ,jj+1,jk,Nnn) - uu (ji  ,jj  ,jk,Nnn) ) * fmask(ji  ,jj  ,jk)   & 
     112                  &               - ( uu (ji  ,jj  ,jk,Nnn) - uu (ji  ,jj-1,jk,Nnn) ) * fmask(ji  ,jj-1,jk) 
     113               zlv_vu(ji,jj,jk,1) = ( vv (ji+1,jj  ,jk,Nnn) - vv (ji  ,jj  ,jk,Nnn) ) * fmask(ji  ,jj  ,jk)   & 
     114                  &               - ( vv (ji  ,jj  ,jk,Nnn) - vv (ji-1,jj  ,jk,Nnn) ) * fmask(ji-1,jj  ,jk) 
    115115               ! 
    116116               zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj  ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
    117117               zlv_vv(ji,jj,jk,2) = ( zfv(ji  ,jj+1,jk) - 2.*zfv(ji,jj,jk) + zfv(ji  ,jj-1,jk) ) * vmask(ji,jj,jk) 
     
    131131      !                                      !  Horizontal advection  ! 
    132132      DO jk = 1, jpkm1                       ! ====================== ! 
    133133         !                                         ! horizontal volume fluxes 
    134          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    135          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     134         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) 
     135         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) 
    136136         ! 
    137137         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point 
    138138            DO ji = 1, fs_jpim1   ! vector opt. 
    139                zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
    140                zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     139               zui = ( uu(ji,jj,jk,Nii) + uu(ji+1,jj  ,jk,Nii) ) 
     140               zvj = ( vv(ji,jj,jk,Nii) + vv(ji  ,jj+1,jk,Nii) ) 
    141141               ! 
    142142               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1) 
    143143               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1) 
     
    163163               ENDIF 
    164164               ! 
    165165               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   & 
    166                   &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u ) 
     166                  &                * ( uu(ji,jj,jk,Nii) + uu(ji  ,jj+1,jk,Nii) - gamma1 * zl_u ) 
    167167               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   & 
    168                   &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v ) 
     168                  &                * ( vv(ji,jj,jk,Nii) + vv(ji+1,jj  ,jk,Nii) - gamma1 * zl_v ) 
    169169            END DO 
    170170         END DO 
    171171         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes 
     
    198198      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface 
    199199         DO jj = 2, jpjm1 
    200200            DO ji = fs_2, fs_jpim1 
    201                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) 
    202                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) 
     201               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1,Nii) + e1e2t(ji+1,jj) * ww(ji+1,jj,1,Nii) ) * uu(ji,jj,1,Nii) 
     202               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1,Nii) + e1e2t(ji,jj+1) * ww(ji,jj+1,1,Nii) ) * vv(ji,jj,1,Nii) 
    203203            END DO 
    204204         END DO 
    205205      ENDIF 
    206206      DO jk = 2, jpkm1                          ! interior fluxes 
    207207         DO jj = 2, jpj 
    208208            DO ji = 2, jpi 
    209                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     209               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk,Nii) 
    210210            END DO 
    211211         END DO 
    212212         DO jj = 2, jpjm1 
    213213            DO ji = fs_2, fs_jpim1   ! vector opt. 
    214                zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
    215                zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     214               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( uu(ji,jj,jk,Nii) + uu(ji,jj,jk-1,Nii) ) 
     215               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vv(ji,jj,jk,Nii) + vv(ji,jj,jk-1,Nii) ) 
    216216            END DO 
    217217         END DO 
    218218      END DO 
  • dynkeg.F90

    ==============================
    old new  
    5656      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 
    5757      !!              * kscheme = nkeg_HW : Hollingsworth correction following 
    5858      !!      Arakawa (2001). The now horizontal kinetic energy is given by: 
    59       !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  ) 
    60       !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ] 
     59      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((uu(j+1,Nii)+uu(j-1,Nii))/2)^2  ) 
     60      !!                    + mj-1(  2 * vn^2 + ((vv(i+1,Nii)+vv(i-1,Nii))/2)^2  ) ] 
    6161      !! 
    6262      !!      Take its horizontal gradient and add it to the general momentum 
    6363      !!      trend (ua,va). 
     
    108108                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    109109                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    110110                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
    111                      un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
     111                     uu(ii-ifu,ij,jk,Nii) = uu(ii,ij,jk,Nii) * umask(ii,ij,jk) 
    112112                  END DO 
    113113               END DO 
    114114               ! 
     
    118118                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    119119                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    120120                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
    121                      vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
     121                     vv(ii,ij-ifv,jk,Nii) = vv(ii,ij,jk,Nii) * vmask(ii,ij,jk) 
    122122                  END DO 
    123123               END DO 
    124124            ENDIF 
     
    131131         DO jk = 1, jpkm1 
    132132            DO jj = 2, jpj 
    133133               DO ji = fs_2, jpi   ! vector opt. 
    134                   zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
    135                      &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) 
    136                   zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   & 
    137                      &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) 
     134                  zu =    uu(ji-1,jj  ,jk,Nii) * uu(ji-1,jj  ,jk,Nii)   & 
     135                     &  + uu(ji  ,jj  ,jk,Nii) * uu(ji  ,jj  ,jk,Nii) 
     136                  zv =    vv(ji  ,jj-1,jk,Nii) * vv(ji  ,jj-1,jk,Nii)   & 
     137                     &  + vv(ji  ,jj  ,jk,Nii) * vv(ji  ,jj  ,jk,Nii) 
    138138                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
    139139               END DO 
    140140            END DO 
     
    144144         DO jk = 1, jpkm1 
    145145            DO jj = 2, jpjm1 
    146146               DO ji = fs_2, jpim1   ! vector opt. 
    147                   zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    & 
    148                      &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  & 
    149                      &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   & 
    150                      &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) 
     147                  zu = 8._wp * ( uu(ji-1,jj  ,jk,Nii) * uu(ji-1,jj  ,jk,Nii)    & 
     148                     &         + uu(ji  ,jj  ,jk,Nii) * uu(ji  ,jj  ,jk,Nii) )  & 
     149                     &   +     ( uu(ji-1,jj-1,jk,Nii) + uu(ji-1,jj+1,jk,Nii) ) * ( uu(ji-1,jj-1,jk,Nii) + uu(ji-1,jj+1,jk,Nii) )   & 
     150                     &   +     ( uu(ji  ,jj-1,jk,Nii) + uu(ji  ,jj+1,jk,Nii) ) * ( uu(ji  ,jj-1,jk,Nii) + uu(ji  ,jj+1,jk,Nii) ) 
    151151                     ! 
    152                   zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    & 
    153                      &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  & 
    154                      &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   & 
    155                      &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) 
     152                  zv = 8._wp * ( vv(ji  ,jj-1,jk,Nii) * vv(ji  ,jj-1,jk,Nii)    & 
     153                     &         + vv(ji  ,jj  ,jk,Nii) * vv(ji  ,jj  ,jk,Nii) )  & 
     154                     &  +      ( vv(ji-1,jj-1,jk,Nii) + vv(ji+1,jj-1,jk,Nii) ) * ( vv(ji-1,jj-1,jk,Nii) + vv(ji+1,jj-1,jk,Nii) )   & 
     155                     &  +      ( vv(ji-1,jj  ,jk,Nii) + vv(ji+1,jj  ,jk,Nii) ) * ( vv(ji-1,jj  ,jk,Nii) + vv(ji+1,jj  ,jk,Nii) ) 
    156156                  zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 
    157157               END DO 
    158158            END DO 
     
    163163 
    164164      IF (ln_bdy) THEN 
    165165         ! restore velocity masks at points outside boundary 
    166          un(:,:,:) = un(:,:,:) * umask(:,:,:) 
    167          vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
     166         uu(:,:,:,Nii) = uu(:,:,:,Nii) * umask(:,:,:) 
     167         vv(:,:,:,Nii) = vv(:,:,:,Nii) * vmask(:,:,:) 
    168168      ENDIF 
    169169 
    170170      ! 
Index: TEST_FILES/dynvor.F90
==============================
  • flo_oce.F90

    ==============================
    old new  
    5959      !!---------------------------------------------------------------------- 
    6060      !!                 ***  FUNCTION flo_oce_alloc  *** 
    6161      !!---------------------------------------------------------------------- 
    62       ALLOCATE( wb(jpi,jpj,jpk) , nfloat(jpnfl) , nisobfl(jpnfl) , ngrpfl(jpnfl) , & 
     62      ALLOCATE( ww(jpi,jpj,jpk,Nnn) , nfloat(jpnfl) , nisobfl(jpnfl) , ngrpfl(jpnfl) , & 
    6363         &      flxx(jpnfl)     , flyy(jpnfl)   , flzz(jpnfl)    ,                 & 
    6464         &      tpifl(jpnfl)    , tpjfl(jpnfl)  , tpkfl(jpnfl)   , STAT=flo_oce_alloc ) 
    6565      ! 

Ok This was wrong.Nnn is not what should go into the ALLOCATE statement

  • floats.F90

    ==============================
    old new  
    6464      ! 
    6565      CALL flo_rst( kt )      ! trajectories restart 
    6666      ! 
    67       wb(:,:,:) = wn(:,:,:)         ! Save the old vertical velocity field 
     67      ww(:,:,:,Nnn) = ww(:,:,:,Nii)         ! Save the old vertical velocity field 
    6868      ! 
    6969      IF( ln_timing )   CALL timing_stop('flo_stp') 
    7070      ! 
     
    131131      ! 
    132132      CALL flo_dom                  ! compute/read initial position of floats 
    133133      ! 
    134       wb(:,:,:) = wn(:,:,:)         ! set wb for computation of floats trajectories at the first time step 
     134      ww(:,:,:,Nnn) = ww(:,:,:,Nii)         ! set wb for computation of floats trajectories at the first time step 
    135135      ! 
    136136   END SUBROUTINE flo_init 
    137137 
Index: TEST_FILES/sbcfwb.F90
==============================
  • trabbl.F90

    ==============================
    old new  
    347347            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 
    348348            ! 
    349349            zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth 
    350             zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity 
    351             zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
     350            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Nii)          ! bottom velocity 
     351            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Nii) 
    352352         END DO 
    353353      END DO 
    354354      ! 

So far so good….

Automating the tiling changes

Here 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:

!cat TESTDO_FILES/testdo.F90
!   some random text
!   followed by a valid loop pair
     DO jj  = 2,jpjm1     ! with comments
        DO ji = 2,jpim1
           some loop content
           more loop content
        END DO
     END DO
!   followed by an invalid loop pair
     DO jj  = 2,jpjm1
        j = jj-1

        DO ji = 2,jpim1
           some loop content
           more loop content
        END DO
     END DO
!   followed by a valid loop pair with a nested do
     DO jj  = 1,jpjm1
        DO ji = 1,jpi
           some loop content
           do jn = 1, jptrc
              yet more loop content
           end do
           more loop content
        END DO
     END DO
!   followed by a valid 3D loop
     DO jk  = 2,jpkm1
        DO jj  = 1,jpjm1
           DO ji = 1,jpi
              some loop content
              more loop content
           END DO
        END DO
     END DO

This file is processed as follows:

  perl do2dfinder.pl TESTDO_FILES/testdo.F90 > TESTDO_FILES_2D/testdo.F90
  perl do3dfinder.pl TESTDO_FILES_2D/testdo.F90 > TESTDO_FILES_3D/testdo.F90
  diff -u TESTDO_FILES/testdo.F90 TESTDO_FILES_3D/testdo.F90 > testdo.patch

and the resulting changes are:

  • testdo.F90

    old new  
    11!   some random text 
    22!   followed by a valid loop pair 
    3      DO jj  = 2,jpjm1     ! with comments 
    4         DO ji = 2,jpim1 
    5            some loop content 
    6            more loop content 
    7         END DO 
    8      END DO 
     3     DO_2D_00_00 
     4        some loop content 
     5        more loop content 
     6     END_2D 
    97!   followed by an invalid loop pair 
    108     DO jj  = 2,jpjm1 
    119        j = jj-1 
     
    1614        END DO 
    1715     END DO 
    1816!   followed by a valid loop pair with a nested do 
    19      DO jj  = 1,jpjm1 
    20         DO ji = 1,jpi 
    21            some loop content 
    22            do jn = 1, jptrc 
    23               yet more loop content 
    24            end do 
    25            more loop content 
    26         END DO 
    27      END DO 
     17     DO_2D_10_11 
     18        some loop content 
     19        do jn = 1, jptrc 
     20           yet more loop content 
     21        end do 
     22        more loop content 
     23     END_2D 
    2824!   followed by a valid 3D loop 
    29      DO jk  = 2,jpkm1 
    30         DO jj  = 1,jpjm1 
    31            DO ji = 1,jpi 
    32               some loop content 
    33               more loop content 
    34            END DO 
    35         END DO 
    36      END DO 
     25     DO_3D_10_11( 2,jpkm1 ) 
     26        some loop content 
     27        more loop content 
     28     END_3D 

traldf_iso.F90 provides a more stringent test:

   perl do2dfinder.pl TESTDO_FILES/traldf_iso.F90 > TESTDO_FILES_2D/traldf_iso.F90
   perl do3dfinder.pl TESTDO_FILES_2D/traldf_iso.F90 > TESTDO_FILES_3D/traldf_iso.F90
   diff -u TESTDO_FILES/traldf_iso.F90 TESTDO_FILES_3D/traldf_iso.F90 > traldf_iso.patch
  • traldf_iso.F90

    old new  
    143143      ! 
    144144      IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
    145145         ! 
    146          DO jk = 2, jpkm1 
    147             DO jj = 2, jpjm1 
    148                DO ji = fs_2, fs_jpim1   ! vector opt. 
    149                   ! 
    150                   zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    151                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    152                   zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    153                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    154                      ! 
    155                   zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    156                      &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    157                   zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    158                      &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
    159                      ! 
    160                   ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    161                      &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
    162                END DO 
    163             END DO 
    164          END DO 
     146         DO_3D_00_00( 2, jpkm1 ) 
     147            ! 
     148            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     149               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     150            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     151               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     152               ! 
     153            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     154               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     155            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     156               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     157               ! 
     158            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     159               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     160         END_3D 
    165161         ! 
    166162         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
    167             DO jk = 2, jpkm1 
    168                DO jj = 2, jpjm1 
    169                   DO ji = fs_2, fs_jpim1 
    170                      akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
    171                         &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
    172                         &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
    173                         &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
    174                         &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
    175                   END DO 
    176                END DO 
    177             END DO 
     163            DO_3D_00_00( 2, jpkm1 ) 
     164               akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     165                  &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     166                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     167                  &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     168                  &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     169            END_3D 
    178170            ! 
    179171            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
    180                DO jk = 2, jpkm1 
    181                   DO jj = 1, jpjm1 
    182                      DO ji = 1, fs_jpim1 
    183                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    184                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
    185                      END DO 
    186                   END DO 
    187                END DO 
     172               DO_3D_10_10( 2, jpkm1 ) 
     173                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     174                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     175               END_3D 
    188176            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    189                DO jk = 2, jpkm1 
    190                   DO jj = 1, jpjm1 
    191                      DO ji = 1, fs_jpim1 
    192                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
    193                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    194                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
    195                      END DO 
    196                   END DO 
    197                END DO 
     177               DO_3D_10_10( 2, jpkm1 ) 
     178                  ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     179                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     180                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     181               END_3D 
    198182           ENDIF 
    199183           ! 
    200184         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
     
    215199         !!end 
    216200 
    217201         ! Horizontal tracer gradient 
    218          DO jk = 1, jpkm1 
    219             DO jj = 1, jpjm1 
    220                DO ji = 1, fs_jpim1   ! vector opt. 
    221                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    222                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    223                END DO 
    224             END DO 
    225          END DO 
     202         DO_3D_10_10( 1, jpkm1 ) 
     203            zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     204            zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     205         END_3D 
    226206         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
    227             DO jj = 1, jpjm1              ! bottom correction (partial bottom cell) 
    228                DO ji = 1, fs_jpim1   ! vector opt. 
    229                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    230                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    231                END DO 
    232             END DO 
     207            DO_2D_10_10 
     208               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
     209               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     210            END_2D 
    233211            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
    234                DO jj = 1, jpjm1 
    235                   DO ji = 1, fs_jpim1   ! vector opt. 
    236                      IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
    237                      IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
    238                   END DO 
    239                END DO 
     212               DO_2D_10_10 
     213                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
     214                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
     215               END_2D 
    240216            ENDIF 
    241217         ENDIF 
    242218         ! 
     
    252228            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
    253229            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
    254230            ENDIF 
    255             DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    256                DO ji = 1, fs_jpim1   ! vector opt. 
    257                   zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
    258                   zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
    259                   ! 
    260                   zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
    261                      &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    262                   ! 
    263                   zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
    264                      &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    265                   ! 
    266                   zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    267                   zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
    268                   ! 
    269                   zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    270                      &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
    271                      &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    272                   zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    273                      &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    274                      &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
    275                END DO 
    276             END DO 
    277             ! 
    278             DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
    279                DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
    281                      &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    282                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    283                END DO 
    284             END DO 
     231            DO_2D_10_10 
     232               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
     233               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     234               ! 
     235               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     236                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     237               ! 
     238               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     239                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     240               ! 
     241               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     242               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     243               ! 
     244               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     245                  &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     246                  &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     247               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     248                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     249                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
     250            END_2D 
     251            ! 
     252            DO_2D_00_00 
     253               pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     254                  &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     255                  &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     256            END_2D 
    285257         END DO                                        !   End of slab 
    286258 
    287259         !!---------------------------------------------------------------------- 
     
    295267         !                          ! Surface and bottom vertical fluxes set to zero 
    296268         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    297269 
    298          DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1) 
    299             DO jj = 2, jpjm1 
    300                DO ji = fs_2, fs_jpim1   ! vector opt. 
    301                   ! 
    302                   zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    303                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    304                   zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    305                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    306                      ! 
    307                   zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    308                      &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    309                   zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    310                      &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
    311                      ! 
    312                   zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
    313                   zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
    314                   ! 
    315                   ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
    316                      &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
    317                      &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
    318                      &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
    319                END DO 
    320             END DO 
    321          END DO 
     270         DO_3D_00_00( 2, jpkm1 ) 
     271            ! 
     272            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     273               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     274            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     275               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     276               ! 
     277            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     278               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     279            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     280               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     281               ! 
     282            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     283            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     284            ! 
     285            ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     286               &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     287               &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     288               &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     289         END_3D 
    322290         !                                !==  add the vertical 33 flux  ==! 
    323291         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    324             DO jk = 2, jpkm1 
    325                DO jj = 1, jpjm1 
    326                   DO ji = fs_2, fs_jpim1 
    327                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
    328                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    329                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
    330                   END DO 
    331                END DO 
    332             END DO 
     292            DO_3D_10_00( 2, jpkm1 ) 
     293               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     294                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
     295                  &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     296            END_3D 
    333297            ! 
    334298         ELSE                                   ! bilaplacian 
    335299            SELECT CASE( kpass ) 
    336300            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    337                DO jk = 2, jpkm1 
    338                   DO jj = 1, jpjm1 
    339                      DO ji = fs_2, fs_jpim1 
    340                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
    341                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    342                            &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    343                      END DO 
    344                   END DO 
    345                END DO 
     301               DO_3D_10_00( 2, jpkm1 ) 
     302                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
     303                     &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
     304                     &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     305               END_3D 
    346306            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
    347                DO jk = 2, jpkm1 
    348                   DO jj = 1, jpjm1 
    349                      DO ji = fs_2, fs_jpim1 
    350                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
    351                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    352                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
    353                      END DO 
    354                   END DO 
    355                END DO 
     307               DO_3D_10_00( 2, jpkm1 ) 
     308                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
     309                     &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
     310                     &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     311               END_3D 
    356312            END SELECT 
    357313         ENDIF 
    358314         ! 
    359          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
    360             DO jj = 2, jpjm1 
    361                DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
    363                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    364                END DO 
    365             END DO 
    366          END DO 
     315         DO_3D_00_00( 1, jpkm1 ) 
     316            pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     317               &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     318         END_3D 
    367319         ! 
    368320         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
    369321             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==! 

And finally the two scripts that achieve this. Note the logic that may need tightening at sections 5 and 6 in the first script:

#cat do2dfinder.pl
#
open(F,$ARGV[0]) || die "Cannot open $ARGV[0]: $!";
while(<F>) {
   if ( $_ =~ /^\s*DO\s* jj/i) {
      # Start processing loop if line contains DO jj (case and whitespace independent)
      #
      # 1. Store the current line
      #
      $jline = $_;
      #
      # 2. Read the next line and check if it contains DO ji
      #
      my $iline = <F> || die "DO jj line at end of file?";
      if ( $iline =~ /^\s*do\s*ji/i) {
         #
         # 3. Initialise a count to track any nested do loops
         #
         my $docount = 0;
         #
         # 4. Store the loop limits from the two lines stored and remove spaces and new-lines
         #
         ($jargs = $jline) =~ s/(^.*?)=([^\!\n]*)(\!*.*)/\2/;
         ($iargs = $iline) =~ s/(^.*?)=([^\!\n]*)(\!*.*)/\2/;
         chomp($iargs);
         chomp($jargs);
         $iargs =~ s/^\s+//; $iargs =~ s/\s+$//;
         $jargs =~ s/^\s+//; $jargs =~ s/\s+$//;
         @ilims = split(/,/,$iargs);
         @jlims = split(/,/,$jargs);
         #
         # 5. Convert i-limits into Gurvan's indexing
         # i.e. lower limit: [2,fs_2] --> 0 ; 1--> 1
         # i.e. upper limit: *m* --> 0 else --> 1
         #
         if( @ilims[0] =~ /2/ ) {
           @ilims[0] = 0;
         } elsif ( @ilims[0] == 1 ) {
           @ilims[0] = 1;
         }
         if( @ilims[1] =~ /m/ ) {
           @ilims[1] = 0;
         } else {
           @ilims[1] = 1;
         }
         $iargs = join('',@ilims[0], @ilims[1]);
         #
         # 6. Convert j-limits into Gurvan's indexing
         # i.e. lower limit: [2,fs_2] --> 0 ; 1--> 1
         # i.e. upper limit: *m* --> 0 else --> 1
         #
         if( @jlims[0] =~ /2/ ) {
           @jlims[0] = 0;
         } elsif ( @jlims[0] == 1 ) {
           @jlims[0] = 1;
         }
         if( @jlims[1] =~ /m/ ) {
           @jlims[1] = 0;
         } else {
           @jlims[1] = 1;
         }
         $jargs = join('',@jlims[0], @jlims[1]);
         #
         # 7. Store the leading indentation for the outer loop
         #
         ($jspac = $jline) =~ s/(^[\s]*)([^\s]*).*/\1/;
         chomp($jspac);
         #
         # 8. Construct a DO_2D line to replace the original statements
         #
         print $jspac,join('_',"DO_2D",$jargs,$iargs),"\n";
         #
         # 9. Now process the loop contents until the matching pair of END DO statements
         #
         while ( $docount >= 0 || ! ( $iline =~ /^\s*end\s*do/i ) ) {
            $iline = <F> || die "reached end of file within do loop?" ;
            #
            # 10. Increment a counter if another DO statement is found
            #
            if ( $iline =~ /^\s*do\s*/i )  { $docount++ };
            #
            # 11. Decrement a counter if a END DO statement is found
            #
            if ( $iline =~ /^\s*end\s*do/i )  { $docount-- };
            #
            # 12. A negative counter means the matching END DO for the ji loop has been reached
            #
            if ( $docount < 0 ) {
               #
               # 13. Check the next line is the expected 2nd END DO.
               #     Output END_2D statement if it is
               #
               $jline = <F> || die "reached end of file looking for end do?" ;
               if ( ! ($jline =~ /^\s*end\s*do/i) )  {
                  die "END DOs are not consecutive";
               } else {
                  print $jspac,"END_2D\n";
               }

            } else {
               #
               # 14. This is a line inside the loop. Remove three leading spaces (if any) and output.
               #
               $iline =~ s/^\s\s\s//;
               print $iline;
            }
         }
      } else {
         #
         # 15. Consecutive DO statements were not found. Do not process these loops.
         #
         print $jline;
         print $iline;
      }
   } else {
      #
      # 16. Code outside of a DO construct. Leave unchanged.
      #
      print $_;
   }
}

and

#cat do3dfinder.pl
#
use IO::Scalar;
open(F,$ARGV[0]) || die "Cannot open $ARGV[0]: $!";
while(<F>) {
   if ( $_ =~ /^\s*DO\s* jk/i) {
      # Start processing loop if line contains DO jk (case and whitespace independent)
      #
      # 1. Store the current line
      #
      $jline = $_;
      #
      # 2. Read the next line and check if it contains DO_2D
      #
      my $iline = <F> || die "DO jk line at end of file?";
      if ( $iline =~ /^\s*DO_2D_.*/i) {
         my $isinavlid = 0;
         #
         # 3. Initialise a count to track any nested do loops
         #
         my $docount = 0;
         #
         # 4. Store the loop limits from the k-loop line and remove spaces and new-lines
         #
         ($jargs = $jline) =~ s/(^.*?)=([^\!\n]*)(\!*.*)/\2/;
         chomp($jargs);
         $jargs =~ s/^\s+//; $jargs =~ s/\s+$//;
         #
         # 5. Store the leading indentation for the outer loop
         #
         ($jspac = $jline) =~ s/(^[\s]*)([^\s]*).*/\1/;
         chomp($jspac);
         #
         # 6. Construct a DO_3D line to replace the original statements
         #
         # Keep two versions of output until we know if it is transformable
         #
         my $ostr = "";
         my $astr = "";
         my $orig = new IO::Scalar \$ostr;
         print $orig $jline;
         print $orig $iline;
         my $altr = new IO::Scalar \$astr;
         $iline =~ s/DO_2D/DO_3D/;
         $iline =~ s/^\s+//; $iline =~ s/\s+$//;
         $iline =~ s/DO_2D/DO_3D/; chomp($iline);
         print $altr $jspac,$iline,"( ",$jargs," )\n";
         #
         # 7. Now process the loop contents until the matching  END_2D statement
         #
         while ( $docount >= 0 || ! ( $iline =~ /^\s*END_2D/i ) ) {
            $iline = <F> || eval{ $isinvalid = 1 };
            #print $orig $iline;
            #
            # 8. Increment a counter if another DO_2D statement is found
            #
            if ( $iline =~ /^\s*do_2d/i )  { $docount++ };
            #
            # 9. Decrement a counter if a END DO statement is found
            #
            if ( $iline =~ /^\s*end_2d/i )  { $docount-- };
            #
            # 10. A negative counter means the matching END_2D for the ji loop has been reached
            #
            if ( $docount < 0 ) {
               #
               # 11. Check the next line is the expected END DO for the jk loop.
               #     Output END_3D statement if it is
               #
               $jline = <F> || eval {$isinvalid = 1} ;
               if ( ! ($jline =~ /^\s*end\s*do/i) )  {
                  $isinvalid = 1 ;
                  print $orig $iline;
                  print $orig $jline;
               } else {
                  print $altr $jspac,"END_3D\n";
               }
               if ( $isinvalid == 0 ) {
                  print $altr;
               } else {
                  print $orig;
                  $isinvalid = 0;
               }

            } else {
               #
               # 12. This is a line inside the loop. Remove three leading spaces (if any) and output.
               #
               print $orig $iline;
               $iline =~ s/^\s\s\s//;
               print $altr $iline;
            }
         }
      } else {
         #
         # 13. Consecutive DO statements were not found. Do not process these loops.
         #
         print $jline;
         print $iline;
      }
   } else {
      #
      # 14. Code outside of a DO construct. Leave unchanged.
      #
      print $_;
   }
}

Sanity Check

Introducing a form of the proposed domain_substitute.h90 file to the final version and running through a preprocssor should recover the equivalent of the original file (barring white space changes and line concatenations). For example:

cat domain_substitute.h90
#define kJs 2
#define kIs 2
#define kJe jpjm1
#define kIe jpim1
#define DO_2D_00_00 DO jj = kJs ,kJe     ; DO ji = kIs ,kIe
#define DO_2D_10_10 DO jj = kJs-1, kJe   ; DO ji = kIs-1, kIe
#define DO_2D_01_01 DO jj = kJs , kJe+1  ; DO ji = kIs , kIe+1
#define DO_2D_11_11 DO jj = kJs-1, kJe+1 ; DO ji = kIs-1, kIe+1

#define DO_3D_00_00(ks,ke) DO jk = ks, ke ; DO_2D_00_00
#define DO_3D_10_10(ks,ke) DO jk = ks, ke ; DO_2D_10_10
#define DO_3D_01_01(ks,ke) DO jk = ks, ke ; DO_2D_01_01
#define DO_3D_11_11(ks,ke) DO jk = ks, ke ; DO_2D_11_11

#define END_2D END DO ; END DO
#define END_3D END DO ; END DO ; END DO

ed - TESTDO_FILES_3D/traldf_iso.F90 << EOF
> 0a
> #include "domain_substute.h90"
> .
> w
> q
> EOF

gfortran -E -P TESTDO_FILES_3D/traldf_iso.F90 > SANITY/traldf_iso.f90
diff -u TESTDO_FILES/traldf_iso.F90 SANITY/traldf_iso.f90 > sanity.patch

And this does appear to be the case:

  • (a) TESTDO_FILES/traldf_iso.F90 vs. (b) SANITY/traldf_iso.f90

    a b  
     1 
     2 
     3 
    14MODULE traldf_iso 
    25   !!====================================================================== 
    36   !!                   ***  MODULE  traldf_iso  *** 
     
    3942   LOGICAL  ::   l_hst   ! flag to compute heat transport 
    4043 
    4144   !! * Substitutions 
    42 #  include "vectopt_loop_substitute.h90" 
     45   !!---------------------------------------------------------------------- 
     46   !!                   ***  vectopt_loop_substitute  *** 
     47   !!---------------------------------------------------------------------- 
     48   !! ** purpose :   substitute the inner loop start/end indices with CPP macro 
     49   !!                allow unrolling of do-loop (useful with vector processors) 
     50   !!---------------------------------------------------------------------- 
     51   !!---------------------------------------------------------------------- 
     52   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     53   !! $Id: vectopt_loop_substitute.h90 10068 2018-08-28 14:09:04Z nicolasmartin $ 
     54   !! Software governed by the CeCILL license (see ./LICENSE) 
     55   !!---------------------------------------------------------------------- 
    4356   !!---------------------------------------------------------------------- 
    4457   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    4558   !! $Id: traldf_iso.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ 
     
    143156      ! 
    144157      IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
    145158         ! 
    146          DO jk = 2, jpkm1 
    147             DO jj = 2, jpjm1 
    148                DO ji = fs_2, fs_jpim1   ! vector opt. 
    149                   ! 
    150                   zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    151                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    152                   zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    153                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    154                      ! 
    155                   zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    156                      &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    157                   zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    158                      &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
    159                      ! 
    160                   ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    161                      &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
    162                END DO 
    163             END DO 
    164          END DO 
     159         DO jk =  2,  jpkm1  ; DO jj = 2 ,jpjm1     ; DO ji = 2 ,jpim1 
     160            ! 
     161            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     162               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     163            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     164               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     165               ! 
     166            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     167               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     168            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     169               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     170               ! 
     171            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     172               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     173         END DO ; END DO ; END DO 
    165174         ! 
    166175         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
    167             DO jk = 2, jpkm1 
    168                DO jj = 2, jpjm1 
    169                   DO ji = fs_2, fs_jpim1 
    170                      akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
    171                         &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
    172                         &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
    173                         &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
    174                         &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
    175                   END DO 
    176                END DO 
    177             END DO 
     176            DO jk =  2,  jpkm1  ; DO jj = 2 ,jpjm1     ; DO ji = 2 ,jpim1 
     177               akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
     178                  &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     179                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
     180                  &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
     181                  &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     182            END DO ; END DO ; END DO 
    178183            ! 
    179184            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
    180                DO jk = 2, jpkm1 
    181                   DO jj = 1, jpjm1 
    182                      DO ji = 1, fs_jpim1 
    183                         akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    184                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
    185                      END DO 
    186                   END DO 
    187                END DO 
     185               DO jk =  2,  jpkm1  ; DO jj = 2-1, jpjm1   ; DO ji = 2-1, jpim1 
     186                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
     187                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     188               END DO ; END DO ; END DO 
    188189            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    189                DO jk = 2, jpkm1 
    190                   DO jj = 1, jpjm1 
    191                      DO ji = 1, fs_jpim1 
    192                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
    193                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    194                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
    195                      END DO 
    196                   END DO 
    197                END DO 
     190               DO jk =  2,  jpkm1  ; DO jj = 2-1, jpjm1   ; DO ji = 2-1, jpim1 
     191                  ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     192                  zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     193                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     194               END DO ; END DO ; END DO 
    198195           ENDIF 
    199196           ! 
    200197         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
     
    215212         !!end 
    216213 
    217214         ! Horizontal tracer gradient 
    218          DO jk = 1, jpkm1 
    219             DO jj = 1, jpjm1 
    220                DO ji = 1, fs_jpim1   ! vector opt. 
    221                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    222                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    223                END DO 
    224             END DO 
    225          END DO 
     215         DO jk =  1,  jpkm1  ; DO jj = 2-1, jpjm1   ; DO ji = 2-1, jpim1 
     216            zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     217            zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     218         END DO ; END DO ; END DO 
    226219         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
    227             DO jj = 1, jpjm1              ! bottom correction (partial bottom cell) 
    228                DO ji = 1, fs_jpim1   ! vector opt. 
    229                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    230                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    231                END DO 
    232             END DO 
     220            DO jj = 2-1, jpjm1   ; DO ji = 2-1, jpim1 
     221               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
     222               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     223            END DO ; END DO 
    233224            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
    234                DO jj = 1, jpjm1 
    235                   DO ji = 1, fs_jpim1   ! vector opt. 
    236                      IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
    237                      IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
    238                   END DO 
    239                END DO 
     225               DO jj = 2-1, jpjm1   ; DO ji = 2-1, jpim1 
     226                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
     227                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
     228               END DO ; END DO 
    240229            ENDIF 
    241230         ENDIF 
    242231         ! 
     
    252241            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
    253242            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
    254243            ENDIF 
    255             DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    256                DO ji = 1, fs_jpim1   ! vector opt. 
    257                   zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
    258                   zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
    259                   ! 
    260                   zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
    261                      &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    262                   ! 
    263                   zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
    264                      &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    265                   ! 
    266                   zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    267                   zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
    268                   ! 
    269                   zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    270                      &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
    271                      &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    272                   zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    273                      &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    274                      &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
    275                END DO 
    276             END DO 
     244            DO jj = 2-1, jpjm1   ; DO ji = 2-1, jpim1 
     245               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
     246               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     247               ! 
     248               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     249                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     250               ! 
     251               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     252                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     253               ! 
     254               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
     255               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
     256               ! 
     257               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
     258                  &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
     259                  &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     260               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
     261                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
     262                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
     263            END DO ; END DO 
    277264            ! 
    278             DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
    279                DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
    281                      &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    282                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    283                END DO 
    284             END DO 
     265            DO jj = 2 ,jpjm1     ; DO ji = 2 ,jpim1 
     266               pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     267                  &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
     268                  &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     269            END DO ; END DO 
    285270         END DO                                        !   End of slab 
    286271 
    287272         !!---------------------------------------------------------------------- 
     
    295280         !                          ! Surface and bottom vertical fluxes set to zero 
    296281         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    297282 
    298          DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1) 
    299             DO jj = 2, jpjm1 
    300                DO ji = fs_2, fs_jpim1   ! vector opt. 
    301                   ! 
    302                   zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    303                      &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    304                   zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    305                      &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    306                      ! 
    307                   zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    308                      &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    309                   zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    310                      &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
    311                      ! 
    312                   zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
    313                   zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
    314                   ! 
    315                   ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
    316                      &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
    317                      &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
    318                      &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
    319                END DO 
    320             END DO 
    321          END DO 
     283         DO jk =  2,  jpkm1  ; DO jj = 2 ,jpjm1     ; DO ji = 2 ,jpim1 
     284            ! 
     285            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     286               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
     287            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     288               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
     289               ! 
     290            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
     291               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
     292            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
     293               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     294               ! 
     295            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked 
     296            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
     297            ! 
     298            ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
     299               &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
     300               &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
     301               &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     302         END DO ; END DO ; END DO 
    322303         !                                !==  add the vertical 33 flux  ==! 
    323304         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    324             DO jk = 2, jpkm1 
    325                DO jj = 1, jpjm1 
    326                   DO ji = fs_2, fs_jpim1 
    327                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
    328                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    329                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
    330                   END DO 
    331                END DO 
    332             END DO 
     305            DO_3D_10_00( 2, jpkm1 ) 
     306               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     307                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
     308                  &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     309            END DO ; END DO ; END DO 
    333310            ! 
    334311         ELSE                                   ! bilaplacian 
    335312            SELECT CASE( kpass ) 
    336313            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    337                DO jk = 2, jpkm1 
    338                   DO jj = 1, jpjm1 
    339                      DO ji = fs_2, fs_jpim1 
    340                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
    341                            &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    342                            &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    343                      END DO 
    344                   END DO 
    345                END DO 
     314               DO_3D_10_00( 2, jpkm1 ) 
     315                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
     316                     &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
     317                     &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     318               END DO ; END DO ; END DO 
    346319            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
    347                DO jk = 2, jpkm1 
    348                   DO jj = 1, jpjm1 
    349                      DO ji = fs_2, fs_jpim1 
    350                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
    351                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    352                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
    353                      END DO 
    354                   END DO 
    355                END DO 
     320               DO_3D_10_00( 2, jpkm1 ) 
     321                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
     322                     &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
     323                     &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     324               END DO ; END DO ; END DO 
    356325            END SELECT 
    357326         ENDIF 
    358327         ! 
    359          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
    360             DO jj = 2, jpjm1 
    361                DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
    363                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    364                END DO 
    365             END DO 
    366          END DO 
     328         DO jk =  1,  jpkm1  ; DO jj = 2 ,jpjm1     ; DO ji = 2 ,jpim1 
     329            pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     330               &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     331         END DO ; END DO ; END DO 
    367332         ! 
    368333         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==! 
    369334             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==! 

Since the preview step must be completed before the PI starts the coding, the previewer(s) answers are expected to be completed within the two weeks after the PI has sent his request.
For each question, an iterative process should take place between PI and previewer(s) in order to reach a "YES" answer for each of the following questions.

Questions Answer Comment
Does the previewer agree with the proposed methodology?
Does the previewer agree with the proposed flowchart and list of routines to be changed?
Does the previewer agree with the proposed new list of variables, including agreement with coding rules?
Does the previewer agree with the proposed summary of updates in reference manual?
… … …

Updated on 08/04/2020 02:31:12 by anonymous

Once all "YES" have been reached, the PI can start the development into his development branch.

Tests

Once the development is done, the PI should complete this section below and ask the reviewers to start their review in the lower section.

Questions Answer Comment
Can this change be shown to produce expected impact? (if option activated)?
Can this change be shown to have a null impact? (if option not activated)
Detailed results of restartability and reproducibility when the option is activated. Please indicate the configuration used for this test
Detailed results of SETTE tests (restartability and reproducibility for each of the reference configuration)
Results of the required bit comparability tests been run: Are there no differences when activating the development?
If some differences appear, is reason for the change valid/understood?
If some differences appear, is the !ticket describing in detail the impact this change will have on model configurations?
Is this change expected to preserve all diagnostics?
If no, is reason for the change valid/understood?
Are there significant changes in run time/memory?
… … …

Updated on 12/03/2019 16:40:57 by davestorkey

Review

A successful review is needed to schedule the merge of this development into the future NEMO release during next Merge Party (usually in November).

Code changes and documentation

Question Answer Comment
Is the proposed methodology now implemented?
Are the code changes in agreement with the flowchart defined at Preview step?
Are the code changes in agreement with list of routines and variables as proposed at Preview step?
If not, are the discrepancies acceptable?
Is the in-line documentation accurate and sufficient?
Do the code changes comply with NEMO coding standards?
Is the !ticket of development documented with sufficient details for others to understand the impact of the change?
Are the reference manual tex files now updated following the proposed summary in preview section?
Is there a need for some documentation on the web pages (in addition to in-line and reference manual)?
If yes, please describe and ask PI. A yes answer must include all documentation available.
… … …

Review Summary

Is the review fully successful?

Updated on 08/04/2020 02:31:13 by anonymous

Once review is successful, the development must be scheduled for merge during next Merge Party Meeting.

Attachments (1)

Download all attachments as: .zip