Index: NEMO/trunk/src/TOP/PISCES/P4Z/p4zflx.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 12277)
+++ NEMO/trunk/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 12377)
@@ -52,4 +52,6 @@
    REAL(wp) ::   xconv  = 0.01_wp / 3600._wp   !: coefficients for conversion 
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -59,5 +61,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_flx ( kt, knt )
+   SUBROUTINE p4z_flx ( kt, knt, Kbb, Kmm, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_flx  ***
@@ -71,4 +73,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   !
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs      ! time level indices
       !
       INTEGER  ::   ji, jj, jm, iind, iindm1
@@ -106,14 +109,12 @@
       IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
-            zfact = rhop(ji,jj,1) / 1000. + rtrn
-            zdic  = trb(ji,jj,1,jpdic)
-            zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
-            ! CALCULATE [H2CO3]
-            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
-         END DO
-      END DO
+      DO_2D_11_11
+         ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
+         zfact = rhop(ji,jj,1) / 1000. + rtrn
+         zdic  = tr(ji,jj,1,jpdic,Kbb)
+         zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
+         ! CALCULATE [H2CO3]
+         zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
+      END_2D
 
       ! --------------
@@ -124,50 +125,46 @@
       ! -------------------------------------------
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
-            ztc2 = ztc * ztc
-            ztc3 = ztc * ztc2 
-            ztc4 = ztc2 * ztc2 
-            ! Compute the schmidt Number both O2 and CO2
-            zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4
-            zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4
-            !  wind speed 
-            zws  = wndm(ji,jj) * wndm(ji,jj)
-            ! Compute the piston velocity for O2 and CO2
-            zkgwan = 0.251 * zws
-            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
-            ! compute gas exchange for CO2 and O2
-            zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
-            zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
-         END DO
-      END DO
-
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ztkel = tempis(ji,jj,1) + 273.15
-            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
-            zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
-            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
-            zxc2      = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2
-            zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
-            &           / ( 82.05736 * ztkel ))
-            zfco2 = zpco2atm(ji,jj) * zfugcoeff
-
-            ! Compute CO2 flux for the sea and air
-            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s)
-            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
-            oce_co2(ji,jj) = ( zfld - zflu ) * tmask(ji,jj,1) 
-            ! compute the trend
-            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + oce_co2(ji,jj) * rfact2 / e3t_n(ji,jj,1)
-
-            ! Compute O2 flux 
-            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
-            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj)
-            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
-            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / e3t_n(ji,jj,1)
-         END DO
-      END DO
+      DO_2D_11_11
+         ztc  = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) )
+         ztc2 = ztc * ztc
+         ztc3 = ztc * ztc2 
+         ztc4 = ztc2 * ztc2 
+         ! Compute the schmidt Number both O2 and CO2
+         zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4
+         zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4
+         !  wind speed 
+         zws  = wndm(ji,jj) * wndm(ji,jj)
+         ! Compute the piston velocity for O2 and CO2
+         zkgwan = 0.251 * zws
+         zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
+         ! compute gas exchange for CO2 and O2
+         zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
+         zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
+      END_2D
+
+
+      DO_2D_11_11
+         ztkel = tempis(ji,jj,1) + 273.15
+         zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
+         zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
+         zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
+         zxc2      = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2
+         zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
+         &           / ( 82.05736 * ztkel ))
+         zfco2 = zpco2atm(ji,jj) * zfugcoeff
+
+         ! Compute CO2 flux for the sea and air
+         zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s)
+         zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
+         oce_co2(ji,jj) = ( zfld - zflu ) * tmask(ji,jj,1) 
+         ! compute the trend
+         tr(ji,jj,1,jpdic,Krhs) = tr(ji,jj,1,jpdic,Krhs) + oce_co2(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm)
+
+         ! Compute O2 flux 
+         zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
+         zflu16 = tr(ji,jj,1,jpoxy,Kbb) * zkgo2(ji,jj)
+         zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
+         tr(ji,jj,1,jpoxy,Krhs) = tr(ji,jj,1,jpoxy,Krhs) + zoflx(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm)
+      END_2D
 
       IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst   &
@@ -178,8 +175,8 @@
       t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
  
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('flx ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
       ENDIF
 
@@ -191,5 +188,5 @@
          CALL iom_put( "Dpco2"   , ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
          CALL iom_put( "pCO2sea" , ( zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
-         CALL iom_put( "Dpo2"    , ( atcox * patm(:,:) - atcox * trb(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
+         CALL iom_put( "Dpo2"    , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
          CALL iom_put( "tcflx"   , t_oce_co2_flx     )   ! molC/s
          CALL iom_put( "tcflxcum", t_oce_co2_flx_cum )   ! molC
@@ -222,9 +219,6 @@
       ENDIF
       !
-      REWIND( numnatp_ref )
       READ  ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisext in reference namelist' )
-
-      REWIND( numnatp_cfg )
       READ  ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisext in configuration namelist' )
@@ -304,9 +298,6 @@
          ENDIF
          !
-         REWIND( numnatp_ref )
          READ  ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist' )
-
-         REWIND( numnatp_cfg )
          READ  ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 )
 902      IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisatm in configuration namelist' )
