Index: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zlys.F90
===================================================================
--- NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 11396)
+++ NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 13463)
@@ -20,5 +20,5 @@
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE p4zche          !  Chemical model
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
 
@@ -35,4 +35,6 @@
    REAL(wp) ::   calcon = 1.03E-2   ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
  
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -43,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_lys( kt, knt )
+   SUBROUTINE p4z_lys( kt, knt, Kbb, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_lys  ***
@@ -54,4 +56,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   ! ocean time step and ???
+      INTEGER, INTENT(in)  ::  Kbb, Krhs ! time level indices
       !
       INTEGER  ::   ji, jj, jk, jn
@@ -64,6 +67,4 @@
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
-      zco3    (:,:,:) = 0.
-      zcaldiss(:,:,:) = 0.
       zhinit  (:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
       !
@@ -72,15 +73,11 @@
       !     -------------------------------------------
 
-      CALL solve_at_general( zhinit, zhi )
+      CALL solve_at_general( zhinit, zhi, Kbb )
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
-                  &             + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
-               hi  (ji,jj,jk) = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000.
-            END DO
-         END DO
-      END DO
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+         zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
+            &             + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
+         hi  (ji,jj,jk) = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000.
+      END_3D
 
       !     ---------------------------------------------------------
@@ -90,47 +87,49 @@
       !     ---------------------------------------------------------
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
-               ! DEVIATION OF [CO3--] FROM SATURATION VALUE
-               ! Salinity dependance in zomegaca and divide by rhop/1000 to have good units
-               zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp )
-               zfact    = rhop(ji,jj,jk) / 1000._wp
-               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
-               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn )
+         ! DEVIATION OF [CO3--] FROM SATURATION VALUE
+         ! Salinity dependance in zomegaca and divide by rhop/1000 to have good units
+         zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp )
+         zfact    = rhop(ji,jj,jk) / 1000._wp
+         zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
+         zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn )
 
-               ! SET DEGREE OF UNDER-/SUPERSATURATION
-               excess(ji,jj,jk) = 1._wp - zomegaca
-               zexcess0 = MAX( 0., excess(ji,jj,jk) )
-               zexcess  = zexcess0**nca
+         ! SET DEGREE OF UNDER-/SUPERSATURATION
+         excess(ji,jj,jk) = 1._wp - zomegaca
+         zexcess0 = MAX( 0., excess(ji,jj,jk) )
+         zexcess  = zexcess0**nca
 
-               ! AMOUNT CACO3 (12C) THAT RE-ENTERS SOLUTION
-               !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
-               !       CACO3 GETS DISSOLVED EVEN IN THE CASE OF OVERSATURATION)
-               zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
-              !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
-              !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
-              zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
-              !
-              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
-              tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) -      zcaldiss(ji,jj,jk)
-              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +      zcaldiss(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+         ! AMOUNT CACO3 (12C) THAT RE-ENTERS SOLUTION
+         !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
+         !       CACO3 GETS DISSOLVED EVEN IN THE CASE OF OVERSATURATION)
+         zdispot = kdca * zexcess * tr(ji,jj,jk,jpcal,Kbb)
+        !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
+        !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
+        zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
+        !
+        tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * zcaldiss(ji,jj,jk)
+        tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) -      zcaldiss(ji,jj,jk)
+        tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) +      zcaldiss(ji,jj,jk)
+      END_3D
       !
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( MAX( hi(:,:,:), rtrn ) ) * tmask(:,:,:) )
-         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:)     * 1.e+3               * tmask(:,:,:) )
-         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:)  * 1.e+3               * tmask(:,:,:) )
-         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r     * tmask(:,:,:) )
+         CALL iom_put( "PH" , -1. * LOG10( MAX( hi(:,:,:), rtrn ) ) * tmask(:,:,:) )
+         IF( iom_use( "CO3" ) ) THEN
+            zco3(:,:,jpk) = 0.    ; CALL iom_put( "CO3"   , zco3(:,:,:)     * 1.e+3           * tmask(:,:,:) )
+         ENDIF
+         IF( iom_use( "CO3sat" ) ) THEN
+           zco3sat(:,:,jpk) = 0.  ; CALL iom_put( "CO3sat", zco3sat(:,:,:)  * 1.e+3           * tmask(:,:,:) )
+         ENDIF
+         IF( iom_use( "DCAL" ) ) THEN
+           zcaldiss(:,:,jpk) = 0. ; CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
+         ENDIF              
       ENDIF
       !
-      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="('lys ')")
-        CALL prt_ctl_trc_info(charout)
-        CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+        CALL prt_ctl_info( charout, cdcomp = 'top' )
+        CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -162,10 +161,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampiscal in reference namelist : Pisces CaCO3 dissolution
       READ  ( numnatp_ref, nampiscal, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampiscal in reference namelist', lwp )
-      REWIND( numnatp_cfg )              ! Namelist nampiscal in configuration namelist : Pisces CaCO3 dissolution
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampiscal in reference namelist' )
       READ  ( numnatp_cfg, nampiscal, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampiscal in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampiscal in configuration namelist' )
       IF(lwm) WRITE( numonp, nampiscal )
       !
