= KERNEL-02_Storkey_Coward_IMMERSE_first_steps The PI is responsible to closely follow the progress of the action, and especially to contact NEMO project manager if the delay on preview (or review) are longer than the 2 weeks expected. == Help The action has to be detailed briefly in the 'Summary' for later inclusion in other pages. To do so, '''edit 'Summary' section as a common wiki page and set links for the development ticket and branch'''. '''Out of this, the rest of the page ('Abstract'|'Preview'|'Tests'|'Review') can be edited on-line inside the form fields''' considering the following color code given hereafter: [[span(PI(S), style=background-color:lightgrey)]], [[span(Previewer(s), style=background-color:lightblue)]] and [[span(Reviewer(s), style=background-color:lightgreen)]].\\ 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 #summary ||= '''Action''' =|| KERNEL-02_Storkey_Coward_IMMERSE_first_steps || ||= '''PI(S)''' =|| Dave Storkey, Andrew Coward || {{{#!th '''Digest''' }}} {{{#!td Reorganisation of code to prepare for implementation of RK3 timestepping and tiling (IMMERSE WP3 and WP4) 1. State variables to be renamed with an extra time dimension and time level swapping achieved by swapping time level indices. 2. Prognostic fields to be passed to tendency routines from stp as arguments. 3. Introduction of DO loop macros in preparation for the implementation of tiling. }}} |- ||= '''Dependencies''' =|| || ||= '''Target''' =|| || ||= '''Trac Ticket''' =|| #2258 || ||= '''SVN branch''' =|| [source:/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps] || ||= '''Previewer(s)''' =|| Gurvan Madec || ||= '''Reviewer(s)''' =|| Gurvan Madec || ||= '''Link''' =|| [[ExtractUrl(.)]] || {{{#!comment DON'T REMOVE THIS VOID HEADING OR IT WILL BREAK THE INCLUDE FEATURE BETWEEN WIKI PAGES }}} == == {{{#!comment DON'T REMOVE THIS VOID HEADING OR IT WILL BREAK THE INCLUDE FEATURE BETWEEN WIKI PAGES }}} {{{#!Fold title=Abstract tag=h2 [=#abstract] This section should be completed before starting to develop the code, in order to find agreement on the method beforehand. {{{#!TracForm #!subcontext abstract #!submit_label 'Save Abstract' #!keep_history yes === Description [tf.textarea:description -id=piform -class=taform 'Describe the goal of development, and the methodology.\n\nAdd reference documents or publications if relevant.' 0 20] === Implementation [tf.textarea:implementation -id=piform -class=taform 'Describe flow chart of the changes in the code.\n\nList the .F90 files and modules to be changed.\n\nDetailed list of new variables (including namelists) to be defined. Give for each the chosen name (following coding rules) and definition.' 0 20] === Reference manual and web pages updates [tf.textarea:manual -id=piform -class=taform 'Using part 1 and 2, define the summary of changes to be done in the NEMO reference manual (tex files), and in the content of web pages.' 0 20] '''''Updated on [tf.form_updated_on:] by [tf.form_updater:]''''' }}} 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. }}} {{{#!Fold title=Preview tag=h2 [=#preview] '''1. Variable names and extra dimensions'''\\ [#step1 The refactoring script] \\ [#step2 The refactoring script explained]\\ [#step3 Notes on testing regular expressions]\\ [#step4 Some contrived tests]\\ [#step5 Results on real files]\\ \\ '''2. Do loop macro substitution'''\\ [#step6 Automating the tiling changes]\\ [#step7 Results on traldf_iso.F90]\\ [#step8 do2dfinder.pl]\\ [#step9 do3dfinder.pl]\\ [#step10 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: {{{#!bash #!/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 }}} [=#step1 '''The refactoring script'''] {{{#!bash #!/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 }}} [=#step2 '''The refactoring script explained'''] {{{#!bash # 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 }}} [=#step3 '''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: [[Image(regex101_example_sm.png)]] [=#step4 ''' Some contrived tests:'''] {{{#!f !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:''' {{{#!f 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 }}} {{{#!diff Index: TEST_FILES/contrived_tests.F90 ============================== --- TEST_FILES_ORG/contrived_tests.F90 2019-02-12 17:16:52.000000000 +0000 +++ TEST_FILES/contrived_tests.F90 2019-02-12 17:17:06.000000000 +0000 @@ -1,14 +1,14 @@ - 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?] + 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 }}} 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: [=#step5 '''The results on the sample set of files (patch.list):'''] {{{#!diff Index: TEST_FILES/dynadv_cen2.F90 ============================== --- TEST_FILES_ORG/dynadv_cen2.F90 2019-02-08 10:51:12.000000000 +0000 +++ TEST_FILES/dynadv_cen2.F90 2019-02-12 17:17:06.000000000 +0000 @@ -66,14 +66,14 @@ ! !== Horizontal advection ==! ! DO jk = 1, jpkm1 ! horizontal transport - zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) - zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) + zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) + zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) DO jj = 1, jpjm1 ! horizontal momentum fluxes (at T- and F-point) DO ji = 1, fs_jpim1 ! vector opt. - zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) - zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji ,jj+1,jk) ) - zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) ) - zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) + 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) ) + 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) ) + 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) ) + 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) ) END DO END DO DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes @@ -105,21 +105,21 @@ IF( ln_linssh ) THEN ! linear free surface: advection through the surface DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 - 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) - 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) + 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) + 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) END DO END DO ENDIF DO jk = 2, jpkm1 ! interior advective fluxes DO jj = 2, jpj ! 1/4 * Vertical transport DO ji = 2, jpi - zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) + zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk,Nii) END DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. - 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) ) + 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) ) + 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) ) END DO END DO END DO }}} {{{#!diff Index: TEST_FILES/dynadv_ubs.F90 ============================== --- TEST_FILES_ORG/dynadv_ubs.F90 2019-02-08 10:51:12.000000000 +0000 +++ TEST_FILES/dynadv_ubs.F90 2019-02-12 17:17:06.000000000 +0000 @@ -101,17 +101,17 @@ DO jk = 1, jpkm1 ! Laplacian of the velocity ! ! ! =========================== ! ! ! horizontal volume fluxes - zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) - zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) + zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) + zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) ! DO jj = 2, jpjm1 ! laplacian DO ji = fs_2, fs_jpim1 ! vector opt. - 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) - 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) - zlu_uv(ji,jj,jk,1) = ( ub (ji ,jj+1,jk) - ub (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & - & - ( ub (ji ,jj ,jk) - ub (ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk) - zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj ,jk) - vb (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & - & - ( vb (ji ,jj ,jk) - vb (ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk) + 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) + 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) + zlu_uv(ji,jj,jk,1) = ( uu (ji ,jj+1,jk,Nnn) - uu (ji ,jj ,jk,Nnn) ) * fmask(ji ,jj ,jk) & + & - ( uu (ji ,jj ,jk,Nnn) - uu (ji ,jj-1,jk,Nnn) ) * fmask(ji ,jj-1,jk) + zlv_vu(ji,jj,jk,1) = ( vv (ji+1,jj ,jk,Nnn) - vv (ji ,jj ,jk,Nnn) ) * fmask(ji ,jj ,jk) & + & - ( vv (ji ,jj ,jk,Nnn) - vv (ji-1,jj ,jk,Nnn) ) * fmask(ji-1,jj ,jk) ! 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) 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,13 +131,13 @@ ! ! Horizontal advection ! DO jk = 1, jpkm1 ! ====================== ! ! ! horizontal volume fluxes - zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) - zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) + zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * uu(:,:,jk,Nii) + zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vv(:,:,jk,Nii) ! DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point DO ji = 1, fs_jpim1 ! vector opt. - zui = ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) - zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) + zui = ( uu(ji,jj,jk,Nii) + uu(ji+1,jj ,jk,Nii) ) + zvj = ( vv(ji,jj,jk,Nii) + vv(ji ,jj+1,jk,Nii) ) ! IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) @@ -163,9 +163,9 @@ ENDIF ! zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) & - & * ( un(ji,jj,jk) + un(ji ,jj+1,jk) - gamma1 * zl_u ) + & * ( uu(ji,jj,jk,Nii) + uu(ji ,jj+1,jk,Nii) - gamma1 * zl_u ) zfu_f(ji ,jj ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji ,jj+1,jk,2) ) ) & - & * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) - gamma1 * zl_v ) + & * ( vv(ji,jj,jk,Nii) + vv(ji+1,jj ,jk,Nii) - gamma1 * zl_v ) END DO END DO DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes @@ -198,21 +198,21 @@ IF( ln_linssh ) THEN ! constant volume : advection through the surface DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 - 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) - 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) + 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) + 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) END DO END DO ENDIF DO jk = 2, jpkm1 ! interior fluxes DO jj = 2, jpj DO ji = 2, jpi - zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) + zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk,Nii) END DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. - 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) ) + 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) ) + 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) ) END DO END DO END DO }}} {{{#!diff Index: TEST_FILES/dynkeg.F90 ============================== --- TEST_FILES_ORG/dynkeg.F90 2019-02-08 10:51:12.000000000 +0000 +++ TEST_FILES/dynkeg.F90 2019-02-12 17:17:06.000000000 +0000 @@ -56,8 +56,8 @@ !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] !! * kscheme = nkeg_HW : Hollingsworth correction following !! Arakawa (2001). The now horizontal kinetic energy is given by: - !! zhke = 1/6 [ mi-1( 2 * un^2 + ((un(j+1)+un(j-1))/2)^2 ) - !! + mj-1( 2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2 ) ] + !! zhke = 1/6 [ mi-1( 2 * un^2 + ((uu(j+1,Nii)+uu(j-1,Nii))/2)^2 ) + !! + mj-1( 2 * vn^2 + ((vv(i+1,Nii)+vv(i-1,Nii))/2)^2 ) ] !! !! Take its horizontal gradient and add it to the general momentum !! trend (ua,va). @@ -108,7 +108,7 @@ ii = idx_bdy(ib_bdy)%nbi(jb,igrd) ij = idx_bdy(ib_bdy)%nbj(jb,igrd) ifu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) - un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) + uu(ii-ifu,ij,jk,Nii) = uu(ii,ij,jk,Nii) * umask(ii,ij,jk) END DO END DO ! @@ -118,7 +118,7 @@ ii = idx_bdy(ib_bdy)%nbi(jb,igrd) ij = idx_bdy(ib_bdy)%nbj(jb,igrd) ifv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) - vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) + vv(ii,ij-ifv,jk,Nii) = vv(ii,ij,jk,Nii) * vmask(ii,ij,jk) END DO END DO ENDIF @@ -131,10 +131,10 @@ DO jk = 1, jpkm1 DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. - zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & - & + un(ji ,jj ,jk) * un(ji ,jj ,jk) - zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & - & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) + zu = uu(ji-1,jj ,jk,Nii) * uu(ji-1,jj ,jk,Nii) & + & + uu(ji ,jj ,jk,Nii) * uu(ji ,jj ,jk,Nii) + zv = vv(ji ,jj-1,jk,Nii) * vv(ji ,jj-1,jk,Nii) & + & + vv(ji ,jj ,jk,Nii) * vv(ji ,jj ,jk,Nii) zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) END DO END DO @@ -144,15 +144,15 @@ DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, jpim1 ! vector opt. - zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & - & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & - & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & - & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) + zu = 8._wp * ( uu(ji-1,jj ,jk,Nii) * uu(ji-1,jj ,jk,Nii) & + & + uu(ji ,jj ,jk,Nii) * uu(ji ,jj ,jk,Nii) ) & + & + ( 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) ) & + & + ( uu(ji ,jj-1,jk,Nii) + uu(ji ,jj+1,jk,Nii) ) * ( uu(ji ,jj-1,jk,Nii) + uu(ji ,jj+1,jk,Nii) ) ! - zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & - & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & - & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & - & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) + zv = 8._wp * ( vv(ji ,jj-1,jk,Nii) * vv(ji ,jj-1,jk,Nii) & + & + vv(ji ,jj ,jk,Nii) * vv(ji ,jj ,jk,Nii) ) & + & + ( 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) ) & + & + ( vv(ji-1,jj ,jk,Nii) + vv(ji+1,jj ,jk,Nii) ) * ( vv(ji-1,jj ,jk,Nii) + vv(ji+1,jj ,jk,Nii) ) zhke(ji,jj,jk) = r1_48 * ( zv + zu ) END DO END DO @@ -163,8 +163,8 @@ IF (ln_bdy) THEN ! restore velocity masks at points outside boundary - un(:,:,:) = un(:,:,:) * umask(:,:,:) - vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) + uu(:,:,:,Nii) = uu(:,:,:,Nii) * umask(:,:,:) + vv(:,:,:,Nii) = vv(:,:,:,Nii) * vmask(:,:,:) ENDIF ! }}} {{{#!diff Index: TEST_FILES/dynvor.F90 ============================== }}} {{{#!diff Index: TEST_FILES/flo_oce.F90 ============================== --- TEST_FILES_ORG/flo_oce.F90 2019-02-08 10:51:12.000000000 +0000 +++ TEST_FILES/flo_oce.F90 2019-02-12 17:17:06.000000000 +0000 @@ -59,7 +59,7 @@ !!---------------------------------------------------------------------- !! *** FUNCTION flo_oce_alloc *** !!---------------------------------------------------------------------- - ALLOCATE( wb(jpi,jpj,jpk) , nfloat(jpnfl) , nisobfl(jpnfl) , ngrpfl(jpnfl) , & + ALLOCATE( ww(jpi,jpj,jpk,Nnn) , nfloat(jpnfl) , nisobfl(jpnfl) , ngrpfl(jpnfl) , & & flxx(jpnfl) , flyy(jpnfl) , flzz(jpnfl) , & & tpifl(jpnfl) , tpjfl(jpnfl) , tpkfl(jpnfl) , STAT=flo_oce_alloc ) ! }}} ''' Ok This was wrong.Nnn is not what should go into the ALLOCATE statement''' {{{#!diff Index: TEST_FILES/floats.F90 ============================== --- TEST_FILES_ORG/floats.F90 2019-02-08 10:51:12.000000000 +0000 +++ TEST_FILES/floats.F90 2019-02-12 17:17:06.000000000 +0000 @@ -64,7 +64,7 @@ ! CALL flo_rst( kt ) ! trajectories restart ! - wb(:,:,:) = wn(:,:,:) ! Save the old vertical velocity field + ww(:,:,:,Nnn) = ww(:,:,:,Nii) ! Save the old vertical velocity field ! IF( ln_timing ) CALL timing_stop('flo_stp') ! @@ -131,7 +131,7 @@ ! CALL flo_dom ! compute/read initial position of floats ! - wb(:,:,:) = wn(:,:,:) ! set wb for computation of floats trajectories at the first time step + ww(:,:,:,Nnn) = ww(:,:,:,Nii) ! set wb for computation of floats trajectories at the first time step ! END SUBROUTINE flo_init }}} {{{#!diff Index: TEST_FILES/sbcfwb.F90 ============================== }}} {{{#!diff Index: TEST_FILES/trabbl.F90 ============================== --- TEST_FILES_ORG/trabbl.F90 2019-02-13 12:01:47.000000000 +0000 +++ TEST_FILES/trabbl.F90 2019-02-13 12:01:55.000000000 +0000 @@ -347,8 +347,8 @@ zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) ! zdep(ji,jj) = gdept_n(ji,jj,ik) ! bottom T-level reference depth - zub (ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity - zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) + zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Nii) ! bottom velocity + zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Nii) END DO END DO ! }}} So far so good.... == [=#step6 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: {{{#!f !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: {{{#!diff --- TESTDO_FILES/testdo.F90 2019-03-04 13:43:28.000000000 +0000 +++ TESTDO_FILES_3D/testdo.F90 2019-03-04 13:43:38.000000000 +0000 @@ -1,11 +1,9 @@ ! 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 + DO_2D_00_00 + some loop content + more loop content + END_2D ! followed by an invalid loop pair DO jj = 2,jpjm1 j = jj-1 @@ -16,21 +14,15 @@ 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 + DO_2D_10_11 + some loop content + do jn = 1, jptrc + yet more loop content + end do + more loop content + END_2D ! 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 + DO_3D_10_11( 2,jpkm1 ) + some loop content + more loop content + END_3D }}} [=#step7 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 }}} {{{#!diff --- TESTDO_FILES/traldf_iso.F90 2019-03-04 12:59:44.000000000 +0000 +++ TESTDO_FILES_3D/traldf_iso.F90 2019-03-04 13:26:08.000000000 +0000 @@ -143,58 +143,42 @@ ! IF( kpass == 1 ) THEN !== first pass only ==! ! - DO jk = 2, jpkm1 - DO jj = 2, jpjm1 - DO ji = fs_2, fs_jpim1 ! vector opt. - ! - zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & - & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) - zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & - & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) - ! - zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & - & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku - zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & - & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv - ! - ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & - & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) - END DO - END DO - END DO + DO_3D_00_00( 2, jpkm1 ) + ! + zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & + & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) + zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & + & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) + ! + zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & + & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku + zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & + & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv + ! + ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & + & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) + END_3D ! IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient - DO jk = 2, jpkm1 - DO jj = 2, jpjm1 - DO ji = fs_2, fs_jpim1 - akz(ji,jj,jk) = 0.25_wp * ( & - & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & - & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & - & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & - & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) - END DO - END DO - END DO + DO_3D_00_00( 2, jpkm1 ) + akz(ji,jj,jk) = 0.25_wp * ( & + & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & + & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & + & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & + & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) + END_3D ! IF( ln_traldf_blp ) THEN ! bilaplacian operator - DO jk = 2, jpkm1 - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 - akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & - & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) - END DO - END DO - END DO + DO_3D_10_10( 2, jpkm1 ) + akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & + & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) + END_3D ELSEIF( ln_traldf_lap ) THEN ! laplacian operator - DO jk = 2, jpkm1 - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 - ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) - zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) - akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt - END DO - END DO - END DO + DO_3D_10_10( 2, jpkm1 ) + ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) + zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) + akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt + END_3D ENDIF ! ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit @@ -215,28 +199,20 @@ !!end ! Horizontal tracer gradient - DO jk = 1, jpkm1 - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 ! vector opt. - zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) - zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) - END DO - END DO - END DO + DO_3D_10_10( 1, jpkm1 ) + zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) + zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) + END_3D IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient - DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) - DO ji = 1, fs_jpim1 ! vector opt. - zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) - zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) - END DO - END DO + DO_2D_10_10 + zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) + zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) + END_2D IF( ln_isfcav ) THEN ! first wet level beneath a cavity - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 ! vector opt. - IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) - IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) - END DO - END DO + DO_2D_10_10 + IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) + IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) + END_2D ENDIF ENDIF ! @@ -252,36 +228,32 @@ IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) ENDIF - DO jj = 1 , jpjm1 !== Horizontal fluxes - DO ji = 1, fs_jpim1 ! vector opt. - zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) - zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) - ! - zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & - & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) - ! - zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & - & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) - ! - zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku - zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv - ! - zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & - & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & - & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) - zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & - & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & - & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) - END DO - END DO - ! - DO jj = 2 , jpjm1 !== horizontal divergence and add to pta - DO ji = fs_2, fs_jpim1 ! vector opt. - pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & - & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & - & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) - END DO - END DO + DO_2D_10_10 + zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) + zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) + ! + zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & + & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) + ! + zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & + & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) + ! + zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku + zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv + ! + zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & + & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & + & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) + zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & + & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & + & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) + END_2D + ! + DO_2D_00_00 + pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & + & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & + & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) + END_2D END DO ! End of slab !!---------------------------------------------------------------------- @@ -295,75 +267,55 @@ ! ! Surface and bottom vertical fluxes set to zero ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp - DO jk = 2, jpkm1 ! interior (2=) { 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 = || 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 = || 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 = || 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 $_; } } }}} [=#step9 and] {{{#!perl #cat do3dfinder.pl # use IO::Scalar; open(F,$ARGV[0]) || die "Cannot open $ARGV[0]: $!"; while() { 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 = || 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 = || 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 = || 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 $_; } } }}} == [=#step10 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: {{{#!diff --- TESTDO_FILES/traldf_iso.F90 2019-03-04 12:59:44.000000000 +0000 +++ SANITY/traldf_iso.f90 2019-03-04 17:17:32.000000000 +0000 @@ -1,3 +1,6 @@ + + + MODULE traldf_iso !!====================================================================== !! *** MODULE traldf_iso *** @@ -39,7 +42,17 @@ LOGICAL :: l_hst ! flag to compute heat transport !! * Substitutions -# include "vectopt_loop_substitute.h90" + !!---------------------------------------------------------------------- + !! *** vectopt_loop_substitute *** + !!---------------------------------------------------------------------- + !! ** purpose : substitute the inner loop start/end indices with CPP macro + !! allow unrolling of do-loop (useful with vector processors) + !!---------------------------------------------------------------------- + !!---------------------------------------------------------------------- + !! NEMO/OCE 4.0 , NEMO Consortium (2018) + !! $Id: vectopt_loop_substitute.h90 10068 2018-08-28 14:09:04Z nicolasmartin $ + !! Software governed by the CeCILL license (see ./LICENSE) + !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id: traldf_iso.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ @@ -143,58 +156,42 @@ ! IF( kpass == 1 ) THEN !== first pass only ==! ! - DO jk = 2, jpkm1 - DO jj = 2, jpjm1 - DO ji = fs_2, fs_jpim1 ! vector opt. - ! - zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & - & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) - zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & - & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) - ! - zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & - & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku - zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & - & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv - ! - ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & - & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) - END DO - END DO - END DO + DO jk = 2, jpkm1 ; DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 + ! + zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & + & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) + zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & + & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) + ! + zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & + & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku + zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & + & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv + ! + ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & + & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) + END DO ; END DO ; END DO ! IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient - DO jk = 2, jpkm1 - DO jj = 2, jpjm1 - DO ji = fs_2, fs_jpim1 - akz(ji,jj,jk) = 0.25_wp * ( & - & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & - & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & - & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & - & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) - END DO - END DO - END DO + DO jk = 2, jpkm1 ; DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 + akz(ji,jj,jk) = 0.25_wp * ( & + & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & + & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & + & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & + & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) + END DO ; END DO ; END DO ! IF( ln_traldf_blp ) THEN ! bilaplacian operator - DO jk = 2, jpkm1 - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 - akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & - & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) - END DO - END DO - END DO + DO jk = 2, jpkm1 ; DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 + akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & + & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) + END DO ; END DO ; END DO ELSEIF( ln_traldf_lap ) THEN ! laplacian operator - DO jk = 2, jpkm1 - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 - ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) - zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) - akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt - END DO - END DO - END DO + DO jk = 2, jpkm1 ; DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 + ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) + zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) + akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt + END DO ; END DO ; END DO ENDIF ! ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit @@ -215,28 +212,20 @@ !!end ! Horizontal tracer gradient - DO jk = 1, jpkm1 - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 ! vector opt. - zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) - zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) - END DO - END DO - END DO + DO jk = 1, jpkm1 ; DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 + zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) + zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) + END DO ; END DO ; END DO IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient - DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) - DO ji = 1, fs_jpim1 ! vector opt. - zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) - zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) - END DO - END DO + DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 + zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) + zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) + END DO ; END DO IF( ln_isfcav ) THEN ! first wet level beneath a cavity - DO jj = 1, jpjm1 - DO ji = 1, fs_jpim1 ! vector opt. - IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) - IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) - END DO - END DO + DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 + IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) + IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) + END DO ; END DO ENDIF ENDIF ! @@ -252,36 +241,32 @@ IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) ENDIF - DO jj = 1 , jpjm1 !== Horizontal fluxes - DO ji = 1, fs_jpim1 ! vector opt. - zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) - zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) - ! - zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & - & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) - ! - zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & - & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) - ! - zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku - zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv - ! - zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & - & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & - & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) - zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & - & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & - & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) - END DO - END DO + DO jj = 2-1, jpjm1 ; DO ji = 2-1, jpim1 + zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) + zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) + ! + zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & + & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) + ! + zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & + & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) + ! + zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku + zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv + ! + zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & + & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & + & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) + zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & + & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & + & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) + END DO ; END DO ! - DO jj = 2 , jpjm1 !== horizontal divergence and add to pta - DO ji = fs_2, fs_jpim1 ! vector opt. - pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & - & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & - & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) - END DO - END DO + DO jj = 2 ,jpjm1 ; DO ji = 2 ,jpim1 + pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & + & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & + & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) + END DO ; END DO END DO ! End of slab !!---------------------------------------------------------------------- @@ -295,75 +280,55 @@ ! ! Surface and bottom vertical fluxes set to zero ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp - DO jk = 2, jpkm1 ! interior (2=