Version 12 (modified by acc, 18 months ago) (diff)

KERNEL-02_DaveStorkey_RK3Preparation

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_DaveStorkey_RK3Preparation
PI(S) Dave Storkey

Digest

Reorganisation of code to prepare for implementation of RK3 timestepping:

  1. Prognostic fields to be passed to tendency routines as arguments.
  2. Tracers to be consistency updated as extensive (volume-weighted) quantities throughout the code.
Dependencies
Target
Trac Ticket #0000
SVN branch branches/$YEAR/dev_r{REV}_{SCOPE}-{NUM}_{PIS}-{KEYWORDS}
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

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

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+$//;
         #
         # 5. Store the leading indentation for the outer loop
         #
         ($jspac = $jline) =~ s/(^[\s]*)([^\s]*).*/\1/;
         chomp($jspac);
         #
         # 6. Construct a DO_2D line to replace the original statements
         #
         print $jspac,"DO_2D( ",$iargs," , ",$jargs," )\n";
         #
         # 7. 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?" ;
            #
            # 8. Increment a counter if another DO statement is found
            #
            if ( $iline =~ /^\s*do\s*/i )  { $docount++ };
            #
            # 9. Decrement a counter if a END DO statement is found
            #
            if ( $iline =~ /^\s*end\s*do/i )  { $docount-- };
            #
            # 10. A negative counter means the matching END DO for the ji loop has been reached
            #
            if ( $docount < 0 ) {
               #
               # 11. 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 {
               #
               # 12. This is a line inside the loop. Remove three leading spaces (if any) and output.
               #
               $iline =~ s/^\s\s\s//;
               print $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 $_;
   }
}

Secondly, here is a simple test file:

!   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 an valid loop pair with a nested do
     DO jj  = 2,jpjm1
        DO ji = 2,jpim1
           some loop content
           do jn = 1, jptrc
              yet more loop content
           end do
           more loop content
        END DO
     END DO

Followed by the result of running: perl do2dfinder.pl testdo.F90 > testdo_after.F90 and the difference:

!   some random text
!   followed by a valid loop pair
     DO_2D( 2,jpim1 , 2,jpjm1 )
        some loop content
        more loop content
     END_2D
!   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 an valid loop pair with a nested do
     DO_2D( 2,jpim1 , 2,jpjm1 )
        some loop content
        do jn = 1, jptrc
           yet more loop content
        end do
        more loop content
     END_2D
  • .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( 2,jpim1 , 2,jpjm1 ) 
     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 an valid loop pair with a nested do 
    19      DO jj  = 2,jpjm1 
    20         DO ji = 2,jpim1 
    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( 2,jpim1 , 2,jpjm1 ) 
     18        some loop content 
     19        do jn = 1, jptrc 
     20           yet more loop content 
     21        end do 
     22        more loop content 
     23     END_2D 

And finally, the results on a real case dynvor.F90 (just a sample of the differences):

  • .F90

    old new  
    225225      SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==! 
    226226      CASE ( np_RVO )                           !* relative vorticity 
    227227         DO jk = 1, jpkm1                                 ! Horizontal slab 
    228             DO jj = 1, jpjm1 
    229                DO ji = 1, jpim1 
    230                   zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
    231                      &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    232                END DO 
    233             END DO 
     228            DO_2D( 1, jpim1 , 1, jpjm1 ) 
     229               zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
     230                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     231            END_2D 
    234232            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    235                DO jj = 1, jpjm1 
    236                   DO ji = 1, jpim1 
    237                      zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    238                   END DO 
    239                END DO 
     233               DO_2D( 1, jpim1 , 1, jpjm1 ) 
     234                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     235               END_2D 
    240236            ENDIF 
    241237         END DO 
    242238 
     
    244240 
    245241      CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    246242         DO jk = 1, jpkm1                                 ! Horizontal slab 
    247             DO jj = 1, jpjm1 
    248                DO ji = 1, jpim1                          ! relative vorticity 
    249                   zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
    250                      &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    251                END DO 
    252             END DO 
     243            DO_2D( 1, jpim1 , 1, jpjm1 ) 
     244               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   & 
     245                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
     246            END_2D 
    253247            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    254                DO jj = 1, jpjm1 
    255                   DO ji = 1, jpim1 
    256                      zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    257                   END DO 
    258                END DO 
     248               DO_2D( 1, jpim1 , 1, jpjm1 ) 
     249                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
     250               END_2D 
    259251            ENDIF 
    260252         END DO 
    261253 
     
    271263         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
    272264            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 
    273265         CASE ( np_RVO )                           !* relative vorticity 
    274             DO jj = 2, jpj 
    275                DO ji = 2, jpi   ! vector opt. 
    276                   zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    277                      &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
    278                END DO 
    279             END DO 
     266            DO_2D( 2, jpi , 2, jpj ) 
     267               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
     268                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     269            END_2D 
    280270         CASE ( np_MET )                           !* metric term 
    281             DO jj = 2, jpj 
    282                DO ji = 2, jpi 
    283                   zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    284                      &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
    285                END DO 
    286             END DO 
     271            DO_2D( 2, jpi , 2, jpj ) 
     272               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
     273                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk) 
     274            END_2D 
    287275         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    288             DO jj = 2, jpj 
    289                DO ji = 2, jpi   ! vector opt. 
    290                   zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    291                      &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
    292                END DO 
    293             END DO 
     276            DO_2D( 2, jpi , 2, jpj ) 
     277               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
     278                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     279            END_2D 
    294280         CASE ( np_CME )                           !* Coriolis + metric 
    295             DO jj = 2, jpj 
    296                DO ji = 2, jpi   ! vector opt. 
    297                   zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    298                        &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    299                        &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
    300                END DO 
    301             END DO 
     281            DO_2D( 2, jpi , 2, jpj ) 
     282               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
     283                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
     284                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk) 
     285            END_2D 
    302286         CASE DEFAULT                                             ! error 
    303287            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    304288         END SELECT 
    305289         ! 
    306290         !                                   !==  compute and add the vorticity term trend  =! 
    307          DO jj = 2, jpjm1 
    308             DO ji = 2, jpim1   ! vector opt. 
    309                pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
    310                   &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
    311                   &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
    312                   ! 
    313                pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
    314                   &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
    315                   &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
    316             END DO 
    317          END DO 
     291         DO_2D( 2, jpim1 , 2, jpjm1 ) 
     292            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    & 
     293               &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   & 
     294               &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   ) 
     295               ! 
     296            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    & 
     297               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
     298               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
     299         END_2D 
    318300         !                                             ! =============== 
    319301      END DO                                           !   End of slab 
    320302      !                                                ! =============== 

Extending for the 3D loops

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:

      DO jk = 2, jpkm1                    ! interior advective fluxes
         DO_2D( 2, jpi , 2, jpj )
            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
         END_2D
         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 )
            zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
            zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
         END_2D
      END DO
      DO jk = 1, jpkm1                    ! divergence of vertical momentum flux divergence
         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 )
            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)
            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)
         END_2D
      END DO

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:

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\s*\(/i) {
         my $isinavlid = 0;
         #
         # 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+$//;
         #
         # 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;
         print $altr $jspac,"DO_3D( ",$iargs," , ",$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 $_;
   }
}

These scripts can be run sequentially; for example:

  cp TEST_FILES_ORG/dynadv_cen2.F90 TESTDO_FILES/
  perl do2dfinder.pl TESTDO_FILES/dynadv_cen2.F90 > TESTDO_FILES_2D/dynadv_cen2.F90
  perl do3dfinder.pl TESTDO_FILES_2D/dynadv_cen2.F90 > TESTDO_FILES_3D/dynadv_cen2.F90

And the difference between the original and final version is:

  • dynadv_cen2.F90

    old new  
    6868      DO jk = 1, jpkm1                    ! horizontal transport 
    6969         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    7070         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    71          DO jj = 1, jpjm1                 ! horizontal momentum fluxes (at T- and F-point) 
    72             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) ) 
    77             END DO 
    78          END DO 
    79          DO jj = 2, jpjm1                 ! divergence of horizontal momentum fluxes 
    80             DO ji = fs_2, fs_jpim1   ! vector opt. 
    81                ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
    82                   &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    83                va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
    84                   &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
    85             END DO 
    86          END DO 
     71         DO_2D( 1, fs_jpim1 , 1, jpjm1 ) 
     72            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) ) 
     73            zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) ) 
     74            zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) ) 
     75            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) ) 
     76         END_2D 
     77         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     78            ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
     79               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     80            va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
     81               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     82         END_2D 
    8783      END DO 
    8884      ! 
    8985      IF( l_trddyn ) THEN           ! trends: send trend to trddyn for diagnostic 
     
    9692      ! 
    9793      !                             !==  Vertical advection  ==! 
    9894      ! 
    99       DO jj = 2, jpjm1                    ! surface/bottom advective fluxes set to zero 
    100          DO ji = fs_2, fs_jpim1 
    101             zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp 
    102             zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp 
    103          END DO 
    104       END DO 
     95      DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     96         zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp 
     97         zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp 
     98      END_2D 
    10599      IF( ln_linssh ) THEN                ! linear free surface: advection through the surface 
    106          DO jj = 2, jpjm1 
    107             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) 
    110             END DO 
    111          END DO 
     100         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     101            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) 
     102            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) 
     103         END_2D 
    112104      ENDIF 
    113105      DO jk = 2, jpkm1                    ! interior advective fluxes 
    114          DO jj = 2, jpj                       ! 1/4 * Vertical transport 
    115             DO ji = 2, jpi 
    116                zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
    117             END DO 
    118          END DO 
    119          DO jj = 2, jpjm1 
    120             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) ) 
    123             END DO 
    124          END DO 
    125       END DO 
    126       DO jk = 1, jpkm1                    ! divergence of vertical momentum flux divergence 
    127          DO jj = 2, jpjm1 
    128             DO ji = fs_2, fs_jpim1   ! vector opt. 
    129                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) 
    130                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) 
    131             END DO 
    132          END DO 
     106         DO_2D( 2, jpi , 2, jpj ) 
     107            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     108         END_2D 
     109         DO_2D( fs_2, fs_jpim1 , 2, jpjm1 ) 
     110            zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 
     111            zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 
     112         END_2D 
    133113      END DO 
     114      DO_3D( fs_2, fs_jpim1 , 2, jpjm1 , 1, jpkm1 ) 
     115         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) 
     116         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) 
     117      END_3D 
    134118      ! 
    135119      IF( l_trddyn ) THEN                 ! trends: send trend to trddyn for diagnostic 
    136120         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 

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/09/2020 23:21:05 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/09/2020 23:21:05 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