| Version 19 (modified by davestorkey, 6 years 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_Storkey_Coward_IMMERSE_first_steps |
|---|---|
| PI(S) | Dave Storkey |
Digest | Reorganisation of code to prepare for implementation of RK3 timestepping:
|
| 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 | ExtractUrl(.)? |
Abstract
This section should be completed before starting to develop the code, in order to find agreement on the method beforehand.
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:
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 u n(:,:,:) ! The simplest test. Should ==> uu(:,:,:,jtn)2 u n ( ji , jj, :) ! Check alternative simple case ==> uu ( ji , jj, :,jtn)3 a+u b(:,:,jk) ! Check preceeding operators are correctly recognised4 a-v b(:,:,jk) ! Check preceeding operators are correctly recognised5 a*v n(:,:,jk) ! Check preceeding operators are correctly recognised6 a/w b(:,:,jk) ! Check preceeding operators are correctly recognised7 a%w n(:,:,jk) ! Check preceeding operators are correctly recognised8 .OR.w n(:,:,jk) ! Check preceeding operators are correctly recognised9 (w n(:,:,jk) + wn(:,:,jk-1)) ! Check preceeding brackets are correctly recognised10 u n( ji+1 , & ! Check that entries over11 & jj, jk - 1 ) ! multiple lines are handled12 w n( 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?] 13 13 pun(:,:,:) ! target preceded by non-whitespace or operator. Should be unchanged 14 14 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 66 66 ! !== Horizontal advection ==! 67 67 ! 68 68 DO jk = 1, jpkm1 ! horizontal transport 69 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * u n(:,:,jk)70 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * v n(:,:,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) 71 71 DO jj = 1, jpjm1 ! horizontal momentum fluxes (at T- and F-point) 72 72 DO ji = 1, fs_jpim1 ! vector opt. 73 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( u n(ji,jj,jk) + un(ji+1,jj ,jk) )74 zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( u n(ji,jj,jk) + un(ji ,jj+1,jk) )75 zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( v n(ji,jj,jk) + vn(ji+1,jj ,jk) )76 zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( v n(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) ) 77 77 END DO 78 78 END DO 79 79 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes … … 105 105 IF( ln_linssh ) THEN ! linear free surface: advection through the surface 106 106 DO jj = 2, jpjm1 107 107 DO ji = fs_2, fs_jpim1 108 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * w n(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) * w n(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) 110 110 END DO 111 111 END DO 112 112 ENDIF 113 113 DO jk = 2, jpkm1 ! interior advective fluxes 114 114 DO jj = 2, jpj ! 1/4 * Vertical transport 115 115 DO ji = 2, jpi 116 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * w n(ji,jj,jk)116 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk,Nii) 117 117 END DO 118 118 END DO 119 119 DO jj = 2, jpjm1 120 120 DO ji = fs_2, fs_jpim1 ! vector opt. 121 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( u n(ji,jj,jk) + un(ji,jj,jk-1) )122 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( v n(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) ) 123 123 END DO 124 124 END DO 125 125 END DO
-
dynadv_ubs.F90
==============================
old new 101 101 DO jk = 1, jpkm1 ! Laplacian of the velocity ! 102 102 ! ! =========================== ! 103 103 ! ! horizontal volume fluxes 104 zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * u n(:,:,jk)105 zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * v n(:,:,jk)104 zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) 105 zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) 106 106 ! 107 107 DO jj = 2, jpjm1 ! laplacian 108 108 DO ji = fs_2, fs_jpim1 ! vector opt. 109 zlu_uu(ji,jj,jk,1) = ( u b (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) = ( v b (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) = ( u b (ji ,jj+1,jk) - ub (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &112 & - ( u b (ji ,jj ,jk) - ub (ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk)113 zlv_vu(ji,jj,jk,1) = ( v b (ji+1,jj ,jk) - vb (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &114 & - ( v b (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) 115 115 ! 116 116 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) 117 117 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) … … 131 131 ! ! Horizontal advection ! 132 132 DO jk = 1, jpkm1 ! ====================== ! 133 133 ! ! horizontal volume fluxes 134 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * u n(:,:,jk)135 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * v n(:,:,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) 136 136 ! 137 137 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point 138 138 DO ji = 1, fs_jpim1 ! vector opt. 139 zui = ( u n(ji,jj,jk) + un(ji+1,jj ,jk) )140 zvj = ( v n(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) ) 141 141 ! 142 142 IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) 143 143 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) … … 163 163 ENDIF 164 164 ! 165 165 zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) & 166 & * ( u n(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 ) 167 167 zfu_f(ji ,jj ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji ,jj+1,jk,2) ) ) & 168 & * ( v n(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 ) 169 169 END DO 170 170 END DO 171 171 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes … … 198 198 IF( ln_linssh ) THEN ! constant volume : advection through the surface 199 199 DO jj = 2, jpjm1 200 200 DO ji = fs_2, fs_jpim1 201 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * w n(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) * w n(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) 203 203 END DO 204 204 END DO 205 205 ENDIF 206 206 DO jk = 2, jpkm1 ! interior fluxes 207 207 DO jj = 2, jpj 208 208 DO ji = 2, jpi 209 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * w n(ji,jj,jk)209 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk,Nii) 210 210 END DO 211 211 END DO 212 212 DO jj = 2, jpjm1 213 213 DO ji = fs_2, fs_jpim1 ! vector opt. 214 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( u n(ji,jj,jk) + un(ji,jj,jk-1) )215 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( v n(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) ) 216 216 END DO 217 217 END DO 218 218 END DO
-
dynkeg.F90
==============================
old new 56 56 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 57 57 !! * kscheme = nkeg_HW : Hollingsworth correction following 58 58 !! Arakawa (2001). The now horizontal kinetic energy is given by: 59 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u n(j+1)+un(j-1))/2)^2 )60 !! + mj-1( 2 * vn^2 + ((v n(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 ) ] 61 61 !! 62 62 !! Take its horizontal gradient and add it to the general momentum 63 63 !! trend (ua,va). … … 108 108 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 109 109 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 110 110 ifu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 111 u n(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) 112 112 END DO 113 113 END DO 114 114 ! … … 118 118 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 119 119 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 120 120 ifv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 121 v n(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) 122 122 END DO 123 123 END DO 124 124 ENDIF … … 131 131 DO jk = 1, jpkm1 132 132 DO jj = 2, jpj 133 133 DO ji = fs_2, jpi ! vector opt. 134 zu = u n(ji-1,jj ,jk) * un(ji-1,jj ,jk) &135 & + u n(ji ,jj ,jk) * un(ji ,jj ,jk)136 zv = v n(ji ,jj-1,jk) * vn(ji ,jj-1,jk) &137 & + v n(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) 138 138 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 139 139 END DO 140 140 END DO … … 144 144 DO jk = 1, jpkm1 145 145 DO jj = 2, jpjm1 146 146 DO ji = fs_2, jpim1 ! vector opt. 147 zu = 8._wp * ( u n(ji-1,jj ,jk) * un(ji-1,jj ,jk) &148 & + u n(ji ,jj ,jk) * un(ji ,jj ,jk) ) &149 & + ( u n(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) &150 & + ( u n(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) ) 151 151 ! 152 zv = 8._wp * ( v n(ji ,jj-1,jk) * vn(ji ,jj-1,jk) &153 & + v n(ji ,jj ,jk) * vn(ji ,jj ,jk) ) &154 & + ( v n(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) &155 & + ( v n(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) ) 156 156 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 157 157 END DO 158 158 END DO … … 163 163 164 164 IF (ln_bdy) THEN 165 165 ! restore velocity masks at points outside boundary 166 u n(:,:,:) = un(:,:,:) * umask(:,:,:)167 v n(:,:,:) = vn(:,:,:) * vmask(:,:,:)166 uu(:,:,:,Nii) = uu(:,:,:,Nii) * umask(:,:,:) 167 vv(:,:,:,Nii) = vv(:,:,:,Nii) * vmask(:,:,:) 168 168 ENDIF 169 169 170 170 !
Index: TEST_FILES/dynvor.F90
==============================
-
flo_oce.F90
==============================
old new 59 59 !!---------------------------------------------------------------------- 60 60 !! *** FUNCTION flo_oce_alloc *** 61 61 !!---------------------------------------------------------------------- 62 ALLOCATE( w b(jpi,jpj,jpk) , nfloat(jpnfl) , nisobfl(jpnfl) , ngrpfl(jpnfl) , &62 ALLOCATE( ww(jpi,jpj,jpk,Nnn) , nfloat(jpnfl) , nisobfl(jpnfl) , ngrpfl(jpnfl) , & 63 63 & flxx(jpnfl) , flyy(jpnfl) , flzz(jpnfl) , & 64 64 & tpifl(jpnfl) , tpjfl(jpnfl) , tpkfl(jpnfl) , STAT=flo_oce_alloc ) 65 65 !
Ok This was wrong.Nnn is not what should go into the ALLOCATE statement
-
floats.F90
==============================
old new 64 64 ! 65 65 CALL flo_rst( kt ) ! trajectories restart 66 66 ! 67 w b(:,:,:) = wn(:,:,:) ! Save the old vertical velocity field67 ww(:,:,:,Nnn) = ww(:,:,:,Nii) ! Save the old vertical velocity field 68 68 ! 69 69 IF( ln_timing ) CALL timing_stop('flo_stp') 70 70 ! … … 131 131 ! 132 132 CALL flo_dom ! compute/read initial position of floats 133 133 ! 134 w b(:,:,:) = wn(:,:,:) ! set wb for computation of floats trajectories at the first time step134 ww(:,:,:,Nnn) = ww(:,:,:,Nii) ! set wb for computation of floats trajectories at the first time step 135 135 ! 136 136 END SUBROUTINE flo_init 137 137
Index: TEST_FILES/sbcfwb.F90
==============================
-
trabbl.F90
==============================
old new 347 347 zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 348 348 ! 349 349 zdep(ji,jj) = gdept_n(ji,jj,ik) ! bottom T-level reference depth 350 zub (ji,jj) = u n(ji,jj,mbku(ji,jj)) ! bottom velocity351 zvb (ji,jj) = v n(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) 352 352 END DO 353 353 END DO 354 354 !
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 1 1 ! some random text 2 2 ! 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 9 7 ! followed by an invalid loop pair 10 8 DO jj = 2,jpjm1 11 9 j = jj-1 … … 16 14 END DO 17 15 END DO 18 16 ! 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 28 24 ! 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 143 143 ! 144 144 IF( kpass == 1 ) THEN !== first pass only ==! 145 145 ! 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 165 161 ! 166 162 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 178 170 ! 179 171 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 188 176 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 198 182 ENDIF 199 183 ! 200 184 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit … … 215 199 !!end 216 200 217 201 ! 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 226 206 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 233 211 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 240 216 ENDIF 241 217 ENDIF 242 218 ! … … 252 228 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 253 229 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 254 230 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 285 257 END DO ! End of slab 286 258 287 259 !!---------------------------------------------------------------------- … … 295 267 ! ! Surface and bottom vertical fluxes set to zero 296 268 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 297 269 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 322 290 ! !== add the vertical 33 flux ==! 323 291 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 333 297 ! 334 298 ELSE ! bilaplacian 335 299 SELECT CASE( kpass ) 336 300 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 346 306 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 356 312 END SELECT 357 313 ENDIF 358 314 ! 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 367 319 ! 368 320 IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==! 369 321 ( 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 1 4 MODULE traldf_iso 2 5 !!====================================================================== 3 6 !! *** MODULE traldf_iso *** … … 39 42 LOGICAL :: l_hst ! flag to compute heat transport 40 43 41 44 !! * 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 !!---------------------------------------------------------------------- 43 56 !!---------------------------------------------------------------------- 44 57 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 45 58 !! $Id: traldf_iso.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ … … 143 156 ! 144 157 IF( kpass == 1 ) THEN !== first pass only ==! 145 158 ! 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 165 174 ! 166 175 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 178 183 ! 179 184 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 188 189 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 198 195 ENDIF 199 196 ! 200 197 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit … … 215 212 !!end 216 213 217 214 ! 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 226 219 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 233 224 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 240 229 ENDIF 241 230 ENDIF 242 231 ! … … 252 241 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 253 242 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 254 243 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 277 264 ! 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 285 270 END DO ! End of slab 286 271 287 272 !!---------------------------------------------------------------------- … … 295 280 ! ! Surface and bottom vertical fluxes set to zero 296 281 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 297 282 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 322 303 ! !== add the vertical 33 flux ==! 323 304 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 333 310 ! 334 311 ELSE ! bilaplacian 335 312 SELECT CASE( kpass ) 336 313 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 346 319 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 356 325 END SELECT 357 326 ENDIF 358 327 ! 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 367 332 ! 368 333 IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==! 369 334 ( 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.
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.
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).
Once review is successful, the development must be scheduled for merge during next Merge Party Meeting.
Attachments (1)
-
regex101_example_sm.png
(390.0 KB) -
added by acc 6 years ago.
regex tester
Download all attachments as: .zip
