Index: NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/fpp_wrapper_BENCHO2M.sh
===================================================================
--- NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/fpp_wrapper_BENCHO2M.sh	(revision 15762)
+++ NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/fpp_wrapper_BENCHO2M.sh	(revision 15762)
@@ -0,0 +1,188 @@
+#! /bin/bash
+#======================================================================
+#
+# fpp_wrapper.sh <FPP command> <wrapping mode> <FPP arguments>
+#
+# Wrapped FPP command for optional additional PSyclone processing of
+# source code
+#
+#======================================================================
+# History : 4.0 ! 2021/2022 (S. Mueller) Original version
+#======================================================================
+
+#----------------------------------------------------------------------
+# NEMO 4.0 , NEMO Consortium (2022)
+# $Id$ 
+# Software governed by the CeCILL license (see ./LICENSE)
+#----------------------------------------------------------------------
+
+# Path to PSyclone transformation script
+PSyclone_TS=../mk_psyclone/psyclone_transformation_BENCHO2M.py
+
+FPP=$1; shift
+# Extended pre-processing mode
+MODE=$1; shift
+# Original arguments
+FPP_ARG=$@; shift $(($# - 1))
+#
+# PSyclone pre-processing of NEMO source-code files
+if [ ${MODE} == 'PSyclone' ]; then
+    # Filename of file to be processed
+    fname=$1
+    # Call wrapped preprocessor, store in temporary file
+    $FPP -o ${fname}.preprocessed.$$ ${FPP_ARG}
+    # PSyclone pre-processing
+    # PSyclone processing, unless ACC directives present
+    if `grep -q '!\$ACC' ${fname}`; then
+	cat ${fname}.preprocessed.$$
+    else
+        case ${fname} in
+
+	    # Do NOT process some files that do not pass the PSyclone
+	    # pass-through test (see also PSyclone file
+	    # examples/nemo/scripts/process_nemo.sh)
+	    ( obs_utils.[fF]90 | par_kind.[fF]90 | storng.[fF]90 | diurnal_bulk.[fF]90 | timing.[fF]90 | obs_fbm.[fF]90 | sbccpl.[fF]90 ) \
+		    cat ${fname}.preprocessed.$$ ;;
+
+	    # Special case for dynspg_ts.F90: enable inclusion of the intrinsic
+	    # function 'sum' with a 'dim' argument (this relates to PSyclone
+	    # issue #840, where support of the 'dim' argument is discussed in
+	    # the context of the 'maxval'/'minval' intrinsic fuctions); as
+	    # PSyclone accepts 'sum' as a binary function, the named argument
+	    # 'dim' can simply be replaced by its unnamed variant
+	    ( dynspg_ts.[fF]90 ) \
+		    sed -i -e 's/\(.*[sS][uU][mM](.*\)[dD][iI][mM][[:space:]]*=[[:space:]]*\([0-9[:space:]]*\))/\1\2)/' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+	    # Special case for stpctl.F90: enable inclusion of the intrinsic
+	    # function 'MAXVAL' in kernels regions (pending resolution of
+	    # PSyclone issue #840)
+	    ( stpctl.[fF]90 ) \
+		    sed -i -e 's/\([[:space:]]*zmax([[:digit:]]) = \)MAXVAL(\(.*\)mask = \([lw][lm][ma]sk\)/\1MATMUL(\2\3/' -e 's/\([[:space:]]*zmax([[:digit:]]) = \)MAXVAL/\1SIN/' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e 's/MATMUL\((.*\)\([lw][lm][ma]sk\)/MAXVAL\1mask = \2/' -e 's/SIN(/MAXVAL(/' ;;
+
+	    # Special case for dommsk.F90: enable processing of the majority of
+	    # the source code despite the reluctance of PSyclone to convert
+	    # procedures with namelists into PSyIR Invoke objects (reported)
+	    ( dommsk.[fF]90 ) \
+		    sed -i -e '/^904/aCALL dom_msk2\nEND SUBROUTINE dom_msk\nSUBROUTINE dom_msk2' -e '/END SUBROUTINE dom_msk/cEND SUBROUTINE dom_msk2' ${fname}.preprocessed.$$ ; \
+	            sed -i -e '/INTEGER *:: *ji,/,/INTEGER *:: *ios/H' -e '/^ *SUBROUTINE dom_msk2/G' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+	    # Special case for mocsy_mainmod.F90: enable processing of
+	    # subroutine 'constants' by removing use of 'data' statement
+	    # (reported), counteract PSyclone bug #1558, add "!$acc routine
+	    # seq" (pending PSyclone issue #1041)
+	    ( mocsy_mainmod.[fF]90 ) \
+		    sed -i -e '/DIMENSION(12) :: a0/d' -r -e 's/DATA ([^\/]*)\//real(kind=wp),dimension(12) :: \1 = (\//' -e 's/\/$/\/)/' -r -e 's/= \(\/[[:space:]]*12[[:space:]]*\*([^\/]*)\/\)/= \1/' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e 's/ :: onehalf/, parameter :: onehalf=0.5d0/' -e 's/ :: two/, parameter :: two=2.0d0/' -e 's/ :: three/, parameter :: three=3.0d0/' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(.*)/a !$acc routine seq' | \
+		    sed -e '/^\s*[fF][uU][nN][cC][tT][iI][oO][nN].*(.*)/a !$acc routine seq' ;;
+
+	    # Special case for lib_fortran.F90: add "!$acc routine seq"
+	    # directive to SIGN function (pending PSyclone issue #1041);
+	    # temporarily rename some 'CMPLX' intrinsic-function calls to 'SUM'
+	    # calls to facilitate the inclusion of such function calls in
+	    # kernels regions (lack of explicit support of 'CMPLX' intrinsic in
+	    # PSyIR, reported) 
+	    ( lib_fortran.[fF]90 ) \
+		    sed -i -e 's/[[:space:]]CMPLX( 0.e0/ SUM( 0.e0/g' -e 's/[[:space:]]CMPLX( ztmp/ SUM( ztmp/g' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*function SIGN_/a !$acc routine seq' -e '/^\s*function DDPDD/a !$acc routine seq' | \
+		    sed -e 's/[[:space:]]SUM(0.e0/ CMPLX(0.e0/g' -e 's/SUM(ztmp/CMPLX(ztmp/g' ;;
+
+	    # Special case for trcini_medusa.F90 and trcrst_medusa.F90: adjust
+	    # non-standard FORMAT statements
+	    ( trcrst_medusa.[fF]90 | trcini_medusa.[fF]90 ) \
+		    sed -i -e "s/element:'a10/element:',a10/" ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+	    # Special cases to enable treatment of 'tiny' and 'huge' intrinsics
+	    # (pending PSyclone issue #1161) Note, includes treatment of
+	    # 'tiny(0.)'
+	    ( trcopt_medusa.[fF]90 | btd_carbchem.[fF]90 | btd_detsink.[fF]90 | btd_ironscav.[fF]90 | btd_phytopkn.[fF]90 | btd_zoopkn.[fF]90 ) \
+		    sed -i -e "/^[mM][oO][dD][uU][lL][eE]/a use loc_intrinsics, only : loc_tiny" -e "s/[tT][iI][nN][yY](0.)/loc_tiny(0)/g" -e "s/[tT][iI][nN][yY](/loc_tiny(/g" ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e "s/use loc_intrinsics, only : loc_tiny//" -e "s/loc_tiny(0)/tiny(0.)/g" -e "s/loc_tiny(/tiny(/g" ;; \
+	    ( mocsy_phsolvers.[fF]90 ) \
+		    sed -i -e "/^[mM][oO][dD][uU][lL][eE]/a use loc_intrinsics, only : loc_huge" -e "s/[hH][uU][gG][eE](/loc_huge(/" ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e "s/use loc_intrinsics, only : loc_huge//" -e "s/loc_huge(/huge(/" | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(.*)/a !$acc routine seq' | \
+		    sed -e '/^\s*[fF][uU][nN][cC][tT][iI][oO][nN].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[fF][uU][nN][cC][tT][iI][oO][nN].*(.*)/a !$acc routine seq' |
+		    sed -r -e 's/(^.*:: (niter_at[a-z]*) =.*$)/\1\n!$acc declare copyin(\2)/' ;;
+
+	    # Special case for trcoxy_medusa.F90: add "!$acc routine seq"
+	    # directive to subroutines trc_oxy_medusa, oxy_sato, and
+	    # oxy_schmidt (pending PSyclone issue #1041)
+	    ( trcoxy_medusa.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] trc_oxy_medusa(/a !$acc routine seq' \
+		        -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] oxy_sato(/a !$acc routine seq' \
+		        -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] oxy_schmidt(/a !$acc routine seq' ;;
+
+	    # Special case for trcdms_medusa.F90: add "!$acc routine seq"
+	    # directive to subroutine trc_dms_medusa and dms_flux_ocn (pending
+	    # PSyclone issue #1041)
+	    ( trcdms_medusa.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] trc_dms_medusa(/a !$acc routine seq' \
+		        -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] dms_flux_ocn(/a !$acc routine seq' ;;
+
+            ( sms_medusa.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/:: dmscut/a !$acc declare copyin(dmscut)' \
+		        -e '/:: dmsmin/a !$acc declare copyin(dmsmin)' \
+		        -e '/:: dmsslp/a !$acc declare copyin(dmsslp)' ;;
+
+	    # Special case for btd_gastrans.F90: add "!$acc routine seq"
+	    # directive to subroutines btd_gastrans (pending PSyclone issue
+	    # #1041)
+	    ( btd_gastrans.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] btd_gastrans(/a !$acc routine seq' ;;
+
+	    # Special case for btd_airsea.F90: append "private" attribute to
+	    # all "!$acc loop independent collapse(2)" directives (for
+	    # simplicity in some instances redundantly)
+	    ( btd_airsea.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e 's/collapse(2)/collapse(2) private(f_kwo2_dum,dms_surf,dms_andr,dms_simo,dms_aran,dms_hall,dms_andm,dms_flux)/' ;;
+
+	    # Special case for mocsy_{wrapper,gasflux,singledouble}.F90: add
+	    # "!$acc routine seq" directive to all procedures (pending PSyclone
+	    # issue #1041), and enable processing of subroutine 'schmidt_co2'
+	    # by removing use of 'data' statement (reported)
+	    ( mocsy_wrapper.[fF]90 | mocsy_gasflux.[fF]90 | mocsy_singledouble.[fF]90 ) \
+		    sed -i -r -e 's/^[[:space:]]*[dD][aA][tT][aA] (a.) \/([^\/]*)\/$/\1 = \2/' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(.*)/a !$acc routine seq' ;;
+
+	    # Regular PSyclone processing of the source code within src/OCE 
+	    ( divhor.[fF]90 | dyn*[fF]90 | eosbn2.F90 | ldf*[fF]90 | sbc*[fF]90 | sshwzv.[fF]90 | tra*[fF]90 | wet_dry.[fF]90 | zdf*[fF]90 | zpshde.[fF]90 | dom*.[fF]90 | closea.[fF]90 | daymod.[fF]90 | depth_e3.[fF]90 | dtatsd.[fF]90 | isc*[fF]90 | istate.[fF]90 | phycst.[fF]90 | usr*[fF]90 | step.[fF]90 | lbc*[fF]90 | mppini.[fF]90 | lib_mpp.[fF]90 | cpl_oasis3.[fF]90 | cyclone.[fF]90 | fldread.[fF]90 | geo2ocean.[fF]90 | ocealb.[fF]90 | tideini.[fF]90 | tide_mod.[fF]90 | updtide.[fF]90 | trd*[fF]90 | sto*[fF]90 | obs*[fF]90 | julian.[fF]90 | mpp_map.[fF]90 | iom*[fF]90 | in_out_manager.[fF]90 | prtctl.[fF]90 | restart.[fF]90 | icb*.[fF]90 | flo*.[fF]90 | cool_skin.[fF]90 | diurnal_bulk.[fF]90 | solfrac_mod.[fF]90 | step_diu.[fF]90 | dia*[fF]90 | crs*[fF]90 | c1d.[fF]90 | dtauvd.[fF]90 | dyncor_c1d.[fF]90 | step_c1d.[fF]90 | bdy*[fF]90 | asm*[fF]90 | lib_cray.[fF]90 | lib_fortran.[fF]90 | nemo*[fF]90 | oce.[fF]90 | par_kind.[fF]90 | par_oce.[fF]90 | step.[fF]90 | step_oce.[fF]90 | timing.[fF]90 | trc_oce.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+            # Regular PSyclone processing of source code within src/TOP used for MEDUSA
+            # setup (note, matches files src/OCE/trc_oce.F90 and
+            # src/OCE/TRD/trdtrc.F90, both of which should have been matched already)) 
+	    ( *medusa*.[fF]90 | btd*.[fF]90 | *trc*.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+	    # Do nothing with other NEMO source-code files
+	    ( * ) \
+		    cat ${fname}.preprocessed.$$ ;;
+
+        esac;
+    fi
+    # Remove temporary file
+    rm ${fname}.preprocessed.$$
+else
+    # Transparent call of wrapped pre-processor
+    $FPP ${FPP_ARG}
+fi
Index: NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/fpp_wrapper_MEDUSA.sh
===================================================================
--- NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/fpp_wrapper_MEDUSA.sh	(revision 15762)
+++ NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/fpp_wrapper_MEDUSA.sh	(revision 15762)
@@ -0,0 +1,143 @@
+#! /bin/bash
+#======================================================================
+#
+# fpp_wrapper.sh <FPP command> <wrapping mode> <FPP arguments>
+#
+# Wrapped FPP command for optional additional PSyclone processing of
+# source code
+#
+#======================================================================
+# History : 4.0 ! 2021/2022 (S. Mueller) Original version
+#======================================================================
+
+#----------------------------------------------------------------------
+# NEMO 4.0 , NEMO Consortium (2022)
+# $Id$ 
+# Software governed by the CeCILL license (see ./LICENSE)
+#----------------------------------------------------------------------
+
+# Path to PSyclone transformation script
+PSyclone_TS=../mk_psyclone/psyclone_transformation_MEDUSA.py
+
+FPP=$1; shift
+# Extended pre-processing mode
+MODE=$1; shift
+# Original arguments
+FPP_ARG=$@; shift $(($# - 1))
+#
+# PSyclone pre-processing of NEMO source-code files
+if [ ${MODE} == 'PSyclone' ]; then
+    # Filename of file to be processed
+    fname=$1
+    # Call wrapped preprocessor, store in temporary file
+    $FPP -o ${fname}.preprocessed.$$ ${FPP_ARG}
+    # PSyclone pre-processing
+    # PSyclone processing, unless ACC directives present
+    if `grep -q '!\$ACC' ${fname}`; then
+	cat ${fname}.preprocessed.$$
+    else
+        case ${fname} in
+
+	    # Do NOT process some files that do not pass the PSyclone pass-through test (see also PSyclone file examples/nemo/scripts/process_nemo.sh)
+	    ( obs_utils.[fF]90 | par_kind.[fF]90 | storng.[fF]90 | diurnal_bulk.[fF]90 | timing.[fF]90 | obs_fbm.[fF]90 | sbccpl.F90) \
+		    cat ${fname}.preprocessed.$$ ;;
+
+	    # Do NOT process source code outside of src/TOP (catch files in src/OCE that match *trc*.[fF]90 pattern used below)
+	    ( trc_oce.[fF]90 | trdtrc.[fF]90 ) \
+		    cat ${fname}.preprocessed.$$ ;;
+
+	    # Special case for mocsy_mainmod.F90: enable processing of subroutine 'constants' by removing use of 'data' statement (reported),
+	    #                                     counteract PSyclone bug #1558,
+	    #                                     add "!$acc routine seq" (pending PSyclone issue #1041)
+	    ( mocsy_mainmod.[fF]90 ) \
+		    sed -i -e '/DIMENSION(12) :: a0/d' -r -e 's/DATA ([^\/]*)\//real(kind=wp),dimension(12) :: \1 = (\//' -e 's/\/$/\/)/' -r -e 's/= \(\/[[:space:]]*12[[:space:]]*\*([^\/]*)\/\)/= \1/' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e 's/ :: onehalf/, parameter :: onehalf=0.5d0/' -e 's/ :: two/, parameter :: two=2.0d0/' -e 's/ :: three/, parameter :: three=3.0d0/' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(.*)/a !$acc routine seq' | \
+		    sed -e '/^\s*[fF][uU][nN][cC][tT][iI][oO][nN].*(.*)/a !$acc routine seq' ;;
+
+	    # Special case for lib_fortran.F90: add "!$acc routine seq" directive to SIGN function (pending PSyclone issue #1041)
+	    #                                   temporarily rename some 'CMPLX' intrinsic-function calls to 'SUM' calls to facilitate
+	    #                                   the inclusion of such function calls in kernels regions (lack of explicit support
+	    #                                   of 'CMPLX' intrinsic in PSyIR, reported) 
+	    ( lib_fortran.[fF]90 ) \
+		    sed -i -e 's/[[:space:]]CMPLX( 0.e0/ SUM( 0.e0/g' -e 's/[[:space:]]CMPLX( ztmp/ SUM( ztmp/g' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*function SIGN_/a !$acc routine seq' -e '/^\s*function DDPDD/a !$acc routine seq' | \
+		    sed -e 's/[[:space:]]SUM(0.e0/ CMPLX(0.e0/g' -e 's/SUM(ztmp/CMPLX(ztmp/g' ;;
+
+	    # Special case for trcini_medusa.F90 and trcrst_medusa.F90: adjust non-standard FORMAT statements
+	    ( trcrst_medusa.[fF]90 | trcini_medusa.[fF]90 ) \
+		    sed -i -e "s/element:'a10/element:',a10/" ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+	    # Special cases to enable treatment of 'tiny' and 'huge' intrinsics (pending PSyclone issue #1161)
+	    # Note, includes treatment of 'tiny(0.)'
+	    ( trcopt_medusa.[fF]90 | btd_carbchem.[fF]90 | btd_detsink.[fF]90 | btd_ironscav.[fF]90 | btd_phytopkn.[fF]90 | btd_zoopkn.[fF]90 ) \
+		    sed -i -e "/^[mM][oO][dD][uU][lL][eE]/a use loc_intrinsics, only : loc_tiny" -e "s/[tT][iI][nN][yY](0.)/loc_tiny(0)/g" -e "s/[tT][iI][nN][yY](/loc_tiny(/g" ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e "s/use loc_intrinsics, only : loc_tiny//" -e "s/loc_tiny(0)/tiny(0.)/g" -e "s/loc_tiny(/tiny(/g" ;; \
+	    ( mocsy_phsolvers.[fF]90 ) \
+		    sed -i -e "/^[mM][oO][dD][uU][lL][eE]/a use loc_intrinsics, only : loc_huge" -e "s/[hH][uU][gG][eE](/loc_huge(/" ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e "s/use loc_intrinsics, only : loc_huge//" -e "s/loc_huge(/huge(/" | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(.*)/a !$acc routine seq' | \
+		    sed -e '/^\s*[fF][uU][nN][cC][tT][iI][oO][nN].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[fF][uU][nN][cC][tT][iI][oO][nN].*(.*)/a !$acc routine seq' |
+		    sed -r -e 's/(^.*:: (niter_at[a-z]*) =.*$)/\1\n!$acc declare copyin(\2)/' ;;
+
+	    # Special case for trcoxy_medusa.F90: add "!$acc routine seq" directive to subroutines trc_oxy_medusa, oxy_sato, and oxy_schmidt (pending PSyclone issue #1041)
+	    ( trcoxy_medusa.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] trc_oxy_medusa(/a !$acc routine seq' \
+		        -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] oxy_sato(/a !$acc routine seq' \
+		        -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] oxy_schmidt(/a !$acc routine seq' ;;
+
+	    # Special case for trcdms_medusa.F90: add "!$acc routine seq" directive to subroutine trc_dms_medusa and dms_flux_ocn (pending PSyclone issue #1041)
+	    ( trcdms_medusa.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] trc_dms_medusa(/a !$acc routine seq' \
+		        -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] dms_flux_ocn(/a !$acc routine seq' ;;
+
+            ( sms_medusa.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/:: dmscut/a !$acc declare copyin(dmscut)' \
+		        -e '/:: dmsmin/a !$acc declare copyin(dmsmin)' \
+		        -e '/:: dmsslp/a !$acc declare copyin(dmsslp)' ;;
+
+	    # Special case for btd_gastrans.F90: add "!$acc routine seq" directive to subroutines btd_gastrans (pending PSyclone issue #1041)
+	    ( btd_gastrans.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE] btd_gastrans(/a !$acc routine seq' ;;
+
+	    # Special case for btd_airsea.F90: append "private" attribute to all "!$acc loop independent collapse(2)" directives (for simplicity in some instances redundantly)
+	    ( btd_airsea.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e 's/collapse(2)/collapse(2) private(f_kwo2_dum,dms_surf,dms_andr,dms_simo,dms_aran,dms_hall,dms_andm,dms_flux)/' ;;
+
+	    # Special case for mocsy_wrapper.F90: add "!$acc routine seq" directive to all procedures (pending PSyclone issue #1041),
+	    #                                     enable processing of subroutine 'schmidt_co2' by removing use of 'data' statement (reported)
+	    ( mocsy_wrapper.[fF]90 | mocsy_gasflux.[fF]90 | mocsy_singledouble.[fF]90 ) \
+		    sed -i -r -e 's/^[[:space:]]*[dD][aA][tT][aA] (a.) \/([^\/]*)\/$/\1 = \2/' ${fname}.preprocessed.$$ ; \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(/N;/.*\n&/a !$acc routine seq' | \
+		    sed -e '/^\s*[sS][uU][bB][rR][oO][uU][tT][iI][nN][eE].*(.*)/a !$acc routine seq' ;;
+
+	    # Regular PSyclone processing of source code within src/TOP (note, src/OCE/trc_oce.F90 andn src/OCE/TRD/trdtrc.F90 caught above)
+	    ( *medusa*.[fF]90 | btd*.[fF]90 | *trc*.[fF]90 | traadv_mus.[fF]90 ) \
+		    psyclone -api nemo -l output -s ${COMPIL_DIR}/${PSyclone_TS} -oalg /dev/null ${fname}.preprocessed.$$ | tail -n +2 ;;
+
+	    # Do nothing with other NEMO source-code files
+	    ( * ) \
+		    cat ${fname}.preprocessed.$$ ;;
+
+        esac;
+    fi
+    # Remove temporary file
+    rm ${fname}.preprocessed.$$
+else
+    # Transparent call of wrapped pre-processor
+    $FPP ${FPP_ARG}
+fi
Index: NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/psyclone_transformation_BENCHO2M.py
===================================================================
--- NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/psyclone_transformation_BENCHO2M.py	(revision 15762)
+++ NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/psyclone_transformation_BENCHO2M.py	(revision 15762)
@@ -0,0 +1,1139 @@
+#======================================================================
+#
+# PSyclone transformation script for the PSyclone processing of the
+# BENCHO2M NEMO-MEDUSA model setup (small NEMO-MEDUSA test case)
+#
+#======================================================================
+# History : 4.0 ! 2022 (S. Mueller) Original version
+#======================================================================
+
+#----------------------------------------------------------------------
+# NEMO 4.0 , NEMO Consortium (2022)
+# $Id$
+# Software governed by the CeCILL license (see ./LICENSE)
+#----------------------------------------------------------------------
+
+from psyclone.transformations         import ACCKernelsTrans, ACCLoopTrans
+from psyclone.psyir.transformations   import ArrayRange2LoopTrans
+from psyclone.psyir.nodes             import CodeBlock, IfBlock, Assignment, Call, Reference, Literal
+from psyclone.psyir.nodes             import BinaryOperation, UnaryOperation, ArrayReference, Range
+from psyclone.psyir.nodes             import ArrayOfStructuresReference
+from psyclone.psyir.nodes             import Loop, Node, ACCKernelsDirective, ACCLoopDirective, Return
+from psyclone.psyir.nodes             import StructureReference, ArrayMember
+from psyclone.psyir.symbols           import DataSymbol, ContainerSymbol, ImportInterface, RoutineSymbol
+from psyclone.psyir.symbols           import ScalarType, UnknownFortranType, ArrayType
+from psyclone.psyir.symbols           import INTEGER_TYPE, CHARACTER_TYPE, BOOLEAN_TYPE, DataTypeSymbol
+from psyclone.psyir.frontend.fparser2 import Fparser2Reader
+from psyclone.nemo                    import NemoLoop, NemoKern
+from fparser.two                      import Fortran2003
+from fparser.two.utils                import walk
+from re                               import match
+
+# Exclude a range from modules from acceleration (regular expressions)
+ACC_MODULES_EXCLUDE_RE = [ '^mocsy_.*' ]
+
+# Exclude a range of invokes from accelerated code as the resulting ACC
+# kernels fail to be compiled for some compiler versions
+ACC_INVOKES_EXCLUDE = [ 'trc_bc',                    # trcbc.F90
+                        'trc_dmp_clo',               # trcdmp.F90
+                        'trc_nam_ice',               # trcice.F90
+                        'p4z_ini',                   # trcini_pisces.F90
+                        'p2z_ini',                   # trcini_pisces.F90
+                        'prt_ctl_trc',               # prtctl_trc.F90, avoid assignment of character array withing kernels region
+                        'zdf_osm',                   # zdfosm.F90
+                        'obs_wmo_init',              # obs_types.F90, string array assignments inside loop
+                        'obs_surf_alloc',            # obs_surf_def.F90, string array assignments inside loop
+                        'obs_dist2corners',          # obs_averg_h2d.F90, function call ('grt_cir_dist') confused with array access
+                        'obs_prof_alloc',            # obs_profile_def.F90, string array assignments inside loop
+                        'obs_rea_prof',              # obs_read_prof.F90
+                        'obs_rea_surf',              # obs_read_surf.F90, collased-loop bounds stored in derived-type structure
+                        'obs_wri_prof',              # obs_write.F90, string array assignments inside loop
+                        'obs_wri_surf',              # obs_write.F90, string array assignments inside loop
+                        'maxdist',                   # obs_grid.F90, function call ('grt_cir_dist') confused with array access
+                        'flo_add_new_ariane_floats', # flodom.F90
+                        'icb_utl_copy',              # icbutl.F90, function call ('ice_var_sshdyn') confused with array access
+                        'icb_rst_write',             # icbrst.F90, function call ('nf90_put_var') confused with array access
+                        'icb_dia_init',              # icbdia.F90, function call ('icb_utl_mass') confused with array access
+                        'flo_rst',                   # florst.F90, compilation fails
+                        'icb_clv_flx',               # icbclv.F90, compilation fails
+                        'icb_init',                  # icbini.F90, compilation fails
+                        'icb_thm',                   # icbthm.F90, compilation fails
+                        'icb_stp',                   # icbstp.F90, compilation fails
+                        'crs_dom_ope_2d',            # crsdom.F90, compilation fails
+                        'crs_dom_def',               # crsdom.F90, compilation fails
+                        'dia_dct_wri',               # diadct.F90, string assignment in kernels region
+                        'bdy_dta',                   # bdydta.F90, compilation fails
+                        'bdy_vol2d'                  # bdyvol.F90, loop bounds stored in derived-type structure
+                      ]
+
+# Exclude code that is conditional to a single logical (e.g., debug
+# output, tracer-trend diagnostics) from accelerated code
+ACC_CONDITIONALS_REMOVE = [ 'ln_msa_debug',   # MEDUSA debug output
+                            'l_trdtrc'        # Tracer-trend diagnostic
+                          ]
+
+# Permit some subroutine calls inside accelerated kernels regions
+ACC_SUBROUTINES_PERMIT_CALLS = [ 'trc_oxy_medusa',
+                                 'oxy_sato',
+                                 'oxy_schmidt',
+                                 'trc_dms_medusa',
+                                 'dms_flux_ocn',
+                                 'btd_gastrans',
+                                 'mocsy_interface',
+                                 'ddpdd'
+                               ]
+
+# Exclude vertical loops from accelerated kernels regions in named invokes
+# Note: invoke 'tra_ldf_iso' will be added to this list if AC15_VARIANT is
+# selected
+ACC_INVOKES_EXCLUDE_VERTICAL_LOOPS = [ 'tra_zdf_imp'
+                                     ]
+
+# Prevent inclusion of if constructs inside kernels regions
+ACC_INVOKES_EXCLUDE_IFBLOCKS = { 'lbc_nfd_3d_ptr', # facilitates parallelism detection by compiler
+                                 'lbc_nfd_4d_ptr'  # facilitates parallelism detection by compiler
+                                 }
+
+# Select variants for some transformations or accelerations
+AC15_VARIANT = False
+
+# Select refined timing of MEDUSA component
+TIMING_MEDUSA_REFINED = True
+
+def kernels_region_grow(kernel, direction=0, accept_loop_levels=True, exclude_if=False):
+    # grow a kernels region by identifying whether a parent or neighbouring
+    # node can be included in a kernels region
+    #
+    # kernel: base node for futher growth of the kernels region
+    # direction:          indicates the direction of growth, dir=0 grows the kernels
+    #                     region upwards in the PSyIR, dir in [-1,1] grow the
+    #                     kernels region in either lateral direction,
+    #                     respectively
+    # accept_loop_levels: accept NemoLoop nodes of loop type 'levels'
+    # exclude_if:         if True, prevent inclusion of if constructs inside
+    #                     kernels regions
+    #
+    # Identify parent or neighbouring node
+    pos = kernel.position
+    if abs(direction) == 1 and pos+direction >= 0 and pos+direction < len(kernel.parent.children):
+        kernel = kernel.parent.children[pos+direction]
+    elif direction == 0:
+        kernel = kernel.parent.parent
+    else:
+        kernel = None
+    # To accept, the new node must be either a NemoLoop, IfBlock or Assignment
+    # object, but it must not include CodeBlock or Call objects, not include a
+    # kernels region, not include a Return statement, and not already be inside
+    # a kernels region
+    if kernel:
+        if not isinstance(kernel, (NemoLoop,IfBlock,Assignment)) or kernel.walk((CodeBlock,ACCKernelsDirective,Return)) or has_excluded_call(kernel) or kernel.ancestor(ACCKernelsDirective):
+            kernel = None
+        elif exclude_if and kernel.walk(IfBlock):
+            kernel = None
+        elif not accept_loop_levels and isinstance(kernel, NemoLoop):
+            for l in kernel.walk(NemoLoop):
+                if l.loop_type == 'levels':
+                    kernel = None
+                    break
+    return kernel
+
+# Test for presence of unallowed call in kernels region 
+def has_excluded_call(kernel):
+    for call in kernel.walk(Call):
+        if not call.routine.name in ACC_SUBROUTINES_PERMIT_CALLS:
+            return True
+    return False
+
+# Create allocation statement for variable represented by symbol 'sym' of the
+# shape supplied in argument 'shape' (if list of stings) or of the
+# shape of the symbol named in argument 'shape' (if string); depending on
+# 'cond', make allocation dependent on the allocation status of the variable.
+def node_allocate(sym, shape, cond=True):
+    if sym.is_array:
+        rank = len(sym.shape)
+    elif isinstance(sym.datatype, UnknownFortranType):
+        rnge = match('^.*\(([ :,]*)\).+$',sym.datatype.declaration)
+        if rnge:
+            rank = rnge.group(1).count(':')
+    else:
+        return None
+    if isinstance(shape, str):
+        shape_ref = shape
+        shape = []
+        for n in range(rank):
+            shape += ['LBOUND('+shape_ref+','+str(n+1)+'):UBOUND('+shape_ref+','+str(n+1)+')']
+    allocate = CodeBlock([Fortran2003.Allocate_Stmt('allocate('+sym.name+'('+','.join(shape)+'))')], CodeBlock.Structure.STATEMENT)
+    if cond:
+        allocatd = CodeBlock([Fortran2003.Intrinsic_Function_Reference('allocated('+sym.name+')')], CodeBlock.Structure.EXPRESSION)
+        allocate = IfBlock.create(UnaryOperation.create(UnaryOperation.Operator.NOT, allocatd), [allocate])
+    return allocate
+
+# Compute simple case-insensitive hash value of up to the five first characters
+# of a string
+#
+# Note, the five character limit is applied in order to be able to represent
+# the hash value as an integer of standard precision (31 bit)
+def opt_hashval(opt):
+    val = 0
+    for nc in range(min(5,len(opt.upper()))):
+        val += (ord(opt[nc])&31)*32**nc
+    return val
+
+#########################################################################
+# Main transformation procedure (called by PSyclone's 'psyclone' command)
+#########################################################################
+
+def trans(psy):
+
+    # Module name
+    module_name = psy.name[4:-4]
+
+##################################################
+# Section 1 - Specific source-code transformations
+##################################################
+
+    # Loop through all invokes
+    for invoke in psy.invokes.invoke_list:
+        sched = invoke.schedule
+        st = sched.symbol_table
+        pst = st.parent_symbol_table()
+
+        # Refined MEDUSA timing
+        #
+        # Additional timing instructions that refine 'trc_bio_medusa' section
+        if TIMING_MEDUSA_REFINED and invoke.name == 'trc_bio_medusa':
+            labels = ['trc_bio_medusa_ini', 'trc_bio_msa_nan', 'trc_bio_msa_atm', \
+                      'trc_bio_msa_cch',    'trc_bio_msa_vlp', 'trc_bio_msa_trc', \
+                      'trc_bio_msa_mip',    'trc_bio_msa_chk', 'trc_bio_msa_as',  \
+                      'trc_bio_msa_pln',    'trc_bio_msa_isc', 'trc_bio_msa_det', \
+                      'trc_bio_msa_upd',    'trc_bio_msa_dia', 'trc_bio_msa_ben'  \
+                      ]
+            for n in range(len(labels)):
+                # Find parent nodes (where they change) and positions (incl.
+                # anticipation of position change following addition of start
+                # statement) for the addition of the various timing statements
+                if n==0:
+                    snode = [ib for ib in sched.walk(IfBlock) if isinstance(ib.condition, Reference) and ib.condition.name=='ln_timing'][0]
+                    sp    = snode.parent
+                    spos  = snode.position+1
+                    enode = [c for c in sched.walk(Call) if c.routine.name=='btd_ini'][0]
+                    ep    = enode.parent
+                    epos = enode.position+2
+                elif n==1:
+                    spos = epos+1
+                    epos = [ib for ib in sched.walk(IfBlock) if isinstance(ib.condition, Reference) and ib.condition.name=='lk_oasis'][0].position+1
+                elif n==2:
+                    spos = epos+1
+                    for ib in sched.walk(IfBlock):
+                        if [r for r in ib.walk(Reference) if r.name == 'ln_rsttr']:
+                            epos = ib.position+2
+                            break
+                elif n in [3,4]:
+                    spos = epos+2
+                    epos = spos+2
+                elif n==5:
+                    sp = sp.children[spos+1].loop_body
+                    spos = 0
+                    ep = sp
+                    epos = [ib for ib in sp.walk(IfBlock) if isinstance(ib.condition, Reference) and ib.condition.name == 'ln_msa_omipdic'][0].position+1
+                elif n==6:
+                    spos = epos+1
+                    epos = [c for c in sp.walk(Call) if c.routine.name == 'trc_bio_check'][0].position+2
+                elif n==7:
+                    spos = epos+1
+                    epos = [c for c in sp.walk(Call) if c.routine.name == 'btd_air_sea'][0].ancestor(IfBlock).position+1
+                elif n<13:
+                    spos = epos+1
+                    epos = spos+2
+                elif n==13:
+                    spos = epos+1
+                    epos = spos+4
+                elif n==14:
+                    sp   = vlte.parent
+                    spos = vlte.position+1
+                    ep   = sp
+                    epos = len(ep.children)
+                # Add start statement
+                ts = IfBlock.create(Reference(st.lookup('ln_timing')), [Call.create(RoutineSymbol('timing_start'), [Literal(labels[n], CHARACTER_TYPE)])])
+                sp.addchild(ts, index=spos)
+                # Add stop statement
+                te = IfBlock.create(Reference(st.lookup('ln_timing')), [Call.create(RoutineSymbol('timing_stop'), [Literal(labels[n], CHARACTER_TYPE)])])
+                ep.addchild(te, index=epos)
+                # Remeber where the main vertical loop ends
+                if n == 4:
+                    vlte = te
+
+        # ######## AC22 ####### (part 2)
+        #
+        # Some individual calls of Fortran intrinsic functions in invoke
+        # 'stp_ctl' can be accelerated when run on the GPU (C. Dearden, pers.
+        # comm.). While these reductions are recognised as NemoKern objects,
+        # they are not put inside kernels regions as they are represented by a
+        # scalar assignment of the result of an individual individual function
+        # call; to trigger their inclusion into a kernels region these
+        # assignments are disguised as array assignments.
+        if invoke.name == 'stp_ctl':
+            for a in sched.walk(Assignment):
+                if isinstance(a.lhs, ArrayReference) and a.lhs.name == 'zmax' and not a.lhs.is_full_range(0):
+                    a.lhs.children[0].replace_with(Range.create(a.lhs.children[0].copy(), a.lhs.children[0].copy()))
+
+        # ######## AC22 ####### (part 3)
+        #
+        # Specific treatment of invoke 'ldf_eiv' that includes calls of the
+        # Fortran instrinsic function 'sin' inside kernels region, which leads
+        # to a difference between results when this region is executed on the
+        # CPU or GPU, repsectively. As the 'sin' function is part of the
+        # comutation of a temporary scalar value inside the loop, that does not
+        # depend on time-varying variables, all of those temoprary values can
+        # be precomputed on the CPU and stored during the initial invokation of
+        # 'ldf_eiv'.
+        if invoke.name == 'ldf_eiv':
+            # Collect nodes that represent the assignment containing the 'sin'
+            # call ('a'), the array argument variable of the call ('v'), and
+            # the encompassing 2-dimension loop nest ('ll')
+            for uo in sched.walk(UnaryOperation):
+                if uo.operator == UnaryOperation.Operator.SIN and uo.walk(ArrayReference) and uo.walk(ArrayReference)[0].name == 'gphit':
+                    a = uo.ancestor(Assignment)
+                    v = uo.walk(ArrayReference)[0]
+                    ll = uo.ancestor(Loop).ancestor(Loop)
+                    break
+            # Create new array module symbol based on the name of the original
+            # temporary variable
+            aname = a.lhs.name
+            sym = pst.new_symbol(aname+'_'+invoke.name, symbol_type=DataSymbol, datatype=ArrayType(a.lhs.symbol.datatype, shape=[ArrayType.Extent.DEFERRED]*len(v.indices)))
+            # Create the corresponding allocation statement
+            alloc = node_allocate(sym, v.name)
+            sched.addchild(alloc, index=0)
+            # Upgrade the assignment to an array-element assignment
+            for r in [ref for ref in ll.walk(Reference) if ref.name == aname]:
+                r.replace_with(ArrayReference.create(sym,v.copy().pop_all_children()))
+            # Create a copy of the 2-dimensional loop nest around the
+            # assignment, delete the whole loop body, and move the assignment
+            # from the original location to the loop copy
+            ll = ll.copy()
+            ll.walk(Loop)[1].loop_body.pop_all_children()
+            ll.walk(Loop)[1].loop_body.addchild(a.detach())
+            # Insert modified loop-nest copy as a sibling of the allocation
+            # statement for the newly created look-up array
+            alloc.walk(IfBlock)[0].if_body.addchild(ll,index=1)
+
+        # ######## AC21 #######
+        #
+        # Special treatment of 'lbc_nfd_3d_ptr' and 'lbc_ndf_4d_ptr' to improve
+        # acceleration of two invokes of the north-fold exchange
+        # implementation: transform array assignments into explicit loops and
+        # move a loop-index computation inside the newly created loop nest to
+        # enable loop collapse (C. Dearden, pers. comm.)
+        if invoke.name in [ 'lbc_nfd_3d_ptr', 'lbc_nfd_4d_ptr' ]:
+            # Loop indices and upper loop boundaries of the newly created
+            # loops
+            lindex = [st.lookup('jk'),  st.lookup('jl')]
+            lupper = [st.lookup('ipk'), st.lookup('ipl')]
+            # Find array assigments (incl. array assignments of structure
+            # members) and store them in a list (the order matters for the neighbour
+            # detection below), and lookup dictionaries that indicates the number transformed dimensions
+            # position of the assignment amongst its siblings (the *original*
+            # positions are required for the neighbour detection below)
+            aa = []
+            aan = {}
+            aap = {}
+            for a in [node for node in sched.walk(Assignment) if isinstance(node.lhs, StructureReference) and  isinstance(node.lhs.member, ArrayMember)]:
+                # Find first of a contiguous set of implicit loops
+                # amongst the rightmost array dimensions
+                ndim = len(a.lhs.member.indices)
+                noff = 0
+                for n in range(ndim, 0, -1):
+                    if a.lhs.member.is_full_range(n-1):
+                        noff += 1
+                    else:
+                        break
+                if (noff > 0):
+                    aa.append(a)
+                    aan[a] = noff
+                    aap[a] = a.position
+            # Replace array ranges with indices for the new loops
+            for a in aa:
+                noff = aan[a]
+                for n in range(noff):
+                    # Transform implicit loop to explicit loop
+                    a.lhs.member.indices[ndim-n-1].replace_with(Reference(lindex[noff-n-1]))
+                    for r in a.rhs.walk((ArrayReference,StructureReference)):
+                        if isinstance(r,StructureReference) and isinstance(r.member, ArrayMember):
+                            r = r.member
+                        else:
+                            continue
+                        ndimr = len(r.indices)
+                        if r.is_full_range(ndimr-n-1):
+                            r.indices[ndim-n-1].replace_with(Reference(lindex[noff-n-1]))
+            # Group clusters of adjacent array assignments and create the new
+            # loops
+            loop_body = []
+            for naa in range(len(aa)):
+                a = aa[naa]
+                noff = aan[a]
+                p = a.parent
+                apos = aap[a]
+                if naa < len(aa)-1:
+                    pnext = aa[naa+1].parent
+                    noffnext = aan[aa[naa+1]]
+                    aposnext = aap[aa[naa+1]]
+                else:
+                    pnext = None
+                    noffnext = -1
+                    aposnext = -1
+                acopy = a.copy()
+                loop_body += [acopy]
+                # Additionally move r.h.s. array-index computation (if
+                # any) inside the newly created loop to facilitate loop
+                # collapse of encompassing loops (if possible)
+                if a.position > 0 and isinstance(p.children[a.position-1], Assignment) and p.children[a.position-1].lhs.name in ['ijt', 'iju']:
+                    loop_body = [p.children[a.position-1].detach()] + loop_body
+                if not ( apos == aposnext-1 and p == pnext and noff == noffnext ):
+                    # Create explicit loop(s) and pack the loop body inside a
+                    # NemoKernel to trigger kernels region detection
+                    for n in range(noff):
+                        l = NemoLoop.create(lindex[noff-n-1], Literal('1', INTEGER_TYPE), Reference(lupper[noff-n-1]), Literal('1', INTEGER_TYPE), [NemoKern(loop_body)])
+                        a.replace_with(l)
+                        a = l.loop_body.children[0]
+                        loop_body = a.children[0].copy().pop_all_children()
+                    loop_body = []
+                else:
+                    a.detach()
+
+        # ##### AC17/AC18 #####
+        #
+        # Special treatment of 'dom_vvl_interpol' subroutine and lbclnk and
+        # lbcnfd modules: replace comparisons with option strings by
+        # comparisons with integer variables
+        if invoke.name == 'dom_vvl_interpol' or module_name.startswith('lbc'):
+            if invoke.name == 'dom_vvl_interpol':
+                opt_name = 'pout'
+            else:
+                opt_name = 'cd_nat'
+            # Test if symbol of that name exists and if it can be a(n array of)
+            # character string(s) (i.e. of datatype UnknownFortranType)
+            if opt_name in st.symbols_dict.keys() and isinstance(st.lookup(opt_name).datatype, UnknownFortranType):
+                # Work out if this is an array (assume rank is 1)
+                is_array = False
+                if st.lookup(opt_name).datatype.declaration.find('(:') >= 0:
+                    is_array = True
+                dtype = ScalarType(ScalarType.Intrinsic.INTEGER, ScalarType.Precision.DOUBLE)
+                if is_array:
+                    dtype = ArrayType(dtype, shape=[ArrayType.Extent.DEFERRED])
+                sym = st.new_symbol(opt_name+'_val', symbol_type=DataSymbol, datatype=dtype)
+                # Find and store all literal values referred in binary
+                # comparisons for equality and replace then with hash values
+                opts = {}
+                for ib in sched.walk(IfBlock):
+                    cnode = ib.condition
+                    for bop in ib.condition.walk(BinaryOperation):
+                        if isinstance(bop.children[0], Reference) and bop.children[0].name == opt_name and isinstance(bop.children[1], Literal):
+                            val = opt_hashval(bop.children[1].value)
+                            opts[val] = bop.children[1].copy()
+                            bop.children[0].symbol = sym
+                            bop.children[1].replace_with(Literal(str(val),INTEGER_TYPE))
+                # Create (array of) translation(s) of the dummy-argument
+                # value(s) used for the comparison
+                ibs = []
+                if is_array:
+                    sym_idx = st.new_symbol(opt_name+'_idx', symbol_type=DataSymbol, datatype=INTEGER_TYPE)
+                for val in opts.keys():
+                    if is_array:
+                        # Note: cannot use ArrayReference.create with symbols
+                        # of datatype UnknownFortranType
+                        oref = ArrayReference(st.lookup(opt_name))
+                        oref.addchild(Reference(sym_idx))
+                        nref = ArrayReference.create(sym, [Reference(sym_idx)])
+                    else:
+                        oref = Reference(st.lookup(opt_name))
+                        nref = Reference(sym)
+                    if is_array:
+                        start = BinaryOperation.create(BinaryOperation.Operator.LBOUND, Reference(st.lookup(opt_name)), Literal('1', INTEGER_TYPE))
+                        stop  = BinaryOperation.create(BinaryOperation.Operator.UBOUND, Reference(st.lookup(opt_name)), Literal('1', INTEGER_TYPE))
+                        step  = Literal('1', INTEGER_TYPE)
+                    ibs += [IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.EQ, oref, opts[val]), [Assignment.create(nref, Literal(str(val), INTEGER_TYPE))])]
+                if is_array:
+                    sched.addchild(Loop.create(sym_idx, start, stop, step, ibs), index=0)
+                    sched.addchild(node_allocate(sym, opt_name), index=0)
+                else:
+                    for ib in ibs:
+                        sched.addchild(ib, index=0)
+
+        # ######## AC15 #######
+        #
+        # Special treatment of 'tra_ldf_iso' subroutine: improve the efficiency of
+        # an accelerated IfBlock containing array assignments inside a vertical
+        # loop (previously "NMOD9"); alternatively, a substantial, although
+        # slightly lower, efficiency improvement can be achieved by preventing
+        # vertical loops from being included in kernels regions in invoke
+        # 'tra_ldf_iso'
+        acc_invokes_exclude_vertical_loops_extra = []
+        if invoke.name == 'tra_ldf_iso' and AC15_VARIANT:
+            acc_invokes_exclude_vertical_loops_extra = ['tra_ldf_iso']
+        elif invoke.name == 'tra_ldf_iso':
+            # Prepare translation of implicit to explicit loop
+            trans = ArrayRange2LoopTrans()
+            # Look for three array assignments of which the second and third
+            # are the alternatives of an if block
+            for ib in sched.walk(IfBlock):
+                pos = ib.position - 1
+                if isinstance(ib.if_body[0], Assignment) and ib.if_body[0].is_array_range and isinstance(ib.else_body[0], Assignment) and ib.else_body[0].is_array_range and ib.position > 0 and isinstance(ib.parent.children[pos], Assignment) and ib.parent.children[pos].is_array_range:
+                    # Translate array assignment to explicit loops and extract
+                    # loop meta data
+                    a1 = ib.parent.children[pos]
+                    trans.apply(a1)
+                    trans.apply(a1)
+                    sa = [ib.parent.children[pos].walk(Loop)[1].variable, ib.parent.children[pos].walk(Loop)[0].variable]
+                    # Move if block inside newly created explicit loops
+                    a1.parent.addchild(ib.detach(),index=a1.position+1)
+                    # Replace array ranges by the newly created loop variables
+                    for a in ib.if_body[0].walk(ArrayReference) + ib.else_body[0].walk(ArrayReference):
+                        n = 0
+                        for ar in a.children:
+                            if isinstance(ar, Range):
+                                ar.replace_with(Reference(sa[n]))
+                                n += 1
+
+        # ######## AC12 #######
+        #
+        # Reduce number of times the a halo exchange (procedure 'lbc_lnk' call)
+        # that are made during the MEDUSA time step: only call 'lbc_lnk' at the
+        # deepest depth layer
+        if invoke.name == 'btd_phyto_pkn':
+            for c in [call for call in sched.walk(Call) if call.routine.name == 'lbc_lnk']:
+                if 'jk' in st.symbols_dict.keys() and ('jpj' in st.symbols_dict.keys() or 'jpk' in st.symbols_dict.keys()):
+                    if not 'jpk' in st.symbols_dict.keys():
+                        st.new_symbol('jpk', symbol_type=DataSymbol, datatype=INTEGER_TYPE).copy_properties(st.lookup('jpj'))
+                    sym_jpk = st.lookup('jpk')
+                    pos = c.position
+                    c.parent.addchild(IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.EQ, Reference(st.lookup('jk')), Reference(st.lookup('jpk'))), [c.detach()]), index=pos)
+
+        # ######## AC2 ########
+        #
+        # Selective removal of conditional code segments
+        # Note, only removes segements that are conditional on
+        # single-statement logical
+        [ifb.detach() for ifb in invoke.schedule.walk(IfBlock) if isinstance(ifb.condition, Reference) and (ifb.condition.name in ACC_CONDITIONALS_REMOVE)]
+
+        # ######## AC8 ########
+        # Special treatment of the glob_sum function call in subroutine 'trc_stp'
+        # of module 'trcstp'
+        #
+        # Extract the multiplication of 3-dimensional arrays from the argument of a
+        # function call, so that it can be automatically placed inside a kernels
+        # region
+        if module_name == 'trcstp' and invoke.name == 'trc_stp':
+            # Create temporary array
+            sym = pst.new_symbol('zfn_'+sched.invoke.name, symbol_type=DataSymbol, \
+                    datatype=ArrayType(ScalarType(ScalarType.Intrinsic.REAL, ScalarType.Precision.DOUBLE), shape=[ArrayType.Extent.DEFERRED]*3))
+            sched.addchild(node_allocate(sym, 'cvol'), index=0)
+            # Separate arithmetic operations within the argument list of function calls
+            # from the function calls to enable acceleration of these operations
+            #
+            # Assume that the expression consists of an assignment of the result of the
+            # function call to a variable; the function call is represented by a
+            # CodeBlock
+            #
+            # Look for function call (i.e., Fortran2003.Structure_Constructor instance,
+            # with Fortran2003.Type_Name child) with arithmetic-expression argument
+            # (Fortran2003.Add_Operand instance present)
+            asgn_node = None
+            for asgn_node in sched.walk(Assignment):
+                for cb in asgn_node.walk(CodeBlock):
+                    for fn_ast in walk(cb.ast, Fortran2003.Structure_Constructor):
+                        if walk(fn_ast.children[0], Fortran2003.Type_Name) and walk(fn_ast.children[1], Fortran2003.Add_Operand):
+                            # AST of function-call argument
+                            arg_ast = walk(fn_ast.children[1], Fortran2003.Add_Operand)
+                            # Create assignment of arithmetic expression found in
+                            # argument, and make sure full array ranges are
+                            # explicitely specified so that the expression can be
+                            # detected as the seed for a kernels region
+                            rnge = []
+                            for n in range(len(sym.shape)):
+                                rnge += [Range.create(BinaryOperation.create(BinaryOperation.Operator.LBOUND, Reference(sym), Literal(str(n+1), INTEGER_TYPE)), \
+                                                      BinaryOperation.create(BinaryOperation.Operator.UBOUND, Reference(sym), Literal(str(n+1), INTEGER_TYPE)))]
+                            asgn_new = Assignment.create(ArrayReference.create(sym, rnge), CodeBlock(arg_ast, CodeBlock.Structure.EXPRESSION))
+                            # Include assignment node in PSyIR
+                            asgn_node.parent.addchild(asgn_new, index=asgn_node.position)
+                            # Translate AST for arithmetic expression into PSyIR sub-tree
+                            asgn_new.children[1] = Fparser2Reader()._binary_op_handler(asgn_new.rhs.ast, asgn_new)
+                            # Replace original argument by name of temporary variable
+                            cb.replace_with(CodeBlock([fn_ast.tostr().replace(arg_ast[0].tostr(), sym.name)], CodeBlock.Structure.EXPRESSION))
+
+        # ######## AC8 ########
+        #
+        # Special treatment of glob_sum_* functions (module
+        # lib_fortran)
+        #
+        # Expand the n-dimensional reduction by reducing one dimension
+        # only in each of the expanded loop nests to expose some
+        # parallelism at the expense of moderate memory requirements
+        #
+        # Note, this section makes a range of assumptions:
+        #  1) The reduction occurs in a nested NemoLoop structure
+        #     where one or more of the outermost loops may be
+        #     irrelevant for the reduction (loop count of 1, relict
+        #     from macro-based code generation)
+        #  2) The reduction operation consists of a) an assignment
+        #     that represents the current summand and b) a call to a
+        #     reduction subroutine with two arguments, of which the
+        #     second is the temporary reduction result
+        #  3) the only dummy argument that is an array represents
+        #     the array to be summed
+        #  4) the horizontal shape of the array to be summed is
+        #     identical in each call of the modified subroutine
+        #     and therefore the newly added module array does only
+        #     have to be allocated in the first call of the subroutine
+        if invoke.name.startswith('glob_sum_'):
+            # The name of the input array
+            sym = [s for s in st.symbols if s.is_array and s.is_argument][0]
+            in_var = sym.name
+            ndim = len(sym.shape)
+            # The orginal loops of the reduction that are relevant
+            # for the reduction (outermost will by first element in
+            # array)
+            loops = sched.walk(NemoLoop)[-ndim:]
+            nloops = len(loops)
+            loop_indices = [loops[n].variable for n in range(len(loops))]
+            # Only transform loop nests, not single loops
+            if len(loops) > 1:
+                # The name of the reduction-operator subroutine
+                op_proc = loops[0].walk(Call)[0].routine.name
+                # The name of the accumulation variable, and the
+                # associated symbol
+                out_var = loops[0].walk(Call)[0].children[1].name
+                sym_ov = st.lookup(out_var)
+                # Create expansion of accumulation variable as an
+                # allocatable module variable
+                if isinstance(sym_ov.datatype, UnknownFortranType):
+                    pst.new_symbol(out_var+'_'+invoke.name, symbol_type=DataSymbol, datatype=UnknownFortranType(sym_ov.datatype.declaration.replace('::',',dimension('+','.join([':']*(nloops-1))+'),allocatable ::').replace(out_var, out_var+'_'+invoke.name)))
+                    sym_ov.datatype.declaration = sym_ov.datatype.declaration.replace(out_var, 'UNUSED___'+out_var)
+                else:
+                    pst.new_symbol(out_var+'_'+invoke.name, symbol_type=DataSymbol, datatype=ArrayType(sym_ov.datatype, shape=[ArrayType.Extent.DEFERRED]*(len(loops)-1)))
+                st.rename_symbol(sym_ov, 'UNUSED___'+out_var)
+                sym_ov = pst.lookup(out_var+'_'+invoke.name)
+                # Insert allocation statement
+                sched.children.insert(0, node_allocate(sym_ov, in_var))
+                # Create new references with indices '(:,:)', '(ji,jj), '(ji,1)',
+                # '(1,1)' (3-d case) or '(:)', '(ji)', '(1)' (2-d case)
+                refs = []
+                for n in range(-1,nloops):
+                    ref_new = ArrayReference(sym_ov)
+                    for nn in range(nloops-1):
+                        if n < 0:
+                            ref_new.addchild(Range.create(BinaryOperation.create(BinaryOperation.Operator.LBOUND, Reference(sym_ov), Literal(str(nn+1), INTEGER_TYPE)),BinaryOperation.create(BinaryOperation.Operator.UBOUND, Reference(sym_ov), Literal(str(nn+1), INTEGER_TYPE))))
+                        elif ((nloops-2)-nn)-n < 0:
+                            ref_new.addchild(Literal('1', INTEGER_TYPE))
+                        else:
+                            ref_new.addchild(Reference(loop_indices[nloops-nn-1]))
+                    refs.append(ref_new)
+                # Replace existing references to the accumulation variable
+                for a in sched.walk(Reference):
+                    if a.name == 'UNUSED___'+out_var:
+                        if isinstance(a.parent, Assignment):
+                            a.replace_with(refs[0].copy())
+                        elif isinstance(a.parent, Call) and a.parent.routine.name == op_proc:
+                            a.replace_with(refs[1].copy())
+                        else:
+                            a.replace_with(refs[-1].copy())
+                        first = False
+                for n in range(nloops-1,0,-1):
+                    l = loops[n].copy()
+                    if len(l.loops()) < 3:
+                        l.loop_type = 'unknown'
+                    l.walk(Assignment)[0].detach()
+                    l.walk(Call)[0].children[0].replace_with(refs[n].copy())
+                    l.walk(Call)[0].children[1].replace_with(refs[n+1].copy())
+                    l.start_expr = Literal('2', INTEGER_TYPE)
+                    loops[0].parent.addchild(l,index=loops[0].position+1)
+    
+        # #### AC7/AC9/AC11 ###
+        #
+        # Removal of warning output from the horizontal loops in the air-sea
+        # gas-exchange and carbonate-chemistry update subroutines (invokes
+        # 'btd_air_sea' and 'btd_carbchem')
+        #
+        # With the carbonate-chemistry code resident on the GPU (i.e., calls to
+        # mocsy_interafce accepted in kernels regions), the air-sea gas-exchange
+        # and carbonate-chemistry updates can be fully parallelised with the
+        # exception of the detection of excessive tracer conditions and the 
+        # associated ermegency stop.
+        #
+        # Here, we replace the warning output and the emergency stop by a detection
+        # mechanism but delay the warning output or model halt until the nested
+        # loop over all water columns is complete; for simplicity, the output of
+	# location and and tracer values where excessive values are encountered
+	# is suppressed
+        #
+        # Identify all nodes that encompasse the test, debug output,
+        # and, if present, the halt statement
+        test_nodes_warning = []
+        if invoke.name in ['btd_air_sea', 'btd_carbchem', 'trc_bio_medusa']:
+            for l in [loop for loop in sched.walk(NemoLoop) if loop.loop_type == 'lat']:
+                for n in l.walk((Call, IfBlock)):
+                    if isinstance(n, Call) and n.routine.name == 'ctl_stop':
+                        if not n.parent.parent in test_nodes_warning:
+                            test_nodes_warning += [n.parent.parent]
+                    elif isinstance(n, IfBlock) and isinstance(n.condition, Reference) and n.condition.name == 'lwp':
+                        if not n.parent.parent in test_nodes_warning:
+                            test_nodes_warning += [n.parent.parent]
+            flag_bio_check = None
+            for n in range(len(test_nodes_warning)):
+                test_node = test_nodes_warning[n]
+                # Create symbol used to record occurrence of excessive
+                # tracer values
+                test_sym = st.new_symbol('ifail'+str(n),symbol_type=DataSymbol,datatype=ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE))
+                # Replace IfBlock that represents the test for excessive values by new
+                # IfBlock that i) detects occurrences of excessive values and ii) no
+                # longer contains any nodes were of IfBlock or Call; if a node of type
+                # Call was removed (assumed to be a halt instruction), then execute the
+                # normal loop content that follows the test conditional on not meeting
+                # the original test condition (i.e., prevent the execution of the loop
+                # content in local conditions that previously would have halted the
+                # model execution
+                test_copy = test_node.copy()
+                ifcond = test_node.condition.copy()
+                ifthen = []
+                warn_nodes = []
+                has_stop = False
+                has_grid_loc = False
+                for node in test_node.if_body.children.copy():
+                    # Remember warning/halt statements
+                    if (len(node.walk(Call)) > 0) or (isinstance(node, IfBlock) and node.condition.walk(Reference)[0].name == 'lwp'):
+                        # Remember current failure counter in specific case for futher code
+                        # modification below
+                        for c in node.walk(Call):
+                            if c.routine.name == 'trc_bio_exceptionnal_fix':
+                                flag_bio_check = 'ifail'+str(n)
+                        if isinstance(node, Call) and node.walk(Call)[0].routine.name == 'ctl_stop':
+                            has_stop = True
+                        warn_nodes += [node.detach()]
+                    else:
+                        ifthen += [node.detach()]
+                ifelse = []
+                if has_stop:
+                    for node in test_node.parent.children[test_node.position+1:]:
+                        if not node in test_nodes_warning:
+                            test_node.parent.children.remove(node)
+                            ifelse += [node]
+                test_node_new = IfBlock.create(ifcond, ifthen, ifelse)
+                test_node.replace_with(test_node_new)
+                test_node = test_node_new
+                # The occurrence of excessive values is counted rather
+                # than flagged in order to permit the compiler to implement
+                # a reduction instead of avoiding parallelisation when
+                # preparing the loop for execution on the GPU
+                #
+                # Remove previous debug output and halt instruction, but
+                # keep a copy for later use (the last node in this list
+                # is the actual halt instruction)
+                counter = Assignment.create(Reference(test_sym), BinaryOperation.create(BinaryOperation.Operator.ADD, Reference(test_sym), Literal('1', INTEGER_TYPE)))
+                test_node.children[1].addchild(counter, index=0)
+                # Seek outside and interior of encompassing horizontal loop
+                # and count the steps taken
+                loop_node = test_node
+                while not ( isinstance(loop_node, NemoLoop) and loop_node.loop_type == 'lat' ):
+                    if isinstance(loop_node, NemoLoop) and loop_node.loop_type == 'lon':
+                        loop_node_interior = loop_node.loop_body
+                    loop_node = loop_node.parent
+                # Save a copy of the loop found
+                loop_copy = loop_node.copy()
+                # Identify parent of copy of test_node
+                npos = loop_node.walk(Node).index(test_node)
+                loop_copy_interior = loop_copy.walk(Node)[npos].parent
+                # Remove all children, and attach the nodes removed from the
+                # original IfBlock
+                loop_copy_interior.pop_all_children()
+                loop_copy_interior.addchild(IfBlock.create(test_copy.condition.copy(), warn_nodes, []))
+                # Tag loop interior as NemoKernel to trigger kernels region detection
+                # of the transformed loop Include initialisation of excessive-value
+                # counter
+                loop_node_interior.addchild(NemoKern(loop_node_interior.pop_all_children()))
+                # Add content of original test statement surrounded by the original
+                # loop in separate test statement outside of horizontal loop
+                counter0 = Assignment.create(Reference(test_sym), Literal('0', INTEGER_TYPE))
+                loop_node.parent.addchild(counter0, index=loop_node.position)
+                new_node = IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.GE, Reference(test_sym), Literal('1', INTEGER_TYPE)), [loop_copy])
+                loop_node.parent.addchild(new_node, index=loop_node.position+1)
+            # Special case for invoke 'trc_bio_medusa' where the detection of
+            # excessive values *in the global domain* should trigger a call to
+            # subroutine 'trc_bio_check' that repeats the same check in the
+            # global domain and provides additional diagnostics
+            #
+            # Note, the flag is duplicated before the global reduction in order
+            # to keep the local version of the flag available
+            if flag_bio_check and 'trc_bio_medusa' in psy.invokes.names:
+                for ib in [node for node in sched.walk(IfBlock) if node.condition.walk(Reference)[0].name == flag_bio_check]:
+                    ig = st.new_symbol('ifail_global',symbol_type=DataSymbol,datatype=ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE))
+                    new_ib = IfBlock.create(Reference(pst.lookup('lk_mpp')), [Call.create(pst.lookup('mpp_max'),[Literal('trcbio_medusa', CHARACTER_TYPE), Reference(ig)])])
+                    pos = ib.position
+                    ib.parent.addchild(Assignment.create(Reference(ig),Reference(st.lookup(flag_bio_check))), index=pos)
+                    ib.parent.addchild(new_ib, index=pos+1)
+                    tbc = None
+                    for c in [call for call in sched.walk(Call) if call.routine.name == 'trc_bio_check']:
+                        tbc = c
+                    if tbc:
+                        tbc.detach()
+                        new_ib = IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.GE, Reference(ig), Literal('1', INTEGER_TYPE)), [tbc])
+                        ib.parent.addchild(new_ib, index=pos+1)
+
+        # ######## AC6 ########
+        # Special treatment of carbonate-chemistry component MOCSY
+        if match('^mocsy_.*', module_name):
+            # Replace options character strings by integer values
+            snames = list(st.symbols_dict.keys()) # static copy of symbol table contents
+            opt_varnames = []
+            for sname in [name for name in snames if name.startswith('op')]:
+                opt_sym = st.lookup(sname)
+                if not isinstance(opt_sym.datatype, UnknownFortranType):
+                    opt_sym.datatype = ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE)
+                else:
+                    repl = match('^CHARACTER[^,:]*(.*)$',opt_sym.datatype.declaration)
+                    if repl:
+                        opt_sym.datatype.declaration = 'INTEGER'+repl.group(1)
+                opt_varnames += [sname]
+            for node in sched.walk((Assignment,CodeBlock)):
+                optval_node = None
+                # Assignments of option variables with option strings
+                if isinstance(node,Assignment) and node.lhs.name in opt_varnames and isinstance(node.rhs, Literal):
+                    optval_node = node.rhs
+                # Comparisons of option variables with option strings
+                elif isinstance(node,CodeBlock) and isinstance(node.ast, Fortran2003.Intrinsic_Function_Reference) and \
+                                                    str(node.ast.children[0]).lower() == 'trim' and \
+                                                    str(node.ast.children[1]).lower() in opt_varnames and \
+                                                    isinstance(node.parent,BinaryOperation) and \
+                                                    node.parent.operator == BinaryOperation.Operator.EQ:
+                    optval_node = node.parent.children[1]
+                    node.replace_with(Reference(st.lookup(str(node.ast.children[1]).lower())))
+                # Remove remaining 'trim' functions applied to option variables
+                elif isinstance(node,CodeBlock) and isinstance(node.ast, Fortran2003.Intrinsic_Function_Reference) and \
+                                                    str(node.ast.children[0]).lower() == 'trim' and \
+                                                    str(node.ast.children[1]).lower() in opt_varnames:
+                    node.replace_with(Reference(st.lookup(str(node.ast.children[1]).lower())))
+                # Catch-all for remaining 'trim' intrinsics in CodeBlocks:
+                # replace whole content of such CodeBlocks with 'STOP'
+                # statement
+                elif isinstance(node, CodeBlock):
+                    for an in node.get_ast_nodes:
+                        for ifr in [ref for ref in walk(an,Fortran2003.Intrinsic_Function_Reference) if str(ref.items[0]).lower() == 'trim']:
+                            node.replace_with(CodeBlock([Fortran2003.Stop_Stmt('STOP')],CodeBlock.Structure.STATEMENT))
+                if optval_node:
+                    # Normalise name
+                    name = optval_node.value.replace('/','').lower()
+                    val = opt_hashval(name)
+                    optval_node.replace_with(Literal(str(val),INTEGER_TYPE))
+            # Replace argument 'n' (if available) by integer parameter of value
+            # 1 in order to prevent the use of automatic arrays of shape
+            # 'dimension(n)' (if any)
+            # Note, in the unmodified source code, a value of 1 is passed to
+            # dummy argument 'n' of subroutine 'mocsy_interface' (which is used
+            # to set up the automatic arrays in this subroutine), but
+            # subsequently the a literal of 1 is used to set the value of 'n'
+            # throughout the MOCSY componenent; here we fix the size of the
+            # corresponding arrays to 1, and ignore the value of 'n' supplied
+            # to subroutine 'mocsy_interface'
+            if 'n' in st.symbols_dict.keys():
+                auto_arrays = []
+                for sym in st.symbols:
+                    if isinstance(sym,DataSymbol) and \
+                       len(sym.shape) == 1 and \
+                       not sym.is_argument and \
+                       isinstance(sym.shape[0][1],Reference) and \
+                       sym.shape[0][1].name == 'n':
+                            auto_arrays += [sym]
+                if len(auto_arrays) > 0:
+                    n_name = 'n_'+invoke.name+'_'+module_name
+                    pst.new_symbol(n_name,symbol_type=DataSymbol,datatype=ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE),constant_value=1)
+                    for auto_array in auto_arrays:
+                        auto_array.datatype = ArrayType(auto_array.datatype,shape=[ArrayType.ArrayBounds(lower=1,upper=Reference(pst.lookup(n_name)))])
+
+#############################################################################################
+# Section 2 - Exclusions of modules from generic transformations and kernels-region detection
+#############################################################################################
+
+    # No acceleration in specified modules
+    if any([match(module_re, module_name) for module_re in ACC_MODULES_EXCLUDE_RE]):
+            return psy
+
+    # Loop through all invokes
+    for invoke in psy.invokes.invoke_list:
+        sched = invoke.schedule
+        st = sched.symbol_table
+        pst = st.parent_symbol_table()
+
+        # Skip exlucded invokes
+        if not invoke.name in ACC_INVOKES_EXCLUDE:
+
+#################################################
+# Section 3 - Generic source-code transformations
+#################################################
+
+            # ######## AC20 #######
+            #
+            # Shield all accessible (i.e., not inside CodeBlocks) calls of the
+            # model-data output subroutine 'iom_put' using an IfBlock to avoid
+            # array accesses or array computations inside the dummy argument
+            # list of such calls when model-data output via the IOM interface
+            # is disabled (i.e., compilation without pre-processor key
+            # 'key_iom_put'); this is a generalisation of the practice of
+            # commenting out unshielded 'iom_put' calls selectively in
+            # PSyclone-processed source code to prevent unnecessary array
+            # accesses by the host when using managed-memory GPU acceleration
+            # (C. Dearden, pers. comm.)
+            replace = []
+            # Schedule calls for shielding
+            for c in [call for call in sched.walk(Call) if call.routine.name == 'iom_put']:
+                replace += [c]
+            if replace:
+                # Create 'lk_iomput' symbol import where missing (not or
+                # implicitely imported)
+                try:
+                    sym = st.lookup('lk_iomput')
+                except KeyError:
+                    try:
+                        csym = st.lookup('iom')
+                    except KeyError:
+                        csym = psy.new_symbol('iom', symbol_type=ContainerSymbol)
+                    sym = DataSymbol('lk_iomput', datatype=BOOLEAN_TYPE, interface=ImportInterface(csym))
+                    if not csym.wildcard_import:
+                        st.add(sym)
+                # Shield calls
+                for c in replace:
+                    c.replace_with(IfBlock.create(Reference(sym), [c.copy()]))
+
+            # ##### AC10/AC13 #####
+            # ##### AC16/AC19 #####
+            #
+            # Replacement of automatic arrays by module arrays in all modules that
+            # will be subject to kernels region detection
+            sym_auto = []
+            # Find automatic arrays: local, non-argument, non-constant arrays
+            # with defined array ranges; rule out arrays of derived-type variables
+            for sym in st.local_datasymbols:
+                # Rule out transformation of arrays that are the return values
+                # of functions (AC16)
+                if sym.name.lower() == invoke.name.lower():
+                    continue
+                deferred = False
+                fixed = False
+                if sym.is_local and not sym.is_argument and not sym.is_constant and sym.is_array and not isinstance(sym.datatype.intrinsic, DataTypeSymbol):
+                    for dim in sym.shape:
+                        if not isinstance(dim, ArrayType.ArrayBounds):
+                            deferred = True
+                        elif ( isinstance(dim.lower, Literal) or ( isinstance(dim.lower, Reference) and dim.lower.symbol.is_constant ) ) and \
+                             ( isinstance(dim.upper, Literal) or ( isinstance(dim.upper, Reference) and dim.upper.symbol.is_constant ) ):
+                            fixed = True
+                    if not deferred and not fixed:
+                        sym_auto += [sym]
+            # Extract ranges used for automatic arrays, and create corresponding
+            # allocatable module arrays
+            sym_dict = {}
+            for sym in sym_auto:
+                rngs = []
+                for dim in sym.shape:
+                    rng = []
+                    for bnd in [dim.lower, dim.upper]:
+                        if isinstance(bnd, Literal):
+                            rng += [bnd.value]
+                        else:
+                            rng += [bnd.name]
+                    rngs += [':'.join(rng)]
+                sym_mod = pst.new_symbol(sym.name+'_'+invoke.name, symbol_type=DataSymbol, \
+                                         datatype=ArrayType(ScalarType(sym.datatype.intrinsic, sym.datatype.precision), shape=[ArrayType.Extent.DEFERRED]*len(rngs)))
+                sched.addchild(node_allocate(sym_mod, rngs), index=0)
+                sym_dict[sym.name] = sym_mod.name
+                # If precision of the variable is imported (e.g., 'use
+                # par_kind, only: wp') and if the import statement is present
+                # at procedure and not at module level, replicate the
+                # corresponding use association at the module scope (assume
+                # that if the precision symbol is imported at the procedure
+                # scope, it is not imported in a blanked module import and it
+                # is not imported at the module scope)
+                if isinstance(sym.datatype.precision, DataSymbol) and sym.datatype.precision.is_import and \
+                   sym.datatype.precision.name in st.symbols_dict.keys() and not sym.datatype.precision.name in pst.symbols_dict.keys():
+                    psym = sym.datatype.precision
+                    csym = pst.new_symbol(psym.interface.container_symbol.name, symbol_type=ContainerSymbol)
+                    pst.new_symbol(psym.name, symbol_type=DataSymbol, datatype=psym.datatype, interface=ImportInterface(csym))
+            # Redirect all references (incl. references inside CodeBlock
+            # objects, possibly with multiple AST nodes) to automatic arrays to
+            # point to the newly created corresponding module arrays; further
+            # tests with source-code that is not currently processed revealed
+            # that the capitalisation of references in code blocks can differ
+            # form that used in the PSyIR, and therefore are noramlised here 
+            for ref in sched.walk((Reference,CodeBlock)):
+                if isinstance(ref, Reference) and ref.name in sym_dict.keys():
+                    ref.symbol = pst.lookup(sym_dict[ref.name])
+                elif isinstance(ref, CodeBlock):
+                    for ast in ref.get_ast_nodes:
+                        for name in walk(ast, Fortran2003.Name):
+                            if name.tostr().lower() in sym_dict.keys():
+                                name.string = sym_dict[name.tostr().lower()]
+            # Rename the now-unused automatic arrays (remove method for data
+            # symbols is not implemented)
+            for sym in sym_auto:
+                st.rename_symbol(sym, 'UNUSED___'+sym.name)
+
+            # ######## AC16 ####### (part 1)
+            #
+            # Prevent three types of function calls that are incorrectly
+            # identified as ArrayReferences to be included in kernels regions:
+            #   1) function calls with array arguments;
+            #   2) function calls of functions defined in the parent module;
+            #      and
+            #   3) function calls with ArrayOfStructuresReference arguments
+            #
+            # Note, as the function calls are packed into CodeBlock objects,
+            # they are no longer directly available in the PSyIR tree; as a
+            # consequence, further transformations that inquire or modify these
+            # function calls may no longer fully work (such as a the transformation
+            # of automatic arrays); therefore such CodeBlock objects should
+            # be created just before the kernels region detection, and then
+            # reverted back
+            #
+            # Pack into CodeBlock objects references that are either
+            # RoutineSymbol objects in the module's symbol table (i.e., module
+            # functions), or ArrayReferences with children that contain Range
+            # nodes and are either ArrayReference or
+            # ArrayOfStructuresReferences; these CodeBlocks will not produce
+            # the intended output, and should be reverted prior to converting
+            # the PSyIR to source-code output
+            for r in sched.walk(Reference):
+                replace = False
+                if ( r.name.lower() in pst.symbols_dict.keys() and not r.name.lower() == invoke.name.lower() ) and isinstance(pst.lookup(r.name.lower()), RoutineSymbol):
+                    replace = True 
+                if isinstance(r, ArrayReference):
+                    for arc in [c for c in r.children if c.walk((ArrayReference, ArrayOfStructuresReference)) and c.walk(Range)]:
+                        replace = True
+                if replace:
+                    r.replace_with(CodeBlock([r], CodeBlock.Structure.EXPRESSION))
+
+##########################################
+# Section 4 - Detection of kernels regions
+##########################################
+
+            # List of all kernels and empty list of nodes taken care of
+            kernels = sched.kernels()
+            # Add to kernels list the list of all array assignments
+            # (implicit loops), but exclude those that contain CodeBlocks
+            for a in [node for node in sched.walk(Assignment) if node.is_array_range and not node.walk(CodeBlock)]:
+                kernels += [a]
+            # Exclude vertical loops (NemoLoop objects with
+            # NemoLoop.loop_type == 'levels') from kernels regions for some
+            # invokes
+            vloops = True
+            if invoke.name in ACC_INVOKES_EXCLUDE_VERTICAL_LOOPS + acc_invokes_exclude_vertical_loops_extra:
+                vloops = False
+            # Exclude if constructs in kernels regions for some invokes
+            exclude_if = False
+            if invoke.name in ACC_INVOKES_EXCLUDE_IFBLOCKS:
+                exclude_if = True
+            # Loop through kernels not yet taken care of
+            for kernel in kernels:
+                if not kernel.ancestor(ACCKernelsDirective) and not kernel.walk(ACCKernelsDirective):
+                    # Go up by as many levels as possible
+                    kernel_up = True
+                    while (kernel_up):
+                        kernel_up = kernels_region_grow(kernel, direction=0, accept_loop_levels=vloops, exclude_if=exclude_if)
+                        if kernel_up:
+                            kernel = kernel_up
+                    # Increase scope of kernels region as widely as possible
+                    kernels_region = [kernel]
+                    kernel_wide = kernel
+                    while (kernel_wide):
+                        kernel_wide = kernels_region_grow(kernel_wide, direction=1, accept_loop_levels=vloops, exclude_if=exclude_if)
+                        if kernel_wide:
+                            kernels_region += [kernel_wide]
+                    kernel_wide = kernel
+                    while (kernel_wide):
+                        kernel_wide = kernels_region_grow(kernel_wide, direction=-1, accept_loop_levels=vloops, exclude_if=exclude_if)
+                        if kernel_wide:
+                            kernels_region = [kernel_wide] + kernels_region
+                    # Add "!$acc kernels" directive around identified
+                    # group of nodes *unless* it is a solitary
+                    # InlinedKern (i.e. that hasn't "grown" and thus been
+                    # tested for unwanted content) node and mark all
+                    # kernels and assignments as taken care of
+                    if isinstance(kernels_region[0], (NemoLoop, IfBlock, Assignment)):
+                        ACCKernelsTrans().apply(kernels_region)
+                    # Add "!$acc loop independent collapse(2)" directives
+                    # around 2-dimensional loops that span the horizontal
+                    # plane (i.e., nested loops of types 'lat' and 'lon',
+                    # respectively)
+                    for kernel in kernels_region:
+                        # Detect 'lat' loops
+                        loop_lat = [loop for loop in kernel.loops() if loop.loop_type == 'lat']
+                        # Detect 'lon' loops
+                        loop_lon = [loop for loop in kernel.loops() if loop.loop_type == 'lon']
+                        # Detect 'lat' loops which contain exactly one node that is a 'lon' loop
+                        for loop in [l for l in loop_lat if (len(l.loop_body.children) == 1) and (l.loop_body.children[0] in loop_lon)]:
+                            # Directives to fully parallelise and
+                            # collapse 2-dimensional loop
+                            ACCLoopTrans().apply(loop, {'independent': True, 'collapse': 2})
+            
+###################################
+# Section 5 - General finalisations
+###################################
+
+            # ######## AC16 ####### (part 2)
+            #
+            # Unpack nodes that have previously been packed into CodeBlock
+            # objects
+            for cb in [node for node in sched.walk(CodeBlock) if isinstance(node.ast, Node)]:
+                cb.replace_with(cb.ast)
+            
+####################################
+# Section 6 - Specific finalisations
+####################################
+
+        # ######## AC22 ####### (part 1)
+        #
+        # Special treatment of 'tra_nxt_vvl' to ensure maximum loop
+        # collapse inside kernels regions in order to avoid partial sequential
+        # execution of nested loop
+        if invoke.name == 'tra_nxt_vvl':
+            for kr in sched.walk(ACCKernelsDirective):
+                # Remove any pre-existing loop directives
+                for ald in kr.walk(ACCLoopDirective):
+                    ald.replace_with(ald.walk(NemoLoop)[0].detach())
+                # Find groups of directly nested loops and their depth
+                loops0 = kr.walk(NemoLoop)
+                loops = {}
+                while loops0:
+                    # Process first loop in list
+                    l = loops0[0]
+                    # Move down to the innermost level of directly nested loops
+                    # (if any)
+                    while(isinstance(l.loop_body.children[0], NemoLoop) and
+                            len(l.loop_body.children)==1):
+                        l = l.loop_body.children[0]
+                    depth = 1
+                    # Move up to the outermost level of directly nested loops
+                    # (if any) and remove them from the loop-nest candidate
+                    # list
+                    while(isinstance(l.parent.parent, NemoLoop)):
+                        loops0.remove(l)
+                        depth += 1
+                        l = l.parent.parent
+                    loops0.remove(l)
+                    # Store loop nest and its depth
+                    if depth > 1:
+                        loops[l] = depth
+                # Add new loop directives
+                for l in loops.keys():
+                    ACCLoopTrans().apply(l, {'independent': True, 'collapse': loops[l]})
+
+        # ######## AC16 ####### (part 3)
+        #
+        # Specific treatment of statement function ('fintegral')
+        # in one invoke ('sbc_dcy') to ensure it appears at the top of the
+        # schedule (it may have been displaced by the automatic-array
+        # transformation); statement functions are not currently explicitelly
+        # recognised in the PSyIR, but rather represented in the form of
+        # Assignment nodes, which are affected by transformations like other
+        # nodes
+        if invoke.name == 'sbc_dcy':
+            for a in [node for node in sched.walk(Assignment) if node.lhs.name == 'fintegral']:
+                a.parent.addchild(a.detach(), index=0)
+
+    return psy
+
Index: NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/psyclone_transformation_MEDUSA.py
===================================================================
--- NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/psyclone_transformation_MEDUSA.py	(revision 15762)
+++ NEMO/branches/NERC/dev_r4.0.4_MEDUSA_PSyclone/mk_psyclone/psyclone_transformation_MEDUSA.py	(revision 15762)
@@ -0,0 +1,710 @@
+#======================================================================
+#
+# PSyclone transformation script for the PSyclone processing of the
+# MEDUSA component of NEMO-MEDUSA model setups
+#
+#======================================================================
+# History : 4.0 ! 2022 (S. Mueller) Original version
+#======================================================================
+
+#----------------------------------------------------------------------
+# NEMO 4.0 , NEMO Consortium (2022)
+# $Id$
+# Software governed by the CeCILL license (see ./LICENSE)
+#----------------------------------------------------------------------
+
+
+from psyclone.transformations         import ACCKernelsTrans, ACCLoopTrans
+from psyclone.psyir.nodes             import CodeBlock, IfBlock, Assignment, Call, Reference, Literal, \
+                                             BinaryOperation, UnaryOperation, ArrayReference, Range, Loop, Node
+from psyclone.psyir.symbols           import DataSymbol, ContainerSymbol, ImportInterface, RoutineSymbol
+from psyclone.psyir.symbols           import ScalarType, UnknownFortranType, ArrayType
+from psyclone.psyir.symbols           import INTEGER_TYPE, CHARACTER_TYPE
+from psyclone.psyir.frontend.fparser2 import Fparser2Reader
+from psyclone.nemo                    import NemoLoop, NemoKern
+from fparser.two                      import Fortran2003
+from fparser.two.utils                import walk
+from re                               import match
+
+# Exclude a range from modules from acceleration (regular expressions)
+ACC_EXCLUDE_MODULE_RE = [ '^mocsy_.*' ]
+
+# Exclude a range of invokes from accelerated code as the resulting ACC
+# kernels fail to be compiled
+ACC_EXCLUDE_INVOKES = [ 'trc_bc',               # trcbc.F90
+                        'trc_dmp_clo',          # trcdmp.F90
+                        'trc_nam_ice',          # trcice.F90
+                        'p4z_ini',              # trcini_pisces.F90
+                        'p2z_ini',              # trcini_pisces.F90
+                        'prt_ctl_trc'           # prtctl_trc.F90, avoid assignment of character array withing kernels region
+                      ]
+
+# Exclude code that is conditional to a single logical (e.g., debug
+# output, tracer-trend diagnostics) from accelerated code
+ACC_REMOVE_CONDITIONAL_CODE = [ 'ln_msa_debug',   # MEDUSA debug output
+                                'l_trdtrc'        # Tracer-trend diagnostic
+                              ]
+
+# Permit some subroutine calls inside accelerated kernels regions
+ACC_PERMIT_SUBROUTINE_CALLS = [ 'trc_oxy_medusa', 'oxy_sato', 'oxy_schmidt', 'trc_dms_medusa', 'dms_flux_ocn', 'btd_gastrans', 'mocsy_interface', 'ddpdd' ]
+
+# Select refined timing of MEDUSA component
+TIMING_MEDUSA_REFINED = True
+
+def widen(kernel, dir=1):
+    pos = kernel.parent.children.index(kernel)
+    if (pos+dir >= 0) and (pos+dir < len(kernel.parent.children)):
+        kernel = kernel.parent.children[pos+dir]
+        # Must be either NemoLoop or IfBlock or Assignment to accept, but must not include CodeBlock, Call
+        if not isinstance(kernel, (NemoLoop,IfBlock,Assignment)) or len(kernel.walk(CodeBlock)) > 0 or has_excluded_call(kernel):
+            kernel = False
+    else:
+        kernel = False
+    return kernel
+
+def up(kernel):
+    if (kernel.parent):
+        kernel = kernel.parent.parent
+        # Must be either NemoLoop or IfBlock or Assignment to accept, but must not include CodeBlock, Call
+        if not isinstance(kernel, (NemoLoop,IfBlock,Assignment)) or len(kernel.walk(CodeBlock)) > 0 or has_excluded_call(kernel):
+            kernel = False
+    else:
+        kernel = False
+    return kernel
+
+# Test for presence of unallowed call in kernels region 
+def has_excluded_call(kernel):
+    for call in kernel.walk(Call):
+        if not call.routine.name in ACC_PERMIT_SUBROUTINE_CALLS:
+            return True
+    return False
+
+# Create allocation statement for variable represented by symbol 'sym' of the
+# shape supplied in argument 'shape' (if list of stings) or of the
+# shape of the symbol named in argument 'shape' (if string); depending on
+# 'cond', make allocation dependent on the allocation status of the variable.
+def node_allocate(sym, shape, cond=True):
+    if sym.is_array:
+        rank = len(sym.shape)
+    elif isinstance(sym.datatype, UnknownFortranType):
+        rnge = match('^.*\(([ :,]*)\).+$',sym.datatype.declaration)
+        if rnge:
+            rank = rnge.group(1).count(':')
+    else:
+        return None
+    if isinstance(shape, str):
+        shape_ref = shape
+        shape = []
+        for n in range(rank):
+            shape += ['LBOUND('+shape_ref+','+str(n+1)+'):UBOUND('+shape_ref+','+str(n+1)+')']
+    allocate = CodeBlock([Fortran2003.Allocate_Stmt('allocate('+sym.name+'('+','.join(shape)+'))')], CodeBlock.Structure.STATEMENT)
+    if cond:
+        allocatd = CodeBlock([Fortran2003.Intrinsic_Function_Reference('allocated('+sym.name+')')], CodeBlock.Structure.EXPRESSION)
+        allocate = IfBlock.create(UnaryOperation.create(UnaryOperation.Operator.NOT, allocatd), [allocate])
+    return allocate
+
+# Transformation
+def trans(psy):
+
+    # Module name
+    module_name = psy.name[4:-4]
+
+##################################################
+# Section 1 - Specific source-code transformations
+##################################################
+
+    # Loop through all invokes
+    for invoke in psy.invokes.invoke_list:
+        sched = invoke.schedule
+        st = sched.symbol_table
+        pst = st.parent_symbol_table()
+
+        # Refined MEDUSA timing
+        #
+        # Additional timing instructions that refine 'trc_bio_medusa' section
+        if TIMING_MEDUSA_REFINED and invoke.name == 'trc_bio_medusa':
+            labels = ['trc_bio_medusa_ini', 'trc_bio_msa_nan', 'trc_bio_msa_atm', \
+                      'trc_bio_msa_cch',    'trc_bio_msa_vlp', 'trc_bio_msa_trc', \
+                      'trc_bio_msa_mip',    'trc_bio_msa_chk', 'trc_bio_msa_as',  \
+                      'trc_bio_msa_pln',    'trc_bio_msa_isc', 'trc_bio_msa_det', \
+                      'trc_bio_msa_upd',    'trc_bio_msa_dia', 'trc_bio_msa_ben'  \
+                      ]
+            for n in range(len(labels)):
+                # Find parent nodes (where they change) and positions (incl.
+                # anticipation of position change following addition of start
+                # statement) for the addition of the various timing statements
+                if n==0:
+                    snode = [ib for ib in sched.walk(IfBlock) if isinstance(ib.condition, Reference) and ib.condition.name=='ln_timing'][0]
+                    sp    = snode.parent
+                    spos  = snode.position+1
+                    enode = [c for c in sched.walk(Call) if c.routine.name=='btd_ini'][0]
+                    ep    = enode.parent
+                    epos = enode.position+2
+                elif n==1:
+                    spos = epos+1
+                    epos = [ib for ib in sched.walk(IfBlock) if isinstance(ib.condition, Reference) and ib.condition.name=='lk_oasis'][0].position+1
+                elif n==2:
+                    spos = epos+1
+                    for ib in sched.walk(IfBlock):
+                        if [r for r in ib.walk(Reference) if r.name == 'ln_rsttr']:
+                            epos = ib.position+2
+                            break
+                elif n in [3,4]:
+                    spos = epos+2
+                    epos = spos+2
+                elif n==5:
+                    sp = sp.children[spos+1].loop_body
+                    spos = 0
+                    ep = sp
+                    epos = [ib for ib in sp.walk(IfBlock) if isinstance(ib.condition, Reference) and ib.condition.name == 'ln_msa_omipdic'][0].position+1
+                elif n==6:
+                    spos = epos+1
+                    epos = [c for c in sp.walk(Call) if c.routine.name == 'trc_bio_check'][0].position+2
+                elif n==7:
+                    spos = epos+1
+                    epos = [c for c in sp.walk(Call) if c.routine.name == 'btd_air_sea'][0].ancestor(IfBlock).position+1
+                elif n<13:
+                    spos = epos+1
+                    epos = spos+2
+                elif n==13:
+                    spos = epos+1
+                    epos = spos+4
+                elif n==14:
+                    sp   = vlte.parent
+                    spos = vlte.position+1
+                    ep   = sp
+                    epos = len(ep.children)
+                # Add start statement
+                ts = IfBlock.create(Reference(st.lookup('ln_timing')), [Call.create(RoutineSymbol('timing_start'), [Literal(labels[n], CHARACTER_TYPE)])])
+                sp.addchild(ts, index=spos)
+                # Add stop statement
+                te = IfBlock.create(Reference(st.lookup('ln_timing')), [Call.create(RoutineSymbol('timing_stop'), [Literal(labels[n], CHARACTER_TYPE)])])
+                ep.addchild(te, index=epos)
+                # Remeber where the main vertical loop ends
+                if n == 4:
+                    vlte = te
+
+        # ######## AC12 #######
+        #
+        # Reduce number of times the a halo exchange (procedure 'lbc_lnk' call)
+        # that are made during the MEDUSA time step: only call 'lbc_lnk' at the
+        # deepest depth layer
+        if invoke.name == 'btd_phyto_pkn':
+            for c in sched.walk(Call):
+                if c.routine.name == 'lbc_lnk':
+                    if 'jk' in st.symbols_dict.keys() and ('jpj' in st.symbols_dict.keys() or 'jpk' in st.symbols_dict.keys()):
+                        if not 'jpk' in st.symbols_dict.keys():
+                            st.new_symbol('jpk', symbol_type=DataSymbol, datatype=INTEGER_TYPE).copy_properties(st.lookup('jpj'))
+                        sym_jpk = st.lookup('jpk')
+                        pos = c.position
+                        c.parent.addchild(IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.EQ, Reference(st.lookup('jk')), Reference(st.lookup('jpk'))), [c.detach()]), index=pos)
+
+        # ######## AC2 ########
+        #
+        # Selective removal of conditional code segments
+        # Note, only removes segements that are conditional on
+        # single-statement logical
+        for ifb in invoke.schedule.walk(IfBlock):
+            if isinstance(ifb.condition, Reference) and (ifb.condition.name in ACC_REMOVE_CONDITIONAL_CODE):
+                ifb.parent.children.remove(ifb)
+
+        # ######## AC8 ########
+        # Special treatment of the glob_sum function call in subroutine 'trc_stp'
+        # of module 'trcstp'
+        #
+        # Extract the multiplication of 3-dimensional arrays from the argument of a
+        # function call, so that it can be automatically placed inside a kernels
+        # region
+        if module_name == 'trcstp' and invoke.name == 'trc_stp':
+            # Create temporary array
+            sym = pst.new_symbol('zfn_'+sched.invoke.name, symbol_type=DataSymbol, datatype=ArrayType(ScalarType(ScalarType.Intrinsic.REAL, ScalarType.Precision.DOUBLE), shape=[ArrayType.Extent.DEFERRED]*3))
+            sched.children.insert(0, node_allocate(sym, 'cvol'))
+            # Separate arithmetic operations within the argument list of function calls
+            # from the function calls to enable acceleration of these operations
+            #
+            # Assume that the expression consists of an assignment of the result of the
+            # function call to a variable; the function call is represented by a
+            # CodeBlock
+            #
+            # Look for function call (i.e., Fortran2003.Structure_Constructor instance,
+            # with Fortran2003.Type_Name child) with arithmetic-expression argument
+            # (Fortran2003.Add_Operand instance present)
+            asgn_node = None
+            for asgn_node in sched.walk(Assignment):
+                for cb in asgn_node.walk(CodeBlock):
+                    for fn_ast in walk(cb.ast, Fortran2003.Structure_Constructor):
+                        if walk(fn_ast.children[0], Fortran2003.Type_Name) and walk(fn_ast.children[1], Fortran2003.Add_Operand):
+                            # AST of function-call argument
+                            arg_ast = walk(fn_ast.children[1], Fortran2003.Add_Operand)
+                            # Create assignment of arithmetic expression found in
+                            # argument, and make sure full array ranges are
+                            # explicitely specified so that the expression can be
+                            # detected as the seed for a kernels region
+                            rnge = []
+                            for n in range(len(sym.shape)):
+                                rnge += [Range.create(BinaryOperation.create(BinaryOperation.Operator.LBOUND, Reference(sym), Literal(str(n+1), INTEGER_TYPE)), BinaryOperation.create(BinaryOperation.Operator.UBOUND, Reference(sym), Literal(str(n+1), INTEGER_TYPE)))]
+                            asgn_new = Assignment.create(ArrayReference.create(sym, rnge), CodeBlock(arg_ast, CodeBlock.Structure.EXPRESSION))
+                            # Include assignment node in PSyIR
+                            asgn_node.parent.children.insert(asgn_node.parent.children.index(asgn_node), asgn_new)
+                            # Translate AST for arithmetic expression into PSyIR sub-tree
+                            asgn_new.children[1] = Fparser2Reader()._binary_op_handler(asgn_new.rhs.ast, asgn_new)
+                            # Replace original argument by name of temporary variable
+                            cb.replace_with(CodeBlock([fn_ast.tostr().replace(arg_ast[0].tostr(), sym.name)], CodeBlock.Structure.EXPRESSION))
+
+        # ######## AC8 ########
+        #
+        # Special treatment of glob_sum_* functions (module
+        # lib_fortran)
+        #
+        # Expand the n-dimensional reduction by reducing one dimension
+        # only in each of the expanded loop nests to expose some
+        # parallelism at the expense of moderate memory requirements
+        #
+        # Note, this section makes a range of assumptions:
+        #  1) The reduction occurs in a nested NemoLoop structure
+        #     where one or more of the outermost loops may be
+        #     irrelevant for the reduction (loop count of 1, relict
+        #     from macro-based code generation)
+        #  2) The reduction operation consists of a) an assignment
+        #     that represents the current summand and b) a call to a
+        #     reduction subroutine with two arguments, of which the
+        #     second is the temporary reduction result
+        #  3) the only dummy argument that is an array represents
+        #     the array to be summed
+        #  4) the horizontal shape of the array to be summed is
+        #     identical in each call of the modified subroutine
+        #     and therefore the newly added module array does only
+        #     have to be allocated in the first call of the subroutine
+        if invoke.name.startswith('glob_sum_'):
+            # The name of the input array
+            for sym in st.symbols:
+                if sym.is_array and sym.is_argument:
+                    in_var = sym.name
+                    ndim = len(sym.shape)
+                    break
+            # The orginal loops of the reduction that are relevant
+            # for the reduction (outermost will by first element in
+            # array)
+            loops = sched.walk(NemoLoop)[-ndim:]
+            nloops = len(loops)
+            loop_indices = []
+            for n in range(len(loops)):
+                loop_indices.append(loops[n].variable)
+            # Only transform loop nests, not single loops
+            if len(loops) > 1:
+                # The name of the reduction-operator subroutine
+                op_proc = loops[0].walk(Call)[0].routine.name
+                # The name of the accumulation variable, and the
+                # associated symbol
+                out_var = loops[0].walk(Call)[0].children[1].name
+                sym_ov = st.lookup(out_var)
+                # Create expansion of accumulation variable as an
+                # allocatable module variable
+                if isinstance(sym_ov.datatype, UnknownFortranType):
+                    pst.new_symbol(out_var+'_'+invoke.name, symbol_type=DataSymbol, datatype=UnknownFortranType(sym_ov.datatype.declaration.replace('::',',dimension('+','.join([':']*(nloops-1))+'),allocatable ::').replace(out_var, out_var+'_'+invoke.name)))
+                    sym_ov.datatype.declaration = sym_ov.datatype.declaration.replace(out_var, 'UNUSED___'+out_var)
+                else:
+                    pst.new_symbol(out_var+'_'+invoke.name, symbol_type=DataSymbol, datatype=ArrayType(sym_ov.datatype, shape=[ArrayType.Extent.DEFERRED]*(len(loops)-1)))
+                st.rename_symbol(sym_ov, 'UNUSED___'+out_var)
+                sym_ov = pst.lookup(out_var+'_'+invoke.name)
+                # Insert allocation statement
+                sched.children.insert(0, node_allocate(sym_ov, in_var))
+                # Create new references with indices '(:,:)', '(ji,jj), '(ji,1)',
+                # '(1,1)' (3-d case) or '(:)', '(ji)', '(1)' (2-d case)
+                refs = []
+                for n in range(-1,nloops):
+                    ref_new = ArrayReference(sym_ov)
+                    for nn in range(nloops-1):
+                        if n < 0:
+                            ref_new.addchild(Range.create(BinaryOperation.create(BinaryOperation.Operator.LBOUND, Reference(sym_ov), Literal(str(nn+1), INTEGER_TYPE)),BinaryOperation.create(BinaryOperation.Operator.UBOUND, Reference(sym_ov), Literal(str(nn+1), INTEGER_TYPE))))
+                        elif ((nloops-2)-nn)-n < 0:
+                            ref_new.addchild(Literal('1', INTEGER_TYPE))
+                        else:
+                            ref_new.addchild(Reference(loop_indices[nloops-nn-1]))
+                    refs.append(ref_new)
+                # Replace existing references to the accumulation variable
+                for a in sched.walk(Reference):
+                    if a.name == 'UNUSED___'+out_var:
+                        if isinstance(a.parent, Assignment):
+                            a.replace_with(refs[0].copy())
+                        elif isinstance(a.parent, Call) and a.parent.routine.name == op_proc:
+                            a.replace_with(refs[1].copy())
+                        else:
+                            a.replace_with(refs[-1].copy())
+                        first = False
+                for n in range(nloops-1,0,-1):
+                    l = loops[n].copy()
+                    if len(l.loops()) < 3:
+                        l.loop_type = 'unknown'
+                    l.walk(Assignment)[0].detach()
+                    l.walk(Call)[0].children[0].replace_with(refs[n].copy())
+                    l.walk(Call)[0].children[1].replace_with(refs[n+1].copy())
+                    l.start_expr = Literal('2', INTEGER_TYPE)
+                    loops[0].parent.children.insert(loops[0].parent.children.index(loops[0])+1,l)
+    
+        # #### AC7/AC9/AC11 ###
+        #
+        # Removal of warning output from the horizontal loops in the air-sea
+        # gas-exchange and carbonate-chemistry update subroutines (invokes
+        # 'btd_air_sea' and 'btd_carbchem')
+        #
+        # With the carbonate-chemistry code resident on the GPU (i.e., calls to
+        # mocsy_interafce accepted in kernels regions), the air-sea gas-exchange
+        # and carbonate-chemistry updates can be fully parallelised with the
+        # exception of the detection of excessive tracer conditions and the the
+        # detection and the associated
+        # ermegency stop.
+        #
+        # Here, we replace the warning output and the emergency stop by a detection
+        # mechanism but delay the warning output or model halt until the nested
+        # loop over all water columns is complete; for simplicity, the output of
+        # location and and tracer
+        # values where excessive values are encountered is suppressed
+        #
+        # Identify all nodes that encompasse the test, debug output,
+        # and, if present, the halt statement
+        test_nodes_warning = []
+        if invoke.name in ['btd_air_sea', 'btd_carbchem', 'trc_bio_medusa']:
+            for l in sched.walk(NemoLoop):
+                if l.loop_type == 'lat':
+                    for n in l.walk((Call, IfBlock)):
+                        if isinstance(n, Call) and n.routine.name == 'ctl_stop':
+                            if not n.parent.parent in test_nodes_warning:
+                                test_nodes_warning += [n.parent.parent]
+                        elif isinstance(n, IfBlock) and isinstance(n.condition, Reference) and n.condition.name == 'lwp':
+                            if not n.parent.parent in test_nodes_warning:
+                                test_nodes_warning += [n.parent.parent]
+            flag_bio_check = None
+            for n in range(len(test_nodes_warning)):
+                test_node = test_nodes_warning[n]
+                # Create symbol used to record occurrence of excessive
+                # tracer values
+                test_sym = st.new_symbol('ifail'+str(n),symbol_type=DataSymbol,datatype=ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE))
+                # Replace IfBlock that represents the test for excessive values by new
+                # IfBlock that i) detects occurrences of excessive values and ii) no
+                # longer contains any nodes were of IfBlock or Call; if a node of type
+                # Call was removed (assumed to be a halt instruction), then execute the
+                # normal loop content that follows the test conditional on not meeting
+                # the original test condition (i.e., prevent the execution of the loop
+                # content in local conditions that previously would have halted the
+                # model execution
+                test_copy = test_node.copy()
+                ifcond = test_node.condition.copy()
+                ifthen = []
+                warn_nodes = []
+                has_stop = False
+                has_grid_loc = False
+                for node in test_node.if_body.children.copy():
+                    # Remember warning/halt statements
+                    if (len(node.walk(Call)) > 0) or (isinstance(node, IfBlock) and node.condition.walk(Reference)[0].name == 'lwp'):
+                        # Remember current failure counter in specific case for futher code
+                        # modification below
+                        for c in node.walk(Call):
+                            if c.routine.name == 'trc_bio_exceptionnal_fix':
+                                flag_bio_check = 'ifail'+str(n)
+                        if isinstance(node, Call) and node.walk(Call)[0].routine.name == 'ctl_stop':
+                            has_stop = True
+                        warn_nodes += [node.detach()]
+                    else:
+                        ifthen += [node.detach()]
+                ifelse = []
+                if has_stop:
+                    for node in test_node.parent.children[test_node.parent.children.index(test_node)+1:]:
+                        if not node in test_nodes_warning:
+                            test_node.parent.children.remove(node)
+                            ifelse += [node]
+                test_node_new = IfBlock.create(ifcond, ifthen, ifelse)
+                test_node.replace_with(test_node_new)
+                test_node = test_node_new
+                # The occurrence of excessive values is counted rather
+                # than flagged in order to permit the compiler to implement
+                # a reduction instead of avoiding parallelisation when
+                # preparing the loop for execution on the GPU
+                #
+                # Remove previous debug output and halt instruction, but
+                # keep a copy for later use (the last node in this list
+                # is the actual halt instruction)
+                counter = Assignment.create(Reference(test_sym), BinaryOperation.create(BinaryOperation.Operator.ADD, Reference(test_sym), Literal('1', INTEGER_TYPE)))
+                test_node.children[1].children.insert(0,counter)
+                # Seek outside and interior of encompassing horizontal loop
+                # and count the steps taken
+                loop_node = test_node
+                while not ( isinstance(loop_node, NemoLoop) and loop_node.loop_type == 'lat' ):
+                    if isinstance(loop_node, NemoLoop) and loop_node.loop_type == 'lon':
+                        loop_node_interior = loop_node.loop_body
+                    loop_node = loop_node.parent
+                # Save a copy of the loop found
+                loop_copy = loop_node.copy()
+                # Identify parent of copy of test_node
+                npos = loop_node.walk(Node).index(test_node)
+                loop_copy_interior = loop_copy.walk(Node)[npos].parent
+                # Remove all children, and attach the nodes removed from the
+                # original IfBlock
+                loop_copy_interior.pop_all_children()
+                loop_copy_interior.addchild(IfBlock.create(test_copy.condition.copy(), warn_nodes, []))
+                # Tag loop interior as NemoKernel to trigger kernels region detection
+                # of the transformed loop Include initialisation of excessive-value
+                # counter
+                loop_node_interior.addchild(NemoKern(loop_node_interior.pop_all_children()))
+                # Add content of original test statement surrounded by the original
+                # loop in separate test statement outside of horizontal loop
+                counter0 = Assignment.create(Reference(test_sym), Literal('0', INTEGER_TYPE))
+                loop_node.parent.children.insert(loop_node.parent.children.index(loop_node), counter0)
+                new_node = IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.GE, Reference(test_sym), Literal('1', INTEGER_TYPE)), [loop_copy])
+                loop_node.parent.children.insert(loop_node.parent.children.index(loop_node)+1,new_node)
+            # Special case for invoke 'trc_bio_medusa' where the detection of
+            # excessive values *in the global domain* should trigger a call to
+            # subroutine 'trc_bio_check' that repeats the same check in the
+            # global domain and provides additional diagnostics
+            #
+            # Note, the flag is duplicated before the global reduction in order
+            # to keep the local version of the flag available
+            if flag_bio_check and 'trc_bio_medusa' in psy.invokes.names:
+                for ib in sched.walk(IfBlock):
+                    if ib.condition.walk(Reference)[0].name == flag_bio_check:
+                        ig = st.new_symbol('ifail_global',symbol_type=DataSymbol,datatype=ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE))
+                        new_ib = IfBlock.create(Reference(pst.lookup('lk_mpp')), [Call.create(pst.lookup('mpp_max'),[Literal('trcbio_medusa', CHARACTER_TYPE), Reference(ig)])])
+                        pos = ib.parent.children.index(ib)
+                        ib.parent.children.insert(pos, Assignment.create(Reference(ig),Reference(st.lookup(flag_bio_check))))
+                        ib.parent.children.insert(pos+1, new_ib)
+                        tbc = None
+                        for c in sched.walk(Call):
+                            if c.routine.name == 'trc_bio_check':
+                                tbc = c
+                        if tbc:
+                            tbc.detach()
+                            new_ib = IfBlock.create(BinaryOperation.create(BinaryOperation.Operator.GE, Reference(ig), Literal('1', INTEGER_TYPE)), [tbc])
+                            ib.parent.children.insert(pos+1, new_ib)
+
+        # ######## AC6 ########
+        # Special treatment of carbonate-chemistry component MOCSY
+        if match('^mocsy_.*', module_name):
+            # Replace options character strings by integer values
+            snames = list(st.symbols_dict.keys()) # static copy of symbol table contents
+            opt_varnames = []
+            for sname in snames:
+                if sname.startswith('op'):
+                    opt_sym = st.lookup(sname)
+                    if not isinstance(opt_sym.datatype, UnknownFortranType):
+                        opt_sym.datatype = ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE)
+                    else:
+                        repl = match('^CHARACTER[^,:]*(.*)$',opt_sym.datatype.declaration)
+                        if repl:
+                            opt_sym.datatype.declaration = 'INTEGER'+repl.group(1)
+                    opt_varnames += [sname]
+            for node in sched.walk((Assignment,CodeBlock)):
+                optval_node = None
+                # Assignments of option variables with option strings
+                if isinstance(node,Assignment) and node.lhs.name in opt_varnames and isinstance(node.rhs, Literal):
+                    optval_node = node.rhs
+                # Comparisons of option variables with option strings
+                elif isinstance(node,CodeBlock) and isinstance(node.ast, Fortran2003.Intrinsic_Function_Reference) and \
+                                                    str(node.ast.children[0]).lower() == 'trim' and \
+                                                    str(node.ast.children[1]).lower() in opt_varnames and \
+                                                    isinstance(node.parent,BinaryOperation) and \
+                                                    node.parent.operator == BinaryOperation.Operator.EQ:
+                    optval_node = node.parent.children[1]
+                    node.replace_with(Reference(st.lookup(str(node.ast.children[1]).lower())))
+                # Remove remaining 'trim' functions applied to option variables
+                elif isinstance(node,CodeBlock) and isinstance(node.ast, Fortran2003.Intrinsic_Function_Reference) and \
+                                                    str(node.ast.children[0]).lower() == 'trim' and \
+                                                    str(node.ast.children[1]).lower() in opt_varnames:
+                    node.replace_with(Reference(st.lookup(str(node.ast.children[1]).lower())))
+                # Catch-all for remaining 'trim' intrinsics in CodeBlocks:
+                # replace whole content of such CodeBlocks with 'STOP'
+                # statement
+                elif isinstance(node, CodeBlock):
+                    for an in node.get_ast_nodes:
+                        for ifr in walk(an,Fortran2003.Intrinsic_Function_Reference):
+                            if str(ifr.items[0]).lower() == 'trim':
+                                node.replace_with(CodeBlock([Fortran2003.Stop_Stmt('STOP')],CodeBlock.Structure.STATEMENT))
+                if optval_node:
+                    # Normalise name
+                    name = optval_node.value.replace('/','').lower()
+                    # Compute simple case-insensitive hash of option string
+                    # Note, only use first 5 characters; this is not fully safe,
+                    # the result fits into standard integer precision
+                    val = 0
+                    for nc in range(min(5,len(name.upper()))):
+                        val += (ord(name[nc])&31)*32**nc
+                    optval_node.replace_with(Literal(str(val),INTEGER_TYPE))
+            # Replace argument 'n' (if available) by integer parameter of value
+            # 1 in order to prevent the use of automatic arrays of shape
+            # 'dimension(n)' (if any)
+            # Note, in the unmodified source code, a value of 1 is passed to
+            # dummy argument 'n' of subroutine 'mocsy_interface' (which is used
+            # to set up the automatic arrays in this subroutine), but
+            # subsequently the a literal of 1 is used to set the value of 'n'
+            # throughout the MOCSY componenent; here we fix the size of the
+            # corresponding arrays to 1, and ignore the value of 'n' supplied
+            # to subroutine 'mocsy_interface'
+            if 'n' in st.symbols_dict.keys():
+                auto_arrays = []
+                for sym in st.symbols:
+                    if isinstance(sym,DataSymbol) and \
+                       len(sym.shape) == 1 and \
+                       not sym.is_argument and \
+                       isinstance(sym.shape[0][1],Reference) and \
+                       sym.shape[0][1].name == 'n':
+                            auto_arrays += [sym]
+                if len(auto_arrays) > 0:
+                    n_name = 'n_'+invoke.name+'_'+module_name
+                    pst.new_symbol(n_name,symbol_type=DataSymbol,datatype=ScalarType(ScalarType.Intrinsic.INTEGER,ScalarType.Precision.DOUBLE),constant_value=1)
+                    for auto_array in auto_arrays:
+                        auto_array.datatype = ArrayType(auto_array.datatype,shape=[ArrayType.ArrayBounds(lower=1,upper=Reference(pst.lookup(n_name)))])
+
+#############################################################################################
+# Section 2 - Exclusions of modules from generic transformations and kernels-region detection
+#############################################################################################
+
+    # No acceleration in specified modules
+    for module_re in ACC_EXCLUDE_MODULE_RE:
+        if match(module_re, module_name):
+            return psy
+
+    # Loop through all invokes
+    for invoke in psy.invokes.invoke_list:
+
+        # Skip exlucded invokes
+        if not invoke.name in ACC_EXCLUDE_INVOKES:
+
+            sched = invoke.schedule
+            st = sched.symbol_table
+            pst = st.parent_symbol_table()
+
+#################################################
+# Section 3 - Generic source-code transformations
+#################################################
+
+            # ######## AC10/AC13 #######
+            # Replacement of automatic arrays by module arrays in all modules that
+            # will be subject to kernels region detection
+            sym_auto = []
+            # Find automatic arrays: local, non-argument, non-constant arrays
+            # with defined array ranges
+            for sym in st.local_datasymbols:
+                deferred = False
+                if sym.is_local and not sym.is_argument and not sym.is_constant and sym.is_array:
+                    for dim in sym.shape:
+                        if not isinstance(dim, ArrayType.ArrayBounds):
+                            deferred = True
+                    if not deferred:
+                        sym_auto += [sym]
+            # Extract ranges used for automatic arrays, and create corresponding
+            # allocatable module arrays
+            sym_dict = {}
+            for sym in sym_auto:
+                rngs = []
+                for dim in sym.shape:
+                    rng = []
+                    for bnd in [dim.lower, dim.upper]:
+                        if isinstance(bnd, Literal):
+                            rng += [bnd.value]
+                        else:
+                            rng += [bnd.name]
+                    rngs += [':'.join(rng)]
+                sym_mod = pst.new_symbol(sym.name+'_'+invoke.name, symbol_type=DataSymbol, datatype=ArrayType(ScalarType(sym.datatype.intrinsic, sym.datatype.precision), shape=[ArrayType.Extent.DEFERRED]*len(rngs)))
+                sched.addchild(node_allocate(sym_mod, rngs), index=0)
+                sym_dict[sym.name] = sym_mod.name
+                # If precision of the variable is imported (e.g., 'use
+                # par_kind, only: wp') and if the import statement is present
+                # at procedure and not at module level, replicate the
+                # corresponding use association at the module scope (assume
+                # that if the precision symbol is imported at the procedure
+                # scope, it is not imported in a blanked module import and it
+                # is not imported at the module scope)
+                if isinstance(sym.datatype.precision, DataSymbol) and sym.datatype.precision.is_import and sym.datatype.precision.name in st.symbols_dict.keys() and not sym.datatype.precision.name in pst.symbols_dict.keys():
+                    psym = sym.datatype.precision
+                    csym = pst.new_symbol(psym.interface.container_symbol.name, symbol_type=ContainerSymbol)
+                    pst.new_symbol(psym.name, symbol_type=DataSymbol, datatype=psym.datatype, interface=ImportInterface(csym))
+            # Redirect all references (incl. references inside CodeBlock
+            # objects, possibly with multiple AST nodes) to automatic arrays to
+            # point to the newly created corresponding module arrays; further
+            # tests with source-code that is not currently processed revealed
+            # that the capitalisation of references in code blocks can differ
+            # form that used in the PSyIR, and therefore are noramlised here 
+            for ref in sched.walk((Reference,CodeBlock)):
+                if isinstance(ref, Reference) and ref.name in sym_dict.keys():
+                    ref.symbol = pst.lookup(sym_dict[ref.name])
+                elif isinstance(ref, CodeBlock):
+                    for ast in ref.get_ast_nodes:
+                        for name in walk(ast, Fortran2003.Name):
+                            if name.tostr().lower() in sym_dict.keys():
+                                name.string = sym_dict[name.tostr().lower()]
+            # Rename the now-unused automatic arrays (remove method for data
+            # symbols is not implemented)
+            for sym in sym_auto:
+                st.rename_symbol(sym, 'UNUSED___'+sym.name)
+
+##########################################
+# Section 4 - Detection of kernels regions
+##########################################
+
+            # List of all kernels and empty list of kernels taken care of
+            kernels = sched.kernels()
+            kernels_done = []
+
+            # Add to kernels list the list of all array assignments
+            # (implicit loops), but exclude those that contain CodeBlocks
+            for a in sched.walk(Assignment):
+                if a.is_array_range and not a.walk(CodeBlock):
+                    kernels += [a]
+
+
+            # Loop through kernels not yet taken care of
+            for kernel in kernels:
+                if not kernel in kernels_done:
+
+                    # Go up by one level if possible
+                    kernel_up = True
+                    while (kernel_up):
+                        kernel_up = up(kernel)
+                        if kernel_up:
+                            kernel = kernel_up
+
+                    kernel_region = [kernel]
+                    kernel_wide = kernel
+                    while (kernel_wide):
+                        kernel_wide = widen(kernel_wide)
+                        if kernel_wide:
+                            kernel_region += [kernel_wide]
+                    kernel_wide = kernel
+                    while (kernel_wide):
+                        kernel_wide = widen(kernel_wide, -1)
+                        if kernel_wide:
+                            kernel_region = [kernel_wide] + kernel_region
+
+                    # Add "!$acc kernels" directive around identified
+                    # group of nodes *unless* it is a solitary
+                    # InlinedKern (i.e. that hasn't "grown" and thus been
+                    # tested for unwanted content) node and mark all
+                    # kernels and assignments as taken care of
+                    if isinstance(kernel_region[0], (NemoLoop, IfBlock, Assignment)):
+                        for kernel in kernel_region:
+                            kernels_done += kernel.kernels()
+                            kernels_done += kernel.walk(Assignment)
+                            kernels_done += kernel.walk(Call)
+                        ACCKernelsTrans().apply(kernel_region)
+
+                    # Add "!$acc loop independent collapse(2)" directives
+                    # around 2-dimensional loops that span the horizontal
+                    # plane (i.e., nested loops of types 'lat' and 'lon',
+                    # respectively)
+                    for kernel in kernel_region:
+                        loop_lat = []
+                        loop_lon = []
+                        for loop in kernel.loops():
+                            # Detect 'lat' loops
+                            if (loop.loop_type == 'lat'):
+                                loop_lat += [ loop ]
+                            # Detect 'lon' loops
+                            elif (loop.loop_type == 'lon'):
+                                loop_lon += [ loop ]
+                        for loop in loop_lat:
+                            # Detect 'lat' loops which contain exactly one node that is a 'lon' loop
+                            if (len(loop.loop_body.children) == 1) and (loop.loop_body.children[0] in loop_lon):
+                                # Directives to fully parallelise and
+                                # collapse 2-dimensional loop
+                                ACCLoopTrans().apply(loop, {'independent': True, 'collapse': 2})
+
+    return psy
