Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 9169)
@@ -13,11 +13,12 @@
    USE oce_trc         !
    USE trc             ! 
-   USE sms_pisces
-   USE p2zopt
+   USE sms_pisces      !
+   USE p2zopt          !
+   USE trd_oce         !
+   USE trdtrc          !
+   !
    USE lbclnk          ! 
    USE prtctl_trc      ! Print control for debbuging
-   USE trd_oce
-   USE trdtrc
-   USE iom
+   USE iom             !
    
    IMPLICIT NONE
@@ -83,7 +84,6 @@
       !!        
       !!---------------------------------------------------------------------
-      !!
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
-      !!
+      !
       INTEGER  ::   ji, jj, jk, jl
       REAL(wp) ::   zdet, zzoo, zphy, zno3, znh4, zdom      ! now concentrations
@@ -96,6 +96,6 @@
       REAL(wp) ::   znh4a, zdeta, zdoma, zzoobod, zboddet, zdomaju
       REAL(wp) ::   ze3t
-      REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:) :: zw2d
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zw2d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zw3d
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -103,5 +103,5 @@
       IF( ln_timing )   CALL timing_start('p2z_bio')
       !
-      IF( lk_iomput ) ALLOCATE( zw2d(jpi,jpj,17), zw3d(jpi,jpj,jpk,3) )
+      IF( lk_iomput )   ALLOCATE( zw2d(jpi,jpj,17), zw3d(jpi,jpj,jpk,3) )
 
       IF( kt == nittrc000 ) THEN
@@ -113,6 +113,6 @@
       xksi(:,:) = 0.e0        ! zooplakton closure ( fbod)
       IF( lk_iomput ) THEN
-         zw2d  (:,:,:) = 0.e0
-         zw3d(:,:,:,:) = 0.e0
+         zw2d  (:,:,:) = 0._wp
+         zw3d(:,:,:,:) = 0._wp
       ENDIF
 
@@ -311,6 +311,5 @@
                tra(ji,jj,jk,jpdom) = tra(ji,jj,jk,jpdom) + zdoma
                !
-                IF( lk_iomput ) THEN
-                  ! convert fluxes in per day
+                IF( lk_iomput ) THEN                  ! convert fluxes in per day
                   ze3t = e3t_n(ji,jj,jk) * 86400._wp
                   zw2d(ji,jj,1)  = zw2d(ji,jj,1)  + zno3phy * ze3t
@@ -335,34 +334,34 @@
                   zw3d(ji,jj,jk,2) = znh4phy * 86400._wp
                   zw3d(ji,jj,jk,3) = znh4no3 * 86400._wp
-                   !
-                ENDIF
+                  !
+               ENDIF
             END DO
          END DO
       END DO
-
+      !
       IF( lk_iomput ) THEN
-        CALL lbc_lnk( zw2d(:,:,:),'T', 1. )
-        CALL lbc_lnk_multi( zw3d(:,:,:,1),'T', 1., zw3d(:,:,:,2),'T', 1., zw3d(:,:,:,3),'T', 1. )
-        ! Save diagnostics
-        CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
-        CALL iom_put( "TNH4PHY", zw2d(:,:,2) )
-        CALL iom_put( "TPHYDOM", zw2d(:,:,3) )
-        CALL iom_put( "TPHYNH4", zw2d(:,:,4) )
-        CALL iom_put( "TPHYZOO", zw2d(:,:,5) )
-        CALL iom_put( "TPHYDET", zw2d(:,:,6) )
-        CALL iom_put( "TDETZOO", zw2d(:,:,7) )
-        CALL iom_put( "TZOODET", zw2d(:,:,8) )
-        CALL iom_put( "TZOOBOD", zw2d(:,:,9) )
-        CALL iom_put( "TZOONH4", zw2d(:,:,10) )
-        CALL iom_put( "TZOODOM", zw2d(:,:,11) )
-        CALL iom_put( "TNH4NO3", zw2d(:,:,12) )
-        CALL iom_put( "TDOMNH4", zw2d(:,:,13) )
-        CALL iom_put( "TDETNH4", zw2d(:,:,14) )
-        CALL iom_put( "TPHYTOT", zw2d(:,:,15) )
-        CALL iom_put( "TZOOTOT", zw2d(:,:,16) )
+         CALL lbc_lnk( zw2d(:,:,:),'T', 1. )
+         CALL lbc_lnk_multi( zw3d(:,:,:,1),'T', 1., zw3d(:,:,:,2),'T', 1., zw3d(:,:,:,3),'T', 1. )
+         ! Save diagnostics
+         CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
+         CALL iom_put( "TNH4PHY", zw2d(:,:,2) )
+         CALL iom_put( "TPHYDOM", zw2d(:,:,3) )
+         CALL iom_put( "TPHYNH4", zw2d(:,:,4) )
+         CALL iom_put( "TPHYZOO", zw2d(:,:,5) )
+         CALL iom_put( "TPHYDET", zw2d(:,:,6) )
+         CALL iom_put( "TDETZOO", zw2d(:,:,7) )
+         CALL iom_put( "TZOODET", zw2d(:,:,8) )
+         CALL iom_put( "TZOOBOD", zw2d(:,:,9) )
+         CALL iom_put( "TZOONH4", zw2d(:,:,10) )
+         CALL iom_put( "TZOODOM", zw2d(:,:,11) )
+         CALL iom_put( "TNH4NO3", zw2d(:,:,12) )
+         CALL iom_put( "TDOMNH4", zw2d(:,:,13) )
+         CALL iom_put( "TDETNH4", zw2d(:,:,14) )
+         CALL iom_put( "TPHYTOT", zw2d(:,:,15) )
+         CALL iom_put( "TZOOTOT", zw2d(:,:,16) )
          ! 
-        CALL iom_put( "FNO3PHY", zw3d(:,:,:,1) )
-        CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
-        CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
+         CALL iom_put( "FNO3PHY", zw3d(:,:,:,1) )
+         CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
+         CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
          !
       ENDIF
@@ -374,5 +373,5 @@
       ENDIF
       !
-      IF( lk_iomput ) DEALLOCATE( zw2d, zw3d )
+      IF( lk_iomput )   DEALLOCATE( zw2d, zw3d )
       !
       IF( ln_timing )  CALL timing_stop('p2z_bio')
@@ -390,101 +389,100 @@
       !!
       !!----------------------------------------------------------------------
+      INTEGER ::   ios   ! Local integer
+      !!
       NAMELIST/namlobphy/ tmumax, rgamma, fphylab, tmminp, aki
       NAMELIST/namlobnut/ akno3, aknh4, taunn, psinut
       NAMELIST/namlobzoo/ rppz, taus, aks, rpnaz, rdnaz, tauzn, fzoolab, fdbod, tmminz
-      NAMELIST/namlobdet/  taudn, fdetlab
+      NAMELIST/namlobdet/ taudn, fdetlab
       NAMELIST/namlobdom/ taudomn
-      INTEGER :: ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
-
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' p2z_bio_init : LOBSTER bio-model initialization'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
+      !
       REWIND( numnatp_ref )              ! Namelist namlobphy in reference namelist : Lobster biological parameters
       READ  ( numnatp_ref, namlobphy, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobphy in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlobphy in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist namlobphy in configuration namelist : Lobster biological parameters
       READ  ( numnatp_cfg, namlobphy, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobphy in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlobphy in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, namlobphy )
-
+      !
       IF(lwp) THEN
-          WRITE(numout,*) ' Namelist namlobphy'
-          WRITE(numout,*) '    phyto max growth rate                                tmumax    =', 86400 * tmumax, ' d'
-          WRITE(numout,*) '    phytoplankton exudation fraction                     rgamma    =', rgamma
-          WRITE(numout,*) '    NH4 fraction of phytoplankton exsudation             fphylab   =', fphylab
-          WRITE(numout,*) '    minimal phyto mortality rate                         tmminp    =', 86400 * tmminp
-          WRITE(numout,*) '    light hlaf saturation constant                       aki       =', aki
-          WRITE(numout,*) ' '
+         WRITE(numout,*) '   Namelist namlobphy'
+         WRITE(numout,*) '      phyto max growth rate                                tmumax    =', 86400 * tmumax, ' d'
+         WRITE(numout,*) '      phytoplankton exudation fraction                     rgamma    =', rgamma
+         WRITE(numout,*) '      NH4 fraction of phytoplankton exsudation             fphylab   =', fphylab
+         WRITE(numout,*) '      minimal phyto mortality rate                         tmminp    =', 86400 * tmminp
+         WRITE(numout,*) '      light hlaf saturation constant                       aki       =', aki
       ENDIF
 
       REWIND( numnatp_ref )              ! Namelist namlobnut in reference namelist : Lobster nutriments parameters
       READ  ( numnatp_ref, namlobnut, IOSTAT = ios, ERR = 903)
-903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobnut in reference namelist', lwp )
-
+903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlobnut in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist namlobnut in configuration namelist : Lobster nutriments parameters
       READ  ( numnatp_cfg, namlobnut, IOSTAT = ios, ERR = 904 )
-904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobnut in configuration namelist', lwp )
+904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlobnut in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, namlobnut )
 
       IF(lwp) THEN
-          WRITE(numout,*) ' Namelist namlobnut'
-          WRITE(numout,*) '    half-saturation nutrient for no3 uptake              akno3     =', akno3
-          WRITE(numout,*) '    half-saturation nutrient for nh4 uptake              aknh4     =', aknh4
-          WRITE(numout,*) '    nitrification rate                                   taunn     =', taunn
-          WRITE(numout,*) '    inhibition of no3 uptake by nh4                      psinut    =', psinut
-          WRITE(numout,*) ' '
+         WRITE(numout,*) 
+         WRITE(numout,*) '   Namelist namlobnut'
+         WRITE(numout,*) '      half-saturation nutrient for no3 uptake              akno3     =', akno3
+         WRITE(numout,*) '      half-saturation nutrient for nh4 uptake              aknh4     =', aknh4
+         WRITE(numout,*) '      nitrification rate                                   taunn     =', taunn
+         WRITE(numout,*) '      inhibition of no3 uptake by nh4                      psinut    =', psinut
       ENDIF
 
       REWIND( numnatp_ref )              ! Namelist namlobzoo in reference namelist : Lobster zooplankton parameters
       READ  ( numnatp_ref, namlobzoo, IOSTAT = ios, ERR = 905)
-905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobzoo in reference namelist', lwp )
-
+905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlobzoo in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist namlobzoo in configuration namelist : Lobster zooplankton parameters
       READ  ( numnatp_cfg, namlobzoo, IOSTAT = ios, ERR = 906 )
-906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobzoo in configuration namelist', lwp )
+906   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlobzoo in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, namlobzoo )
 
       IF(lwp) THEN
-          WRITE(numout,*) ' Namelist namlobzoo'
-          WRITE(numout,*) '    zoo preference for phyto                             rppz      =', rppz
-          WRITE(numout,*) '    maximal zoo grazing rate                             taus      =', 86400 * taus, ' d'
-          WRITE(numout,*) '    half saturation constant for zoo food                aks       =', aks
-          WRITE(numout,*) '    non-assimilated phyto by zoo                         rpnaz     =', rpnaz
-          WRITE(numout,*) '    non-assimilated detritus by zoo                      rdnaz     =', rdnaz
-          WRITE(numout,*) '    zoo specific excretion rate                          tauzn     =', 86400 * tauzn
-          WRITE(numout,*) '    minimal zoo mortality rate                           tmminz    =', 86400 * tmminz
-          WRITE(numout,*) '    NH4 fraction of zooplankton excretion                fzoolab   =', fzoolab
-          WRITE(numout,*) '    Zooplankton mortality fraction that goes to detritus fdbod     =', fdbod
-          WRITE(numout,*) ' '
+         WRITE(numout,*) 
+         WRITE(numout,*) '   Namelist namlobzoo'
+         WRITE(numout,*) '      zoo preference for phyto                             rppz      =', rppz
+         WRITE(numout,*) '      maximal zoo grazing rate                             taus      =', 86400 * taus, ' d'
+         WRITE(numout,*) '      half saturation constant for zoo food                aks       =', aks
+         WRITE(numout,*) '      non-assimilated phyto by zoo                         rpnaz     =', rpnaz
+         WRITE(numout,*) '      non-assimilated detritus by zoo                      rdnaz     =', rdnaz
+         WRITE(numout,*) '      zoo specific excretion rate                          tauzn     =', 86400 * tauzn
+         WRITE(numout,*) '      minimal zoo mortality rate                           tmminz    =', 86400 * tmminz
+         WRITE(numout,*) '      NH4 fraction of zooplankton excretion                fzoolab   =', fzoolab
+         WRITE(numout,*) '      Zooplankton mortality fraction that goes to detritus fdbod     =', fdbod
       ENDIF
 
       REWIND( numnatp_ref )              ! Namelist namlobdet in reference namelist : Lobster detritus parameters
       READ  ( numnatp_ref, namlobdet, IOSTAT = ios, ERR = 907)
-907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobdet in reference namelist', lwp )
-
+907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlobdet in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist namlobdet in configuration namelist : Lobster detritus parameters
       READ  ( numnatp_cfg, namlobdet, IOSTAT = ios, ERR = 908 )
-908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobdet in configuration namelist', lwp )
+908   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlobdet in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, namlobdet )
 
       IF(lwp) THEN
-          WRITE(numout,*) ' Namelist namlobdet'
-          WRITE(numout,*) '    detrital breakdown rate                              taudn     =', 86400 * taudn , ' d'
-          WRITE(numout,*) '    NH4 fraction of detritus dissolution                 fdetlab   =', fdetlab
-          WRITE(numout,*) ' '
+          WRITE(numout,*) 
+          WRITE(numout,*) '   Namelist namlobdet'
+          WRITE(numout,*) '      detrital breakdown rate                              taudn     =', 86400 * taudn , ' d'
+          WRITE(numout,*) '      NH4 fraction of detritus dissolution                 fdetlab   =', fdetlab
       ENDIF
 
       REWIND( numnatp_ref )              ! Namelist namlobdom in reference namelist : Lobster DOM breakdown rate
       READ  ( numnatp_ref, namlobdom, IOSTAT = ios, ERR = 909)
-909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobdom in reference namelist', lwp )
-
+909   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlobdom in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist namlobdom in configuration namelist : Lobster DOM breakdown rate
       READ  ( numnatp_cfg, namlobdom, IOSTAT = ios, ERR = 910 )
-910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobdom in configuration namelist', lwp )
+910   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlobdom in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, namlobdom )
 
       IF(lwp) THEN
-          WRITE(numout,*) ' Namelist namlobdom'
-          WRITE(numout,*) '    DOM breakdown rate                                 taudomn     =', 86400 * taudn , ' d'
-          WRITE(numout,*) ' '
+          WRITE(numout,*) 
+          WRITE(numout,*) '   Namelist namlobdom'
+          WRITE(numout,*) '      DOM breakdown rate                                 taudomn     =', 86400 * taudn , ' d'
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90	(revision 9169)
@@ -11,10 +11,11 @@
    !!----------------------------------------------------------------------
    USE oce_trc         !
-   USE trc
-   USE sms_pisces
-   USE lbclnk
-   USE trd_oce
-   USE trdtrc
-   USE iom
+   USE trd_oce         !
+   USE trdtrc          !
+   USE trc             !
+   USE sms_pisces      !
+   !
+   USE lbclnk          !
+   USE iom             !
    USE prtctl_trc      ! Print control for debbuging
 
@@ -26,7 +27,7 @@
 
    REAL(wp), PUBLIC ::   sedlam      !: time coefficient of POC remineralization in sediments
-   REAL(wp), PUBLIC ::   sedlostpoc  ! mass of POC lost in sediments 
-   REAL(wp), PUBLIC ::   vsed        ! detritus sedimentation speed [m/s] 
-   REAL(wp), PUBLIC ::   xhr         ! coeff for martin''s remineralisation profile
+   REAL(wp), PUBLIC ::   sedlostpoc  !: mass of POC lost in sediments 
+   REAL(wp), PUBLIC ::   vsed        !: detritus sedimentation speed [m/s] 
+   REAL(wp), PUBLIC ::   xhr         !: coeff for martin''s remineralisation profile
 
    !!----------------------------------------------------------------------
@@ -126,22 +127,23 @@
       !!
       !!----------------------------------------------------------------------
+      INTEGER ::   ios   ! Local integer
+      !!
       NAMELIST/namlobsed/ sedlam, sedlostpoc, vsed, xhr
-      INTEGER :: ios                 ! Local integer output status for namelist read
-
+      !!----------------------------------------------------------------------
+      !
       REWIND( numnatp_ref )              ! Namelist namlobsed in reference namelist : Lobster sediments
       READ  ( numnatp_ref, namlobsed, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlosed in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlosed in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist namlobsed in configuration namelist : Lobster sediments
       READ  ( numnatp_cfg, namlobsed, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlobsed in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlobsed in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, namlobsed )
-
+      !
       IF(lwp) THEN
-          WRITE(numout,*) ' Namelist namlobsed'
-          WRITE(numout,*) '    time coeff of POC in sediments                       sedlam    =', sedlam
-          WRITE(numout,*) '    Sediment geol loss for POC                           sedlostpoc=', sedlostpoc
-          WRITE(numout,*) '    detritus sedimentation speed                         vsed      =', 86400 * vsed  , ' d'
-          WRITE(numout,*) '    coeff for martin''s remineralistion                  xhr       =', xhr
+          WRITE(numout,*) '   Namelist namlobsed'
+          WRITE(numout,*) '      time coeff of POC in sediments                sedlam    =', sedlam
+          WRITE(numout,*) '      Sediment geol loss for POC                    sedlostpoc=', sedlostpoc
+          WRITE(numout,*) '      detritus sedimentation speed                  vsed      =', 86400 * vsed  , ' d'
+          WRITE(numout,*) '      coeff for martin''s remineralistion           xhr       =', xhr
           WRITE(numout,*) ' '
       ENDIF
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 9169)
@@ -7,32 +7,31 @@
    !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-   !!   p4z_fechem       :  Compute remineralization/scavenging of iron
-   !!   p4z_fechem_init  :  Initialisation of parameters for remineralisation
-   !!   p4z_fechem_alloc :  Allocate remineralisation variables
+   !!   p4z_fechem       : Compute remineralization/scavenging of iron
+   !!   p4z_fechem_init  : Initialisation of parameters for remineralisation
+   !!   p4z_fechem_alloc : Allocate remineralisation variables
    !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zche          !  chemical model
-   USE p4zsbc          !  Boundary conditions from sediments
-   USE prtctl_trc      !  print control for debugging
-   USE iom             !  I/O manager
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE p4zche          ! chemical model
+   USE p4zsbc          ! Boundary conditions from sediments
+   USE prtctl_trc      ! print control for debugging
+   USE iom             ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   p4z_fechem      ! called in p4zbio.F90
-   PUBLIC   p4z_fechem_init ! called in trcsms_pisces.F90
-
-   !! * Shared module variables
-   LOGICAL          ::  ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
-   LOGICAL          ::  ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
-   LOGICAL          ::  ln_fecolloid !: boolean for variable colloidal fraction
-   REAL(wp), PUBLIC ::  xlam1        !: scavenging rate of Iron 
-   REAL(wp), PUBLIC ::  xlamdust     !: scavenging rate of Iron by dust 
-   REAL(wp), PUBLIC ::  ligand       !: ligand concentration in the ocean 
-   REAL(wp), PUBLIC ::  kfep         !: rate constant for nanoparticle formation
-
-   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth
+   PUBLIC   p4z_fechem        ! called in p4zbio.F90
+   PUBLIC   p4z_fechem_init   ! called in trcsms_pisces.F90
+
+   LOGICAL          ::   ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
+   LOGICAL          ::   ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
+   LOGICAL          ::   ln_fecolloid !: boolean for variable colloidal fraction
+   REAL(wp), PUBLIC ::   xlam1        !: scavenging rate of Iron 
+   REAL(wp), PUBLIC ::   xlamdust     !: scavenging rate of Iron by dust 
+   REAL(wp), PUBLIC ::   ligand       !: ligand concentration in the ocean 
+   REAL(wp), PUBLIC ::   kfep         !: rate constant for nanoparticle formation
+
+   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth      !!gm  <<<== DOCTOR names SVP !!!
 
    !!----------------------------------------------------------------------
@@ -56,6 +55,5 @@
       !!                    and one particulate form (ln_fechem)
       !!---------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
       !
       INTEGER  ::   ji, jj, jk, jic, jn
@@ -67,14 +65,14 @@
       REAL(wp) ::   zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2
       REAL(wp) ::   zrum, zcodel, zargu, zlight
-      REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
-      REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
-      REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
-      REAL(wp) :: ztfe, zoxy, zhplus
-      REAL(wp) :: zaggliga, zaggligb
-      REAL(wp) :: dissol, zligco
+      REAL(wp) ::   zkox, zkph1, zkph2, zph, zionic, ztligand
+      REAL(wp) ::   za, zb, zc, zkappa1, zkappa2, za0, za1, za2
+      REAL(wp) ::   zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
+      REAL(wp) ::   ztfe, zoxy, zhplus
+      REAL(wp) ::   zaggliga, zaggligb
+      REAL(wp) ::   dissol, zligco
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zTL1, zFe3, ztotlig, precip, zFeL1
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zFeL2, zTL2, zFe2, zFeP
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zstrn, zstrn2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, zFeL1
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zFeL2, zTL2, zFe2, zFeP
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::   zstrn, zstrn2
       !!---------------------------------------------------------------------
       !
@@ -384,30 +382,32 @@
       NAMELIST/nampisfer/ ln_fechem, ln_ligvar, ln_fecolloid, xlam1, xlamdust, ligand, kfep 
       !!----------------------------------------------------------------------
-
-      REWIND( numnatp_ref )              ! Namelist nampisfer in reference namelist : Pisces iron chemistry
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_rem_init : Initialization of iron chemistry parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
+      REWIND( numnatp_ref )            ! Namelist nampisfer in reference namelist : Pisces iron chemistry
       READ  ( numnatp_ref, nampisfer, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisfer in reference namelist', lwp )
-
-      REWIND( numnatp_cfg )              ! Namelist nampisfer in configuration namelist : Pisces iron chemistry
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisfer in reference namelist', lwp )
+      REWIND( numnatp_cfg )            ! Namelist nampisfer in configuration namelist : Pisces iron chemistry
       READ  ( numnatp_cfg, nampisfer, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampisfer in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampisfer )
-
-      IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem    =', ln_fechem
-         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar    =', ln_ligvar
-         WRITE(numout,*) '    Variable colloidal fraction of Fe3+       ln_fecolloid =', ln_fecolloid
-         WRITE(numout,*) '    scavenging rate of Iron                   xlam1        =', xlam1
-         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust     =', xlamdust
-         WRITE(numout,*) '    ligand concentration in the ocean         ligand       =', ligand
-         WRITE(numout,*) '    rate constant for nanoparticle formation  kfep         =', kfep
-      ENDIF
-      !
-      IF( ln_fechem ) THEN
-         ! initialization of some constants used by the complexe chemistry scheme
-         ! ----------------------------------------------------------------------
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisfer in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampisfer )
+
+      IF(lwp) THEN                     ! control print
+         WRITE(numout,*) '   Namelist : nampisfer'
+         WRITE(numout,*) '      enable complex iron chemistry scheme      ln_fechem    =', ln_fechem
+         WRITE(numout,*) '      variable concentration of ligand          ln_ligvar    =', ln_ligvar
+         WRITE(numout,*) '      Variable colloidal fraction of Fe3+       ln_fecolloid =', ln_fecolloid
+         WRITE(numout,*) '      scavenging rate of Iron                   xlam1        =', xlam1
+         WRITE(numout,*) '      scavenging rate of Iron by dust           xlamdust     =', xlamdust
+         WRITE(numout,*) '      ligand concentration in the ocean         ligand       =', ligand
+         WRITE(numout,*) '      rate constant for nanoparticle formation  kfep         =', kfep
+      ENDIF
+      !
+      IF( ln_fechem ) THEN             ! set some constants used by the complexe chemistry scheme
+         !
          spd = 3600. * 24.
          con = 1.E9
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 9169)
@@ -4,10 +4,10 @@
    !! TOP :   PISCES CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
    !!======================================================================
-   !! History :    -   !  1988-07  (E. MAIER-REIMER) Original code
-   !!              -   !  1998     (O. Aumont) additions
-   !!              -   !  1999     (C. Le Quere) modifications
-   !!             1.0  !  2004     (O. Aumont) modifications
-   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
-   !!                  !  2011-02  (J. Simeon, J. Orr) Include total atm P correction 
+   !! History :   -   !  1988-07  (E. MAIER-REIMER) Original code
+   !!             -   !  1998     (O. Aumont) additions
+   !!             -   !  1999     (C. Le Quere) modifications
+   !!            1.0  !  2004     (O. Aumont) modifications
+   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!                 !  2011-02  (J. Simeon, J. Orr) Include total atm P correction 
    !!----------------------------------------------------------------------
    !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
@@ -15,11 +15,11 @@
    !!   p4z_patm      :   Read sfc atm pressure [atm] for each grid cell
    !!----------------------------------------------------------------------
-   USE oce_trc                      !  shared variables between ocean and passive tracers 
-   USE trc                          !  passive tracers common variables
-   USE sms_pisces                   !  PISCES Source Minus Sink variables
-   USE p4zche                       !  Chemical model
-   USE prtctl_trc                   !  print control for debugging
-   USE iom                          !  I/O manager
-   USE fldread                      !  read input fields
+   USE oce_trc        !  shared variables between ocean and passive tracers 
+   USE trc            !  passive tracers common variables
+   USE sms_pisces     !  PISCES Source Minus Sink variables
+   USE p4zche         !  Chemical model
+   USE prtctl_trc     !  print control for debugging
+   USE iom            !  I/O manager
+   USE fldread        !  read input fields
 
    IMPLICIT NONE
@@ -30,25 +30,25 @@
    PUBLIC   p4z_flx_alloc  
 
-   !                               !!** Namelist  nampisext  **
-   REAL(wp)          ::  atcco2     !: pre-industrial atmospheric [co2] (ppm) 	
-   LOGICAL           ::  ln_co2int  !: flag to read in a file and interpolate atmospheric pco2 or not
-   CHARACTER(len=34) ::  clname     !: filename of pco2 values
-   INTEGER           ::  nn_offset  !: Offset model-data start year (default = 0) 
+   !                                 !!** Namelist  nampisext  **
+   REAL(wp)          ::   atcco2      !: pre-industrial atmospheric [co2] (ppm) 	
+   LOGICAL           ::   ln_co2int   !: flag to read in a file and interpolate atmospheric pco2 or not
+   CHARACTER(len=34) ::   clname      !: filename of pco2 values
+   INTEGER           ::   nn_offset   !: Offset model-data start year (default = 0) 
 
    !!  Variables related to reading atmospheric CO2 time history    
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: atcco2h, years
-   INTEGER  :: nmaxrec, numco2
-
-   !                               !!* nampisatm namelist (Atmospheric PRessure) *
+   INTEGER                                   ::   nmaxrec, numco2   !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   atcco2h, years    !
+
+   !                                  !!* nampisatm namelist (Atmospheric PRessure) *
    LOGICAL, PUBLIC ::   ln_presatm     !: ref. pressure: global mean Patm (F) or a constant (F)
    LOGICAL, PUBLIC ::   ln_presatmco2  !: accounting for spatial atm CO2 in the compuation of carbon flux (T) or not (F)
 
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) ::  patm      ! atmospheric pressure at kt                 [N/m2]
-   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::  sf_patm   ! structure of input fields (file informations, fields read)
-   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::  sf_atmco2 ! structure of input fields (file informations, fields read)
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: satmco2   !: atmospheric pco2 
-
-   REAL(wp) ::  xconv  = 0.01_wp / 3600._wp !: coefficients for conversion 
+   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) ::   patm      ! atmospheric pressure at kt                 [N/m2]
+   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::   sf_patm   ! structure of input fields (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::   sf_atmco2 ! structure of input fields (file informations, fields read)
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  satmco2   !: atmospheric pco2 
+
+   REAL(wp) ::   xconv  = 0.01_wp / 3600._wp   !: coefficients for conversion 
 
    !!----------------------------------------------------------------------
@@ -70,5 +70,4 @@
       !!              - Add option for time-interpolation of atcco2.txt  
       !!---------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt, knt   !
       !
@@ -79,17 +78,16 @@
       REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
       REAL(wp) ::   zyr_dec, zdco2dt
-      CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj) :: zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm  
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
+      CHARACTER (len=25) ::   charout
+      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm  
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_flx')
       !
-
       ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
       !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
       !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
 
-      IF( kt /= nit000 .AND. .NOT.l_co2cpl .AND. knt == 1 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
+      IF( kt /= nit000 .AND. .NOT.l_co2cpl .AND. knt == 1 )   CALL p4z_patm( kt )   ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
 
       IF( ln_co2int .AND. .NOT.ln_presatmco2 .AND. .NOT.l_co2cpl ) THEN 
@@ -226,40 +224,43 @@
       !! ** Method  :   Read the nampisext namelist and check the parameters
       !!      called at the first timestep (nittrc000)
+      !!
       !! ** input   :   Namelist nampisext
       !!----------------------------------------------------------------------
-      INTEGER ::   jm
-      INTEGER ::   ios   ! Local integer 
-      !
+      INTEGER ::   jm, ios   ! Local integer 
+      !!
       NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset
       !!----------------------------------------------------------------------
-      !
-
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' p4z_flx_init : atmospheric conditions for air-sea flux calculation'
+         WRITE(numout,*) ' ~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampisext in reference namelist : Pisces atm. conditions
       READ  ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisext in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisext in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisext in configuration namelist : Pisces atm. conditions
       READ  ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampisext in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisext in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, nampisext )
       !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for air-sea exchange, nampisext'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    Choice for reading in the atm pCO2 file or constant value, ln_co2int =', ln_co2int
-         WRITE(numout,*) ' '
-      ENDIF
+         WRITE(numout,*) '   Namelist : nampisext --- parameters for air-sea exchange'
+         WRITE(numout,*) '      reading in the atm pCO2 file or constant value   ln_co2int =', ln_co2int
+      ENDIF
+
+!!gm  BUG !!!   ===>>>  ln_presatm and ln_presatmco2 are used below, but read in namelist 
+!!gm                    at the end of the routine via a CALL to CALL p4z_patm( nit000 )
+
      IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
          IF(lwp) THEN                         ! control print
-            WRITE(numout,*) '    Constant Atmospheric pCO2 value  atcco2    =', atcco2
-            WRITE(numout,*) ' '
+            WRITE(numout,*) '         Constant Atmospheric pCO2 value               atcco2    =', atcco2
          ENDIF
          satmco2(:,:)  = atcco2      ! Initialisation of atmospheric pco2
       ELSEIF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
          IF(lwp)  THEN
-            WRITE(numout,*) '    Atmospheric pCO2 value  from file clname      =', TRIM( clname )
-            WRITE(numout,*) '    Offset model-data start year      nn_offset   =', nn_offset
-            WRITE(numout,*) ' '
+            WRITE(numout,*) '         Constant Atmospheric pCO2 value               atcco2    =', atcco2
+            WRITE(numout,*) '         Atmospheric pCO2 value  from file             clname    =', TRIM( clname )
+            WRITE(numout,*) '         Offset model-data start year                  nn_offset =', nn_offset
          ENDIF
          CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp )
@@ -270,7 +271,7 @@
          END DO
  100     nmaxrec = jm - 1 
-         ALLOCATE( years  (nmaxrec) )     ;      years  (:) = 0._wp
-         ALLOCATE( atcco2h(nmaxrec) )     ;      atcco2h(:) = 0._wp
-
+         ALLOCATE( years  (nmaxrec) )   ;   years  (:) = 0._wp
+         ALLOCATE( atcco2h(nmaxrec) )   ;   atcco2h(:) = 0._wp
+         !
          REWIND(numco2)
          DO jm = 1, nmaxrec          ! get  xCO2 data
@@ -282,10 +283,8 @@
          IF(lwp)  THEN
             WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
-            WRITE(numout,*) ' '
          ENDIF
       ELSE
          IF(lwp)  THEN
             WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
-            WRITE(numout,*) ' '
          ENDIF
       ENDIF
@@ -304,38 +303,39 @@
       !!                  ***  ROUTINE p4z_atm  ***
       !!
-      !! ** Purpose :   Read and interpolate the external atmospheric sea-levl pressure
+      !! ** Purpose :   Read and interpolate the external atmospheric sea-level pressure
       !! ** Method  :   Read the files and interpolate the appropriate variables
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
-      !
-      INTEGER            ::  ierr
-      INTEGER            ::  ios      ! Local integer output status for namelist read
-      CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files
-      TYPE(FLD_N)        ::  sn_patm  ! informations about the fields to be read
-      TYPE(FLD_N)        ::  sn_atmco2 ! informations about the fields to be read
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !
+      INTEGER            ::   ierr, ios   ! Local integer
+      CHARACTER(len=100) ::   cn_dir      ! Root directory for location of ssr files
+      TYPE(FLD_N)        ::   sn_patm     ! informations about the fields to be read
+      TYPE(FLD_N)        ::   sn_atmco2   ! informations about the fields to be read
       !!
       NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
       !!----------------------------------------------------------------------
-
-      !                                         ! ----------------------- !
-      IF( kt == nit000 ) THEN                   ! First call kt=nittrc000 !
-
+      !
+      IF( kt == nit000 ) THEN    !==  First call kt=nittrc000  ==!
+         !
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) ' p4z_patm : sea-level atmospheric pressure'
+            WRITE(numout,*) ' ~~~~~~~~'
+         ENDIF
+         !
          REWIND( numnatp_ref )              ! Namelist nampisatm in reference namelist : Pisces atm. sea level pressure file
          READ  ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist', lwp )
-
          REWIND( numnatp_cfg )              ! Namelist nampisatm in configuration namelist : Pisces atm. sea level pressure file 
          READ  ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in configuration namelist', lwp )
+902      IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisatm in configuration namelist', lwp )
          IF(lwm) WRITE ( numonp, nampisatm )
          !
          !
          IF(lwp) THEN                                 !* control print
-            WRITE(numout,*)
-            WRITE(numout,*) '   Namelist nampisatm : Atmospheric Pressure as external forcing'
-            WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm = ', ln_presatm
-            WRITE(numout,*) '      spatial atmopsheric CO2 for flux calcs  ln_presatmco2 = ', ln_presatmco2
-            WRITE(numout,*)
+            WRITE(numout,*) '   Namelist : nampisatm --- Atmospheric Pressure as external forcing'
+            WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm    = ', ln_presatm
+            WRITE(numout,*) '      spatial atmopsheric CO2 for flux calcs                ln_presatmco2 = ', ln_presatmco2
          ENDIF
          !
@@ -358,5 +358,5 @@
          ENDIF
          !
-         IF( .NOT.ln_presatm )   patm(:,:) = 1.e0    ! Initialize patm if no reading from a file
+         IF( .NOT.ln_presatm )   patm(:,:) = 1._wp    ! Initialize patm if no reading from a file
          !
       ENDIF
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 9169)
@@ -6,11 +6,11 @@
    !! History :   3.6  !  2016-03  (O. Aumont, A. Tagliabue) Quota model and reorganization
    !!----------------------------------------------------------------------
-   !!   p4z_ligand       :  Compute remineralization/dissolution of organic ligands
-   !!   p4z_ligand_init  :  Initialisation of parameters for remineralisation
+   !!   p4z_ligand     :  Compute remineralization/dissolution of organic ligands
+   !!   p4z_ligand_init:  Initialisation of parameters for remineralisation
    !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE prtctl_trc      !  print control for debugging
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE prtctl_trc      ! print control for debugging
 
    IMPLICIT NONE
@@ -20,5 +20,4 @@
    PUBLIC   p4z_ligand_init    ! called in trcsms_pisces.F90
 
-   !! * Shared module variables
    REAL(wp), PUBLIC ::  rlgw     !: lifetime (years) of weak ligands
    REAL(wp), PUBLIC ::  rlgs     !: lifetime (years) of strong ligands
@@ -39,23 +38,21 @@
       !!
       !! ** Purpose :   Compute remineralization/scavenging of organic ligands
-      !!
-      !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       !
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zlgwp, zlgwpr, zlgwr, zlablgw, zrfepa, zfepr
-      CHARACTER (len=25) :: charout
+      CHARACTER (len=25) ::   charout
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_ligand')
       !
-      ! ------------------------------------------------------------------
-      ! Remineralization of iron ligands
-      ! ------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
+               !
+               ! ------------------------------------------------------------------
+               ! Remineralization of iron ligands
+               ! ------------------------------------------------------------------
                ! production from remineralisation of organic matter
                zlgwp  = orem(ji,jj,jk) * rlig
@@ -68,14 +65,8 @@
                zlgwpr = prlgw * xstep * etot(ji,jj,jk) * trb(ji,jj,jk,jplgw) * (1. - fr_i(ji,jj))
                tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zlgwp - zlgwr - zlgwpr
-            END DO
-         END DO
-      END DO
-
-      ! ----------------------------------------------------------
-      ! Dissolution of nanoparticle Fe
-      ! ----------------------------------------------------------
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+               !
+               ! ----------------------------------------------------------
+               ! Dissolution of nanoparticle Fe
+               ! ----------------------------------------------------------
                ! dissolution rate is maximal in the presence of light and 
                ! lower in the aphotici zone
@@ -86,13 +77,14 @@
                tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) - zfepr
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zfepr
+               !
             END DO
          END DO
       END DO
-
+      !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('ligand1')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-       ENDIF
+      ENDIF
       !
       IF( ln_timing )   CALL timing_stop('p4z_ligand')
@@ -108,8 +100,6 @@
       !!
       !! ** Method  :   Read the nampislig namelist and check the parameters
-      !!      called at the first timestep
       !!
       !! ** input   :   Namelist nampislig
-      !!
       !!----------------------------------------------------------------------
       INTEGER ::   ios   ! Local integer 
@@ -117,23 +107,25 @@
       NAMELIST/nampislig/ rlgw, prlgw, rlgs, rfep, rlig
       !!----------------------------------------------------------------------
-
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_ligand_init : remineralization/scavenging of organic ligands'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+      ENDIF
       REWIND( numnatp_ref )              ! Namelist nampislig in reference namelist : Pisces remineralization
       READ  ( numnatp_ref, nampislig, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampislig in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampislig in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampislig in configuration namelist : Pisces remineralization
       READ  ( numnatp_cfg, nampislig, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampislig in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampislig in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, nampislig )
-
+      !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for ligands, nampislig'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    Dissolution rate of FeP                        rfep =', rfep
-         WRITE(numout,*) '    Lifetime (years) of weak ligands               rlgw =', rlgw
-         WRITE(numout,*) '    Remin ligand production per unit C             rlig =', rlig
-         WRITE(numout,*) '    Photolysis of weak ligand                     prlgw =', prlgw
-         WRITE(numout,*) '    Lifetime (years) of strong ligands             rlgs =', rlgs
+         WRITE(numout,*) '   Namelist : nampislig'
+         WRITE(numout,*) '      Dissolution rate of FeP                      rfep  =', rfep
+         WRITE(numout,*) '      Lifetime (years) of weak ligands             rlgw  =', rlgw
+         WRITE(numout,*) '      Remin ligand production per unit C           rlig  =', rlig
+         WRITE(numout,*) '      Photolysis of weak ligand                    prlgw =', prlgw
+         WRITE(numout,*) '      Lifetime (years) of strong ligands           rlgs  =', rlgs
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 9169)
@@ -215,9 +215,9 @@
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
-        IF( iom_use( "xfracal" ) ) CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
-        IF( iom_use( "LNnut"   ) ) CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        IF( iom_use( "LDnut"   ) ) CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        IF( iom_use( "LNFe"    ) ) CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        IF( iom_use( "LDFe"    ) ) CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "xfracal" ) )   CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
+        IF( iom_use( "LNnut"   ) )   CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LDnut"   ) )   CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LNFe"    ) )   CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "LDFe"    ) )   CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
       ENDIF
       !
@@ -246,37 +246,40 @@
       !!----------------------------------------------------------------------
       !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_lim_init : initialization of nutrient limitations'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampislim in reference namelist : Pisces nutrient limitation parameters
       READ  ( numnatp_ref, namp4zlim, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp4zlim in reference namelist', lwp )
-      !
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zlim in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampislim in configuration namelist : Pisces nutrient limitation parameters 
       READ  ( numnatp_cfg, namp4zlim, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp4zlim in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, namp4zlim )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zlim in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, namp4zlim )
       !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for nutrient limitations, namp4zlim'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    mean rainratio                           caco3r    = ', caco3r
-         WRITE(numout,*) '    NO3 half saturation of nanophyto         concnno3  = ', concnno3
-         WRITE(numout,*) '    NO3 half saturation of diatoms           concdno3  = ', concdno3
-         WRITE(numout,*) '    NH4 half saturation for phyto            concnnh4  = ', concnnh4
-         WRITE(numout,*) '    NH4 half saturation for diatoms          concdnh4  = ', concdnh4
-         WRITE(numout,*) '    half saturation constant for Si uptake   xksi1     = ', xksi1
-         WRITE(numout,*) '    half saturation constant for Si/C        xksi2     = ', xksi2
-         WRITE(numout,*) '    half-sat. of DOC remineralization        xkdoc     = ', xkdoc
-         WRITE(numout,*) '    Iron half saturation for nanophyto       concnfer  = ', concnfer
-         WRITE(numout,*) '    Iron half saturation for diatoms         concdfer  = ', concdfer
-         WRITE(numout,*) '    size ratio for nanophytoplankton         xsizern   = ', xsizern
-         WRITE(numout,*) '    size ratio for diatoms                   xsizerd   = ', xsizerd
-         WRITE(numout,*) '    NO3 half saturation of bacteria          concbno3  = ', concbno3
-         WRITE(numout,*) '    NH4 half saturation for bacteria         concbnh4  = ', concbnh4
-         WRITE(numout,*) '    Minimum size criteria for diatoms        xsizedia  = ', xsizedia
-         WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy
-         WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe
-         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin
-         WRITE(numout,*) '    optimal Fe quota for nano.               qnfelim   = ', qnfelim
-         WRITE(numout,*) '    Optimal Fe quota for diatoms             qdfelim   = ', qdfelim
+         WRITE(numout,*) '   Namelist : namp4zlim'
+         WRITE(numout,*) '      mean rainratio                           caco3r    = ', caco3r
+         WRITE(numout,*) '      NO3 half saturation of nanophyto         concnno3  = ', concnno3
+         WRITE(numout,*) '      NO3 half saturation of diatoms           concdno3  = ', concdno3
+         WRITE(numout,*) '      NH4 half saturation for phyto            concnnh4  = ', concnnh4
+         WRITE(numout,*) '      NH4 half saturation for diatoms          concdnh4  = ', concdnh4
+         WRITE(numout,*) '      half saturation constant for Si uptake   xksi1     = ', xksi1
+         WRITE(numout,*) '      half saturation constant for Si/C        xksi2     = ', xksi2
+         WRITE(numout,*) '      half-sat. of DOC remineralization        xkdoc     = ', xkdoc
+         WRITE(numout,*) '      Iron half saturation for nanophyto       concnfer  = ', concnfer
+         WRITE(numout,*) '      Iron half saturation for diatoms         concdfer  = ', concdfer
+         WRITE(numout,*) '      size ratio for nanophytoplankton         xsizern   = ', xsizern
+         WRITE(numout,*) '      size ratio for diatoms                   xsizerd   = ', xsizerd
+         WRITE(numout,*) '      NO3 half saturation of bacteria          concbno3  = ', concbno3
+         WRITE(numout,*) '      NH4 half saturation for bacteria         concbnh4  = ', concbnh4
+         WRITE(numout,*) '      Minimum size criteria for diatoms        xsizedia  = ', xsizedia
+         WRITE(numout,*) '      Minimum size criteria for nanophyto      xsizephy  = ', xsizephy
+         WRITE(numout,*) '      Fe half saturation for bacteria          concbfe   = ', concbfe
+         WRITE(numout,*) '      halk saturation constant for anoxia       oxymin   =' , oxymin
+         WRITE(numout,*) '      optimal Fe quota for nano.               qnfelim   = ', qnfelim
+         WRITE(numout,*) '      Optimal Fe quota for diatoms             qdfelim   = ', qdfelim
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 9169)
@@ -29,13 +29,10 @@
    PUBLIC   p4z_lys_init    ! called in trcsms_pisces.F90
 
-   !! * Shared module variables
-   REAL(wp), PUBLIC :: kdca !: diss. rate constant calcite
-   REAL(wp), PUBLIC :: nca  !: order of reaction for calcite dissolution
+   REAL(wp), PUBLIC ::   kdca   !: diss. rate constant calcite
+   REAL(wp), PUBLIC ::   nca    !: order of reaction for calcite dissolution
 
-   !! * Module variables
-   REAL(wp) :: calcon = 1.03E-2           !: mean calcite concentration [Ca2+] in sea water [mole/kg solution]
+   INTEGER  ::   rmtss              ! number of seconds per month 
+   REAL(wp) ::   calcon = 1.03E-2   ! mean calcite concentration [Ca2+] in sea water [mole/kg solution]
  
-   INTEGER  :: rmtss                      !: number of seconds per month 
-
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -56,11 +53,11 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step and ???
       !
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zdispot, zfact, zcalcon
       REAL(wp) ::   zomegaca, zexcess, zexcess0
-      CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zco3, zcaldiss, zhinit, zhi, zco3sat
+      CHARACTER (len=25) ::   charout
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zco3, zcaldiss, zhinit, zhi, zco3sat
       !!---------------------------------------------------------------------
       !
@@ -69,10 +66,11 @@
       zco3    (:,:,:) = 0.
       zcaldiss(:,:,:) = 0.
-      zhinit(:,:,:)   = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
+      zhinit  (:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
+      !
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
       !     -------------------------------------------
 
-      CALL solve_at_general(zhinit, zhi)
+      CALL solve_at_general( zhinit, zhi )
 
       DO jk = 1, jpkm1
@@ -80,6 +78,6 @@
             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.
+                  &             + 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
@@ -158,20 +156,22 @@
       NAMELIST/nampiscal/ kdca, nca
       !!----------------------------------------------------------------------
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_lys_init : initialization of CaCO3 dissolution'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      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 )
-      !
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampiscal in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampiscal in configuration namelist : Pisces CaCO3 dissolution
       READ  ( numnatp_cfg, nampiscal, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampiscal in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampiscal )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampiscal in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampiscal )
       !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for CaCO3 dissolution, nampiscal'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    diss. rate constant calcite (per month)   kdca      =', kdca
-         WRITE(numout,*) '    order of reaction for calcite dissolution nca       =', nca
+         WRITE(numout,*) '   Namelist : nampiscal'
+         WRITE(numout,*) '      diss. rate constant calcite (per month)        kdca =', kdca
+         WRITE(numout,*) '      order of reaction for calcite dissolution      nca  =', nca
       ENDIF
       !
@@ -180,4 +180,5 @@
       !
    END SUBROUTINE p4z_lys_init
+
    !!======================================================================
 END MODULE p4zlys
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 9169)
@@ -8,13 +8,13 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Quota model for iron
    !!----------------------------------------------------------------------
-   !!   p4z_meso       :   Compute the sources/sinks for mesozooplankton
-   !!   p4z_meso_init  :   Initialization of the parameters for mesozooplankton
-   !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zprod         !  production
-   USE prtctl_trc      !  print control for debugging
-   USE iom             !  I/O manager
+   !!   p4z_meso       : Compute the sources/sinks for mesozooplankton
+   !!   p4z_meso_init  : Initialization of the parameters for mesozooplankton
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE p4zprod         ! production
+   USE prtctl_trc      ! print control for debugging
+   USE iom             ! I/O manager
 
    IMPLICIT NONE
@@ -24,5 +24,4 @@
    PUBLIC   p4z_meso_init         ! called in trcsms_pisces.F90
 
-   !! * Shared module variables
    REAL(wp), PUBLIC ::  part2        !: part of calcite not dissolved in mesozoo guts
    REAL(wp), PUBLIC ::  xprefc       !: mesozoo preference for POC 
@@ -49,5 +48,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -60,5 +58,6 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step and ???
+      !
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz, zcompam
@@ -73,6 +72,5 @@
       CHARACTER (len=25) :: charout
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
-
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
       !!---------------------------------------------------------------------
       !
@@ -122,5 +120,4 @@
 
                !  Mesozooplankton flux feeding on GOC
-               !  ----------------------------------
                !  ----------------------------------
                zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
@@ -253,5 +250,4 @@
       !!
       !! ** input   :   Namelist nampismes
-      !!
       !!----------------------------------------------------------------------
       INTEGER ::   ios   ! Local integer
@@ -262,35 +258,38 @@
       !!----------------------------------------------------------------------
       !
+      IF(lwp) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) 'p4z_meso_init : Initialization of mesozooplankton parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampismes in reference namelist : Pisces mesozooplankton
       READ  ( numnatp_ref, namp4zmes, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp4zmes in reference namelist', lwp )
-      !
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zmes in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampismes in configuration namelist : Pisces mesozooplankton
       READ  ( numnatp_cfg, namp4zmes, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp4zmes in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, namp4zmes )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zmes in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, namp4zmes )
       !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' ' 
-         WRITE(numout,*) ' Namelist parameters for mesozooplankton, namp4zmes'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    part of calcite not dissolved in mesozoo guts  part2        =', part2
-         WRITE(numout,*) '    mesozoo preference for phyto                   xprefc       =', xprefc
-         WRITE(numout,*) '    mesozoo preference for POC                     xprefp       =', xprefp
-         WRITE(numout,*) '    mesozoo preference for zoo                     xprefz       =', xprefz
-         WRITE(numout,*) '    mesozoo preference for poc                     xprefpoc     =', xprefpoc
-         WRITE(numout,*) '    microzoo feeding threshold  for mesozoo        xthresh2zoo  =', xthresh2zoo
-         WRITE(numout,*) '    diatoms feeding threshold  for mesozoo         xthresh2dia  =', xthresh2dia
-         WRITE(numout,*) '    nanophyto feeding threshold for mesozoo        xthresh2phy  =', xthresh2phy
-         WRITE(numout,*) '    poc feeding threshold for mesozoo              xthresh2poc  =', xthresh2poc
-         WRITE(numout,*) '    feeding threshold for mesozooplankton          xthresh2     =', xthresh2
-         WRITE(numout,*) '    exsudation rate of mesozooplankton             resrat2      =', resrat2
-         WRITE(numout,*) '    mesozooplankton mortality rate                 mzrat2       =', mzrat2
-         WRITE(numout,*) '    maximal mesozoo grazing rate                   grazrat2     =', grazrat2
-         WRITE(numout,*) '    mesozoo flux feeding rate                      grazflux     =', grazflux
-         WRITE(numout,*) '    non assimilated fraction of P by mesozoo       unass2       =', unass2
-         WRITE(numout,*) '    Efficicency of Mesozoo growth                  epsher2      =', epsher2
-         WRITE(numout,*) '    Fraction of mesozoo excretion as DOM           sigma2       =', sigma2
-         WRITE(numout,*) '    half sturation constant for grazing 2          xkgraz2      =', xkgraz2
+         WRITE(numout,*) '   Namelist : namp4zmes'
+         WRITE(numout,*) '      part of calcite not dissolved in mesozoo guts  part2        =', part2
+         WRITE(numout,*) '      mesozoo preference for phyto                   xprefc       =', xprefc
+         WRITE(numout,*) '      mesozoo preference for POC                     xprefp       =', xprefp
+         WRITE(numout,*) '      mesozoo preference for zoo                     xprefz       =', xprefz
+         WRITE(numout,*) '      mesozoo preference for poc                     xprefpoc     =', xprefpoc
+         WRITE(numout,*) '      microzoo feeding threshold  for mesozoo        xthresh2zoo  =', xthresh2zoo
+         WRITE(numout,*) '      diatoms feeding threshold  for mesozoo         xthresh2dia  =', xthresh2dia
+         WRITE(numout,*) '      nanophyto feeding threshold for mesozoo        xthresh2phy  =', xthresh2phy
+         WRITE(numout,*) '      poc feeding threshold for mesozoo              xthresh2poc  =', xthresh2poc
+         WRITE(numout,*) '      feeding threshold for mesozooplankton          xthresh2     =', xthresh2
+         WRITE(numout,*) '      exsudation rate of mesozooplankton             resrat2      =', resrat2
+         WRITE(numout,*) '      mesozooplankton mortality rate                 mzrat2       =', mzrat2
+         WRITE(numout,*) '      maximal mesozoo grazing rate                   grazrat2     =', grazrat2
+         WRITE(numout,*) '      mesozoo flux feeding rate                      grazflux     =', grazflux
+         WRITE(numout,*) '      non assimilated fraction of P by mesozoo       unass2       =', unass2
+         WRITE(numout,*) '      Efficicency of Mesozoo growth                  epsher2      =', epsher2
+         WRITE(numout,*) '      Fraction of mesozoo excretion as DOM           sigma2       =', sigma2
+         WRITE(numout,*) '      half sturation constant for grazing 2          xkgraz2      =', xkgraz2
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 9169)
@@ -8,14 +8,14 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Quota model for iron
    !!----------------------------------------------------------------------
-   !!   p4z_micro       :   Compute the sources/sinks for microzooplankton
-   !!   p4z_micro_init  :   Initialize and read the appropriate namelist
-   !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zlim          !  Co-limitations
-   USE p4zprod         !  production
-   USE iom             !  I/O manager
-   USE prtctl_trc      !  print control for debugging
+   !!   p4z_micro      : Compute the sources/sinks for microzooplankton
+   !!   p4z_micro_init : Initialize and read the appropriate namelist
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE p4zlim          ! Co-limitations
+   USE p4zprod         ! production
+   USE iom             ! I/O manager
+   USE prtctl_trc      ! print control for debugging
 
    IMPLICIT NONE
@@ -25,20 +25,19 @@
    PUBLIC   p4z_micro_init    ! called in trcsms_pisces.F90
 
-   !! * Shared module variables
-   REAL(wp), PUBLIC ::  part        !: part of calcite not dissolved in microzoo guts
-   REAL(wp), PUBLIC ::  xpref2c     !: microzoo preference for POC 
-   REAL(wp), PUBLIC ::  xpref2p     !: microzoo preference for nanophyto
-   REAL(wp), PUBLIC ::  xpref2d     !: microzoo preference for diatoms
-   REAL(wp), PUBLIC ::  xthreshdia  !: diatoms feeding threshold for microzooplankton 
-   REAL(wp), PUBLIC ::  xthreshphy  !: nanophyto threshold for microzooplankton 
-   REAL(wp), PUBLIC ::  xthreshpoc  !: poc threshold for microzooplankton 
-   REAL(wp), PUBLIC ::  xthresh     !: feeding threshold for microzooplankton 
-   REAL(wp), PUBLIC ::  resrat      !: exsudation rate of microzooplankton
-   REAL(wp), PUBLIC ::  mzrat       !: microzooplankton mortality rate 
-   REAL(wp), PUBLIC ::  grazrat     !: maximal microzoo grazing rate
-   REAL(wp), PUBLIC ::  xkgraz      !: non assimilated fraction of P by microzoo 
-   REAL(wp), PUBLIC ::  unass       !: Efficicency of microzoo growth 
-   REAL(wp), PUBLIC ::  sigma1      !: Fraction of microzoo excretion as DOM 
-   REAL(wp), PUBLIC ::  epsher      !: half sturation constant for grazing 1 
+   REAL(wp), PUBLIC ::   part        !: part of calcite not dissolved in microzoo guts
+   REAL(wp), PUBLIC ::   xpref2c     !: microzoo preference for POC 
+   REAL(wp), PUBLIC ::   xpref2p     !: microzoo preference for nanophyto
+   REAL(wp), PUBLIC ::   xpref2d     !: microzoo preference for diatoms
+   REAL(wp), PUBLIC ::   xthreshdia  !: diatoms feeding threshold for microzooplankton 
+   REAL(wp), PUBLIC ::   xthreshphy  !: nanophyto threshold for microzooplankton 
+   REAL(wp), PUBLIC ::   xthreshpoc  !: poc threshold for microzooplankton 
+   REAL(wp), PUBLIC ::   xthresh     !: feeding threshold for microzooplankton 
+   REAL(wp), PUBLIC ::   resrat      !: exsudation rate of microzooplankton
+   REAL(wp), PUBLIC ::   mzrat       !: microzooplankton mortality rate 
+   REAL(wp), PUBLIC ::   grazrat     !: maximal microzoo grazing rate
+   REAL(wp), PUBLIC ::   xkgraz      !: non assimilated fraction of P by microzoo 
+   REAL(wp), PUBLIC ::   unass       !: Efficicency of microzoo growth 
+   REAL(wp), PUBLIC ::   sigma1      !: Fraction of microzoo excretion as DOM 
+   REAL(wp), PUBLIC ::   epsher      !: half sturation constant for grazing 1 
 
    !!----------------------------------------------------------------------
@@ -57,6 +56,6 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::  kt  ! ocean time step
-      INTEGER, INTENT(in) ::  knt 
+      INTEGER, INTENT(in) ::   kt    ! ocean time step
+      INTEGER, INTENT(in) ::   knt   ! ??? 
       !
       INTEGER  :: ji, jj, jk
@@ -185,5 +184,5 @@
       ENDIF
       !
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+      IF(ln_ctl) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('micro')")
          CALL prt_ctl_trc_info(charout)
@@ -215,32 +214,35 @@
       !!----------------------------------------------------------------------
       !
+      IF(lwp) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) 'p4z_micro_init : Initialization of microzooplankton parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampiszoo in reference namelist : Pisces microzooplankton
       READ  ( numnatp_ref, namp4zzoo, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp4zzoo in reference namelist', lwp )
-      !
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zzoo in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampiszoo in configuration namelist : Pisces microzooplankton
       READ  ( numnatp_cfg, namp4zzoo, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp4zzoo in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, namp4zzoo )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zzoo in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, namp4zzoo )
       !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for microzooplankton, namp4zzoo'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    part of calcite not dissolved in microzoo guts  part        =', part
-         WRITE(numout,*) '    microzoo preference for POC                     xpref2c     =', xpref2c
-         WRITE(numout,*) '    microzoo preference for nano                    xpref2p     =', xpref2p
-         WRITE(numout,*) '    microzoo preference for diatoms                 xpref2d     =', xpref2d
-         WRITE(numout,*) '    diatoms feeding threshold  for microzoo         xthreshdia  =', xthreshdia
-         WRITE(numout,*) '    nanophyto feeding threshold for microzoo        xthreshphy  =', xthreshphy
-         WRITE(numout,*) '    poc feeding threshold for microzoo              xthreshpoc  =', xthreshpoc
-         WRITE(numout,*) '    feeding threshold for microzooplankton          xthresh     =', xthresh
-         WRITE(numout,*) '    exsudation rate of microzooplankton             resrat      =', resrat
-         WRITE(numout,*) '    microzooplankton mortality rate                 mzrat       =', mzrat
-         WRITE(numout,*) '    maximal microzoo grazing rate                   grazrat     =', grazrat
-         WRITE(numout,*) '    non assimilated fraction of P by microzoo       unass       =', unass
-         WRITE(numout,*) '    Efficicency of microzoo growth                  epsher      =', epsher
-         WRITE(numout,*) '    Fraction of microzoo excretion as DOM           sigma1      =', sigma1
-         WRITE(numout,*) '    half sturation constant for grazing 1           xkgraz      =', xkgraz
+         WRITE(numout,*) '   Namelist : namp4zzoo'
+         WRITE(numout,*) '      part of calcite not dissolved in microzoo guts  part        =', part
+         WRITE(numout,*) '      microzoo preference for POC                     xpref2c     =', xpref2c
+         WRITE(numout,*) '      microzoo preference for nano                    xpref2p     =', xpref2p
+         WRITE(numout,*) '      microzoo preference for diatoms                 xpref2d     =', xpref2d
+         WRITE(numout,*) '      diatoms feeding threshold  for microzoo         xthreshdia  =', xthreshdia
+         WRITE(numout,*) '      nanophyto feeding threshold for microzoo        xthreshphy  =', xthreshphy
+         WRITE(numout,*) '      poc feeding threshold for microzoo              xthreshpoc  =', xthreshpoc
+         WRITE(numout,*) '      feeding threshold for microzooplankton          xthresh     =', xthresh
+         WRITE(numout,*) '      exsudation rate of microzooplankton             resrat      =', resrat
+         WRITE(numout,*) '      microzooplankton mortality rate                 mzrat       =', mzrat
+         WRITE(numout,*) '      maximal microzoo grazing rate                   grazrat     =', grazrat
+         WRITE(numout,*) '      non assimilated fraction of P by microzoo       unass       =', unass
+         WRITE(numout,*) '      Efficicency of microzoo growth                  epsher      =', epsher
+         WRITE(numout,*) '      Fraction of microzoo excretion as DOM           sigma1      =', sigma1
+         WRITE(numout,*) '      half sturation constant for grazing 1           xkgraz      =', xkgraz
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 9169)
@@ -7,13 +7,13 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
    !!----------------------------------------------------------------------
-   !!   p4z_mort       :   Compute the mortality terms for phytoplankton
-   !!   p4z_mort_init  :   Initialize the mortality params for phytoplankton
-   !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zprod         !  Primary productivity 
-   USE p4zlim          !  Phytoplankton limitation terms
-   USE prtctl_trc      !  print control for debugging
+   !!   p4z_mort       : Compute the mortality terms for phytoplankton
+   !!   p4z_mort_init  : Initialize the mortality params for phytoplankton
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE p4zprod         ! Primary productivity 
+   USE p4zlim          ! Phytoplankton limitation terms
+   USE prtctl_trc      ! print control for debugging
 
    IMPLICIT NONE
@@ -23,10 +23,9 @@
    PUBLIC   p4z_mort_init    
 
-   !! * Shared module variables
-   REAL(wp), PUBLIC :: wchl    !:
-   REAL(wp), PUBLIC :: wchld   !:
-   REAL(wp), PUBLIC :: wchldm  !:
-   REAL(wp), PUBLIC :: mprat   !:
-   REAL(wp), PUBLIC :: mprat2  !:
+   REAL(wp), PUBLIC ::   wchl     !:
+   REAL(wp), PUBLIC ::   wchld    !:
+   REAL(wp), PUBLIC ::   wchldm   !:
+   REAL(wp), PUBLIC ::   mprat    !:
+   REAL(wp), PUBLIC ::   mprat2   !:
 
    !!----------------------------------------------------------------------
@@ -35,5 +34,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -49,9 +47,9 @@
       INTEGER, INTENT(in) ::   kt ! ocean time step
       !!---------------------------------------------------------------------
-
+      !
       CALL p4z_nano            ! nanophytoplankton
-
+      !
       CALL p4z_diat            ! diatoms
-
+      !
    END SUBROUTINE p4z_mort
 
@@ -65,14 +63,14 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  :: ji, jj, jk
-      REAL(wp) :: zsizerat, zcompaph
-      REAL(wp) :: zfactfe, zfactch, zprcaca, zfracal
-      REAL(wp) :: ztortp , zrespp , zmortp 
-      CHARACTER (len=25) :: charout
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zsizerat, zcompaph
+      REAL(wp) ::   zfactfe, zfactch, zprcaca, zfracal
+      REAL(wp) ::   ztortp , zrespp , zmortp 
+      CHARACTER (len=25) ::   charout
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_nano')
       !
-      prodcal(:,:,:) = 0.  !: calcite production variable set to zero
+      prodcal(:,:,:) = 0._wp   ! calcite production variable set to zero
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -139,14 +137,13 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zfactfe,zfactsi,zfactch, zcompadi
-      REAL(wp) ::  zrespp2, ztortp2, zmortp2
-      REAL(wp) ::  zlim2, zlim1
-      CHARACTER (len=25) :: charout
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zfactfe,zfactsi,zfactch, zcompadi
+      REAL(wp) ::   zrespp2, ztortp2, zmortp2
+      REAL(wp) ::   zlim2, zlim1
+      CHARACTER (len=25) ::   charout
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_diat')
       !
-
       !    Aggregation term for diatoms is increased in case of nutrient
       !    stress as observed in reality. The stressed cells become more
@@ -196,5 +193,5 @@
       END DO
       !
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+      IF(ln_ctl) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('diat')")
          CALL prt_ctl_trc_info(charout)
@@ -214,5 +211,5 @@
       !!
       !! ** Method  :   Read the nampismort namelist and check the parameters
-      !!      called at the first timestep
+      !!              called at the first timestep
       !!
       !! ** input   :   Namelist nampismort
@@ -224,22 +221,25 @@
       !!----------------------------------------------------------------------
       !
+      IF(lwp) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) 'p4z_mort_init : Initialization of phytoplankton mortality parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampismort in reference namelist : Pisces phytoplankton
       READ  ( numnatp_ref, namp4zmort, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp4zmort in reference namelist', lwp )
-      !
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zmort in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampismort in configuration namelist : Pisces phytoplankton
       READ  ( numnatp_cfg, namp4zmort, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp4zmort in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, namp4zmort )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zmort in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, namp4zmort )
       !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for phytoplankton mortality, namp4zmort'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    quadratic mortality of phytoplankton      wchl      =', wchl
-         WRITE(numout,*) '    maximum quadratic mortality of diatoms    wchld     =', wchld
-         WRITE(numout,*) '    maximum quadratic mortality of diatoms    wchldm    =', wchldm
-         WRITE(numout,*) '    phytoplankton mortality rate              mprat     =', mprat
-         WRITE(numout,*) '    Diatoms mortality rate                    mprat2    =', mprat2
+         WRITE(numout,*) '   Namelist : namp4zmort'
+         WRITE(numout,*) '      quadratic mortality of phytoplankton        wchl   =', wchl
+         WRITE(numout,*) '      maximum quadratic mortality of diatoms      wchld  =', wchld
+         WRITE(numout,*) '      maximum quadratic mortality of diatoms      wchldm =', wchldm
+         WRITE(numout,*) '      phytoplankton mortality rate                mprat  =', mprat
+         WRITE(numout,*) '      Diatoms mortality rate                      mprat2 =', mprat2
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 9169)
@@ -4,8 +4,8 @@
    !! TOP - PISCES : Compute the light availability in the water column
    !!======================================================================
-   !! History :   1.0  !  2004     (O. Aumont) Original code
-   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
-   !!             3.2  !  2009-04  (C. Ethe, G. Madec)  optimisation
-   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Improve light availability of nano & diat
+   !! History :  1.0  !  2004     (O. Aumont) Original code
+   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!            3.2  !  2009-04  (C. Ethe, G. Madec)  optimisation
+   !!            3.4  !  2011-06  (O. Aumont, C. Ethe) Improve light availability of nano & diat
    !!----------------------------------------------------------------------
    !!   p4z_opt       : light availability in the water column
@@ -15,6 +15,6 @@
    USE sms_pisces     ! Source Minus Sink of PISCES
    USE iom            ! I/O manager
-   USE fldread         !  time interpolation
-   USE prtctl_trc      !  print control for debugging
+   USE fldread        !  time interpolation
+   USE prtctl_trc     !  print control for debugging
 
    IMPLICIT NONE
@@ -27,18 +27,18 @@
    !! * Shared module variables
 
-   LOGICAL  :: ln_varpar   !: boolean for variable PAR fraction
-   REAL(wp) :: parlux      !: Fraction of shortwave as PAR
-   REAL(wp) :: xparsw                 !: parlux/3
-   REAL(wp) :: xsi0r                 !:  1. /rn_si0
+   LOGICAL  ::   ln_varpar   ! boolean for variable PAR fraction
+   REAL(wp) ::   parlux      ! Fraction of shortwave as PAR
+   REAL(wp) ::   xparsw      ! parlux/3
+   REAL(wp) ::   xsi0r       ! 1. /rn_si0
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
    INTEGER , PARAMETER :: nbtimes = 366  !: maximum number of times record in a file
    INTEGER  :: ntimes_par                ! number of time steps in a file
-   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:) :: par_varsw    !: PAR fraction of shortwave
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ekb, ekg, ekr  !: wavelength (Red-Green-Blue)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   par_varsw      ! PAR fraction of shortwave
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ekb, ekg, ekr  ! wavelength (Red-Green-Blue)
 
    INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
 
-   REAL(wp), DIMENSION(3,61) ::   xkrgb   !: tabulated attenuation coefficients for RGB absorption
+   REAL(wp), DIMENSION(3,61) ::   xkrgb   ! tabulated attenuation coefficients for RGB absorption
    
    !!----------------------------------------------------------------------
@@ -70,10 +70,8 @@
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )  CALL timing_start('p4z_opt')
-      !
-      ! Allocate temporary workspace
-      IF( ln_p5z ) ALLOCATE( zetmp5(jpi,jpj) )
-
-      IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt )
+      IF( ln_timing )   CALL timing_start('p4z_opt')
+      IF( ln_p5z    )   ALLOCATE( zetmp5(jpi,jpj) )
+
+      IF( knt == 1 .AND. ln_varpar )   CALL p4z_opt_sbc( kt )
 
       !     Initialisation of variables used to compute PAR
@@ -84,7 +82,7 @@
       !
       !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
-                                               !  --------------------------------------------------------
-                    zchl3d(:,:,:) = trb(:,:,:,jpnch) + trb(:,:,:,jpdch)
-      IF( ln_p5z )  zchl3d(:,:,:) = zchl3d(:,:,:) + trb(:,:,:,jppch)
+      !                                        !  --------------------------------------------------------
+                     zchl3d(:,:,:) = trb(:,:,:,jpnch) + trb(:,:,:,jpdch)
+      IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)    + trb(:,:,:,jppch)
       !
       DO jk = 1, jpkm1   
@@ -105,5 +103,5 @@
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
          !
-         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
@@ -120,5 +118,5 @@
          ENDIF
          !
-         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
@@ -130,5 +128,5 @@
       ELSE
          !
-         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
@@ -240,7 +238,6 @@
       ENDIF
       !
-      IF( ln_p5z ) DEALLOCATE( zetmp5 )
-      !
-      IF( ln_timing )  CALL timing_stop('p4z_opt')
+      IF( ln_p5z    )   DEALLOCATE( zetmp5 )
+      IF( ln_timing )   CALL timing_stop('p4z_opt')
       !
    END SUBROUTINE p4z_opt
@@ -255,13 +252,12 @@
       !!
       !!----------------------------------------------------------------------
-      !! * arguments
-      INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0  
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(out)  , OPTIONAL  ::  pqsr100  
-      !! * local variables
+      INTEGER                         , INTENT(in)              ::   kt                ! ocean time-step
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   )           ::   pqsr              ! shortwave
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   pe1 , pe2 , pe3   ! PAR ( R-G-B)
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::   pe0               !
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out), OPTIONAL ::   pqsr100           !
+      !
       INTEGER    ::   ji, jj, jk     ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zqsr          !   shortwave
+      REAL(wp), DIMENSION(jpi,jpj) ::  zqsr   ! shortwave
       !!----------------------------------------------------------------------
 
@@ -272,5 +268,5 @@
       
       !  Light at the euphotic depth 
-      IF( PRESENT( pqsr100 ) )  pqsr100(:,:) = 0.01 * 3. * zqsr(:,:)
+      IF( PRESENT( pqsr100 ) )   pqsr100(:,:) = 0.01 * 3. * zqsr(:,:)
 
       IF( PRESENT( pe0 ) ) THEN     !  W-level
@@ -285,7 +281,7 @@
                DO ji = 1, jpi
                   pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t_n(ji,jj,jk-1) * xsi0r )
-                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb(ji,jj,jk-1 ) )
-                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg(ji,jj,jk-1 ) )
-                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr(ji,jj,jk-1 ) )
+                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
+                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
+                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr  (ji,jj,jk-1 )        )
                END DO
               !
@@ -327,5 +323,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER ,                INTENT(in) ::   kt     ! ocean time step
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
       !
       INTEGER  :: ji,jj
@@ -357,8 +353,6 @@
       !! ** Input   :   external ascii and netcdf files
       !!----------------------------------------------------------------------
-      INTEGER :: numpar
-      INTEGER :: ierr
-      INTEGER :: ios                 ! Local integer output status for namelist read
-      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
+      INTEGER :: numpar, ierr, ios   ! Local integer 
+      REAL(wp), DIMENSION(nbtimes) ::   zsteps   ! times records
       !
       CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
@@ -367,20 +361,21 @@
       NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux
       !!----------------------------------------------------------------------
-
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_opt_init : '
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
       REWIND( numnatp_ref )              ! Namelist nampisopt in reference namelist : Pisces attenuation coef. and PAR
       READ  ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisopt in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisopt in configuration namelist : Pisces attenuation coef. and PAR
       READ  ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisopt in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, nampisopt )
 
       IF(lwp) THEN
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' namelist : nampisopt '
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
-         WRITE(numout,*) '    PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
-         WRITE(numout,*) '    Default value for the PAR fraction   parlux         = ', parlux
+         WRITE(numout,*) '   Namelist : nampisopt '
+         WRITE(numout,*) '      PAR as a variable fraction of SW     ln_varpar      = ', ln_varpar
+         WRITE(numout,*) '      Default value for the PAR fraction   parlux         = ', parlux
       ENDIF
       !
@@ -391,6 +386,6 @@
       ! ----------------------------------------
       IF( ln_varpar ) THEN
-         IF(lwp) WRITE(numout,*) '    initialize variable par fraction '
-         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   initialize variable par fraction (ln_varpar=T)'
          !
          ALLOCATE( par_varsw(jpi,jpj) )
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 9169)
@@ -23,17 +23,16 @@
    PUBLIC   p4z_poc         ! called in p4zbio.F90
    PUBLIC   p4z_poc_init    ! called in trcsms_pisces.F90
-   PUBLIC   alngam
-   PUBLIC   gamain
-
-   !! * Shared module variables
-   REAL(wp), PUBLIC ::  xremip     !: remineralisation rate of DOC
-   REAL(wp), PUBLIC ::  xremipc    !: remineralisation rate of DOC
-   REAL(wp), PUBLIC ::  xremipn    !: remineralisation rate of DON
-   REAL(wp), PUBLIC ::  xremipp    !: remineralisation rate of DOP
-   INTEGER , PUBLIC ::  jcpoc      !: number of lability classes
-   REAL(wp), PUBLIC ::  rshape     !: shape factor of the gamma distribution
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   alphan, reminp
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: alphap
+   PUBLIC   alngam          !
+   PUBLIC   gamain          !
+
+   REAL(wp), PUBLIC ::   xremip     !: remineralisation rate of DOC
+   REAL(wp), PUBLIC ::   xremipc    !: remineralisation rate of DOC
+   REAL(wp), PUBLIC ::   xremipn    !: remineralisation rate of DON
+   REAL(wp), PUBLIC ::   xremipp    !: remineralisation rate of DOP
+   INTEGER , PUBLIC ::   jcpoc      !: number of lability classes
+   REAL(wp), PUBLIC ::   rshape     !: shape factor of the gamma distribution
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   alphan, reminp   !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   alphap           !:
 
 
@@ -53,6 +52,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step and ???
       !
       INTEGER  ::   ji, jj, jk, jn
@@ -187,6 +185,6 @@
       END DO
 
-      IF( ln_p4z ) THEN   ;  zremigoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
-      ELSE                ;  zremigoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
+      IF( ln_p4z ) THEN   ;   zremigoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
+      ELSE                ;   zremigoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
       ENDIF
 
@@ -260,7 +258,7 @@
      ! -------------------------------------------------------------------
      !
-     totprod(:,:) = 0.
+     totprod (:,:) = 0.
      totthick(:,:) = 0.
-     totcons(:,:) = 0.
+     totcons (:,:) = 0.
      ! intregrated production and consumption of POC in the mixed layer
      ! ----------------------------------------------------------------
@@ -396,6 +394,6 @@
       END DO
 
-     IF( ln_p4z ) THEN   ;  zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
-     ELSE                ;  zremipoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
+     IF( ln_p4z ) THEN   ;   zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
+     ELSE                ;   zremipoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
      ENDIF
 
@@ -473,39 +471,41 @@
       !!
       !! ** Method  :   Read the nampispoc namelist and check the parameters
-      !!      called at the first timestep
+      !!              called at the first timestep
       !!
       !! ** input   :   Namelist nampispoc
-      !!
       !!----------------------------------------------------------------------
+      INTEGER ::   jn            ! dummy loop index
       INTEGER ::   ios, ifault   ! Local integer
-      INTEGER ::   jn
-      REAL(wp) :: remindelta, reminup, remindown
+      REAL(wp)::   remindelta, reminup, remindown
       !!
       NAMELIST/nampispoc/ xremip , jcpoc  , rshape,  &
          &                xremipc, xremipn, xremipp
       !!----------------------------------------------------------------------
-
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_poc_init : Initialization of remineralization parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampisrem in reference namelist : Pisces remineralization
       READ  ( numnatp_ref, nampispoc, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampispoc in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampispoc in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisrem in configuration namelist : Pisces remineralization
       READ  ( numnatp_cfg, nampispoc, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampispoc in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampispoc )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampispoc in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampispoc )
 
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for remineralization, nampispoc'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist : nampispoc'
          IF( ln_p4z ) THEN
-            WRITE(numout,*) '    remineralisation rate of POC              xremip    =', xremip
+            WRITE(numout,*) '      remineralisation rate of POC              xremip    =', xremip
          ELSE
-            WRITE(numout,*) '    remineralisation rate of POC              xremipc   =', xremipc
-            WRITE(numout,*) '    remineralisation rate of PON              xremipn   =', xremipn
-            WRITE(numout,*) '    remineralisation rate of POP              xremipp   =', xremipp
+            WRITE(numout,*) '      remineralisation rate of POC              xremipc   =', xremipc
+            WRITE(numout,*) '      remineralisation rate of PON              xremipn   =', xremipn
+            WRITE(numout,*) '      remineralisation rate of POP              xremipp   =', xremipp
          ENDIF
-         WRITE(numout,*) '    Number of lability classes for POC        jcpoc     =', jcpoc
-         WRITE(numout,*) '    Shape factor of the gamma distribution    rshape    =', rshape
+         WRITE(numout,*) '      Number of lability classes for POC        jcpoc     =', jcpoc
+         WRITE(numout,*) '      Shape factor of the gamma distribution    rshape    =', rshape
       ENDIF
       !
@@ -513,6 +513,5 @@
       ! ---------------------------------------
       !
-      ALLOCATE( alphan(jcpoc), reminp(jcpoc) )
-      ALLOCATE( alphap(jpi,jpj,jpk,jcpoc) )
+      ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
       !
       IF (jcpoc > 1) THEN
@@ -551,38 +550,32 @@
    END SUBROUTINE p4z_poc_init
 
+
    REAL FUNCTION alngam( xvalue, ifault )
-
-!*****************************************************************************80
-!
-!! ALNGAM computes the logarithm of the gamma function.
-!
-!  Modified:
-!
-!    13 January 2008
-!
-!  Author:
-!
-!    Allan Macleod
-!    FORTRAN90 version by John Burkardt
-!
-!  Reference:
-!
-!    Allan Macleod,
-!    Algorithm AS 245,
-!    A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
-!    Applied Statistics,
-!    Volume 38, Number 2, 1989, pages 397-402.
-!
-!  Parameters:
-!
-!    Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
-!
-!    Output, integer ( kind = 4 ) IFAULT, error flag.
-!    0, no error occurred.
-!    1, XVALUE is less than or equal to 0.
-!    2, XVALUE is too big.
-!
-!    Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
-!
+      !*****************************************************************************80
+      !
+      !! ALNGAM computes the logarithm of the gamma function.
+      !
+      !  Modified:    13 January 2008
+      !
+      !  Author  :    Allan Macleod
+      !               FORTRAN90 version by John Burkardt
+      !
+      !  Reference:
+      !    Allan Macleod, Algorithm AS 245,
+      !    A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
+      !    Applied Statistics,
+      !    Volume 38, Number 2, 1989, pages 397-402.
+      !
+      !  Parameters:
+      !
+      !    Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
+      !
+      !    Output, integer ( kind = 4 ) IFAULT, error flag.
+      !    0, no error occurred.
+      !    1, XVALUE is less than or equal to 0.
+      !    2, XVALUE is too big.
+      !
+      !    Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
+      !*****************************************************************************80
   implicit none
 
@@ -746,6 +739,6 @@
    END FUNCTION alngam
 
+
    REAL FUNCTION gamain( x, p, ifault )
-
 !*****************************************************************************80
 !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 9169)
@@ -8,14 +8,14 @@
    !!             3.4  !  2011-05  (O. Aumont, C. Ethe) New parameterization of light limitation
    !!----------------------------------------------------------------------
-   !!   p4z_prod       :   Compute the growth Rate of the two phytoplanktons groups
-   !!   p4z_prod_init  :   Initialization of the parameters for growth
-   !!   p4z_prod_alloc :   Allocate variables for growth
+   !!   p4z_prod       : Compute the growth Rate of the two phytoplanktons groups
+   !!   p4z_prod_init  : Initialization of the parameters for growth
+   !!   p4z_prod_alloc : Allocate variables for growth
    !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zlim          !  Co-limitations of differents nutrients
-   USE prtctl_trc      !  print control for debugging
-   USE iom             !  I/O manager
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE p4zlim          ! Co-limitations of differents nutrients
+   USE prtctl_trc      ! print control for debugging
+   USE iom             ! I/O manager
 
    IMPLICIT NONE
@@ -26,18 +26,17 @@
    PUBLIC   p4z_prod_alloc
 
-   !! * Shared module variables
-   LOGICAL , PUBLIC ::  ln_newprod      !:
-   REAL(wp), PUBLIC ::  pislopen         !:
-   REAL(wp), PUBLIC ::  pisloped        !:
-   REAL(wp), PUBLIC ::  xadap           !:
-   REAL(wp), PUBLIC ::  excretn          !:
-   REAL(wp), PUBLIC ::  excretd         !:
-   REAL(wp), PUBLIC ::  bresp           !:
-   REAL(wp), PUBLIC ::  chlcnm          !:
-   REAL(wp), PUBLIC ::  chlcdm          !:
-   REAL(wp), PUBLIC ::  chlcmin         !:
-   REAL(wp), PUBLIC ::  fecnm           !:
-   REAL(wp), PUBLIC ::  fecdm           !:
-   REAL(wp), PUBLIC ::  grosip          !:
+   LOGICAL , PUBLIC ::   ln_newprod   !:
+   REAL(wp), PUBLIC ::   pislopen     !:
+   REAL(wp), PUBLIC ::   pisloped     !:
+   REAL(wp), PUBLIC ::   xadap        !:
+   REAL(wp), PUBLIC ::   excretn      !:
+   REAL(wp), PUBLIC ::   excretd      !:
+   REAL(wp), PUBLIC ::   bresp        !:
+   REAL(wp), PUBLIC ::   chlcnm       !:
+   REAL(wp), PUBLIC ::   chlcdm       !:
+   REAL(wp), PUBLIC ::   chlcmin      !:
+   REAL(wp), PUBLIC ::   fecnm        !:
+   REAL(wp), PUBLIC ::   fecdm        !:
+   REAL(wp), PUBLIC ::   grosip       !:
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax    !: optimal production = f(temperature)
@@ -45,7 +44,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatomee
    
-   REAL(wp) :: r1_rday                !: 1 / rday
-   REAL(wp) :: texcretn               !: 1 - excretn 
-   REAL(wp) :: texcretd               !: 1 - excretd        
+   REAL(wp) ::   r1_rday    ! 1 / rday
+   REAL(wp) ::   texcretn   ! 1 - excretn 
+   REAL(wp) ::   texcretd   ! 1 - excretd        
 
    !!----------------------------------------------------------------------
@@ -65,5 +64,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt, knt
+      INTEGER, INTENT(in) ::   kt, knt   !
       !
       INTEGER  ::   ji, jj, jk
@@ -475,7 +474,7 @@
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
      ENDIF
-     !
-     IF( ln_timing )  CALL timing_stop('p4z_prod')
-     !
+      !
+      IF( ln_timing )  CALL timing_stop('p4z_prod')
+      !
    END SUBROUTINE p4z_prod
 
@@ -492,38 +491,41 @@
       !! ** input   :   Namelist nampisprod
       !!----------------------------------------------------------------------
-      INTEGER :: ios                 ! Local integer output status for namelist read
+      INTEGER ::   ios   ! Local integer
       !
       NAMELIST/namp4zprod/ pislopen, pisloped, xadap, ln_newprod, bresp, excretn, excretd,  &
          &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
       !!----------------------------------------------------------------------
-
+      !
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_prod_init : phytoplankton growth'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampisprod in reference namelist : Pisces phytoplankton production
       READ  ( numnatp_ref, namp4zprod, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp4zprod in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zprod in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisprod in configuration namelist : Pisces phytoplankton production
       READ  ( numnatp_cfg, namp4zprod, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp4zprod in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, namp4zprod )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zprod in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, namp4zprod )
 
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for phytoplankton growth, namp4zprod'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    Enable new parame. of production (T/F)   ln_newprod    =', ln_newprod
-         WRITE(numout,*) '    mean Si/C ratio                           grosip       =', grosip
-         WRITE(numout,*) '    P-I slope                                 pislopen     =', pislopen
-         WRITE(numout,*) '    Acclimation factor to low light           xadap        =', xadap
-         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excretn      =', excretn
-         WRITE(numout,*) '    excretion ratio of diatoms                excretd      =', excretd
+         WRITE(numout,*) '   Namelist : namp4zprod'
+         WRITE(numout,*) '      Enable new parame. of production (T/F)   ln_newprod    =', ln_newprod
+         WRITE(numout,*) '      mean Si/C ratio                           grosip       =', grosip
+         WRITE(numout,*) '      P-I slope                                 pislopen     =', pislopen
+         WRITE(numout,*) '      Acclimation factor to low light           xadap        =', xadap
+         WRITE(numout,*) '      excretion ratio of nanophytoplankton      excretn      =', excretn
+         WRITE(numout,*) '      excretion ratio of diatoms                excretd      =', excretd
          IF( ln_newprod )  THEN
-            WRITE(numout,*) '    basal respiration in phytoplankton        bresp        =', bresp
-            WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
+            WRITE(numout,*) '      basal respiration in phytoplankton        bresp        =', bresp
+            WRITE(numout,*) '      Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
          ENDIF
-         WRITE(numout,*) '    P-I slope  for diatoms                    pisloped     =', pisloped
-         WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm       =', chlcnm
-         WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm       =', chlcdm
-         WRITE(numout,*) '    Maximum Fe/C in nanophytoplankton         fecnm        =', fecnm
-         WRITE(numout,*) '    Minimum Fe/C in diatoms                   fecdm        =', fecdm
+         WRITE(numout,*) '      P-I slope  for diatoms                    pisloped     =', pisloped
+         WRITE(numout,*) '      Minimum Chl/C in nanophytoplankton        chlcnm       =', chlcnm
+         WRITE(numout,*) '      Minimum Chl/C in diatoms                  chlcdm       =', chlcdm
+         WRITE(numout,*) '      Maximum Fe/C in nanophytoplankton         fecnm        =', fecnm
+         WRITE(numout,*) '      Minimum Fe/C in diatoms                   fecdm        =', fecdm
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 9169)
@@ -29,17 +29,16 @@
    PUBLIC   p4z_rem_alloc
 
-   !! * Shared module variables
-   REAL(wp), PUBLIC ::  xremikc    !: remineralisation rate of DOC 
-   REAL(wp), PUBLIC ::  xremikn    !: remineralisation rate of DON 
-   REAL(wp), PUBLIC ::  xremikp    !: remineralisation rate of DOP 
-   REAL(wp), PUBLIC ::  xremik     !: remineralisation rate of POC 
-   REAL(wp), PUBLIC ::  nitrif     !: NH4 nitrification rate 
-   REAL(wp), PUBLIC ::  xsirem     !: remineralisation rate of POC 
-   REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of POC 
-   REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica 
-   REAL(wp), PUBLIC ::  feratb     !: Fe/C quota in bacteria
-   REAL(wp), PUBLIC ::  xkferb     !: Half-saturation constant for bacteria Fe/C
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr     !: denitrification array
+   REAL(wp), PUBLIC ::   xremikc    !: remineralisation rate of DOC 
+   REAL(wp), PUBLIC ::   xremikn    !: remineralisation rate of DON 
+   REAL(wp), PUBLIC ::   xremikp    !: remineralisation rate of DOP 
+   REAL(wp), PUBLIC ::   xremik     !: remineralisation rate of POC 
+   REAL(wp), PUBLIC ::   nitrif     !: NH4 nitrification rate 
+   REAL(wp), PUBLIC ::   xsirem     !: remineralisation rate of POC 
+   REAL(wp), PUBLIC ::   xsiremlab  !: fast remineralisation rate of POC 
+   REAL(wp), PUBLIC ::   xsilab     !: fraction of labile biogenic silica 
+   REAL(wp), PUBLIC ::   feratb     !: Fe/C quota in bacteria
+   REAL(wp), PUBLIC ::   xkferb     !: Half-saturation constant for bacteria Fe/C
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr   !: denitrification array
 
    !!----------------------------------------------------------------------
@@ -303,34 +302,37 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
-
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_rem_init : Initialization of remineralization parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
       REWIND( numnatp_ref )              ! Namelist nampisrem in reference namelist : Pisces remineralization
       READ  ( numnatp_ref, nampisrem, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisrem in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisrem in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisrem in configuration namelist : Pisces remineralization
       READ  ( numnatp_cfg, nampisrem, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampisrem in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampisrem )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisrem in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampisrem )
 
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameters for remineralization, nampisrem'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist parameters for remineralization, nampisrem'
          IF( ln_p4z ) THEN
-            WRITE(numout,*) '    remineralization rate of DOC              xremik    =', xremik
+            WRITE(numout,*) '      remineralization rate of DOC              xremik    =', xremik
          ELSE
-            WRITE(numout,*) '    remineralization rate of DOC              xremikc   =', xremikc
-            WRITE(numout,*) '    remineralization rate of DON              xremikn   =', xremikn
-            WRITE(numout,*) '    remineralization rate of DOP              xremikp   =', xremikp
+            WRITE(numout,*) '      remineralization rate of DOC              xremikc   =', xremikc
+            WRITE(numout,*) '      remineralization rate of DON              xremikn   =', xremikn
+            WRITE(numout,*) '      remineralization rate of DOP              xremikp   =', xremikp
          ENDIF
-         WRITE(numout,*) '    remineralization rate of Si               xsirem    =', xsirem
-         WRITE(numout,*) '    fast remineralization rate of Si          xsiremlab =', xsiremlab
-         WRITE(numout,*) '    fraction of labile biogenic silica        xsilab    =', xsilab
-         WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif
-         WRITE(numout,*) '    Bacterial Fe/C ratio                      feratb    =', feratb
-         WRITE(numout,*) '    Half-saturation constant for bact. Fe/C   xkferb    =', xkferb
+         WRITE(numout,*) '      remineralization rate of Si               xsirem    =', xsirem
+         WRITE(numout,*) '      fast remineralization rate of Si          xsiremlab =', xsiremlab
+         WRITE(numout,*) '      fraction of labile biogenic silica        xsilab    =', xsilab
+         WRITE(numout,*) '      NH4 nitrification rate                    nitrif    =', nitrif
+         WRITE(numout,*) '      Bacterial Fe/C ratio                      feratb    =', feratb
+         WRITE(numout,*) '      Half-saturation constant for bact. Fe/C   xkferb    =', xkferb
       ENDIF
       !
-      denitr  (:,:,:) = 0._wp
+      denitr(:,:,:) = 0._wp
       !
    END SUBROUTINE p4z_rem_init
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 9169)
@@ -21,28 +21,26 @@
    PUBLIC   p4z_sbc_init   
 
-   LOGICAL , PUBLIC  :: ln_dust     !: boolean for dust input from the atmosphere
-   LOGICAL , PUBLIC  :: ln_solub    !: boolean for variable solubility of atmospheric iron
-   LOGICAL , PUBLIC  :: ln_river    !: boolean for river input of nutrients
-   LOGICAL , PUBLIC  :: ln_ndepo    !: boolean for atmospheric deposition of N
-   LOGICAL , PUBLIC  :: ln_ironsed  !: boolean for Fe input from sediments
-   LOGICAL , PUBLIC  :: ln_hydrofe  !: boolean for Fe input from hydrothermal vents
-   LOGICAL , PUBLIC  :: ln_ironice  !: boolean for Fe input from sea ice
-   REAL(wp), PUBLIC  :: sedfeinput  !: Coastal release of Iron
-   REAL(wp), PUBLIC  :: dustsolub   !: Solubility of the dust
-   REAL(wp), PUBLIC  :: mfrac       !: Mineral Content of the dust
-   REAL(wp), PUBLIC  :: icefeinput  !: Iron concentration in sea ice
-   REAL(wp), PUBLIC  :: wdust       !: Sinking speed of the dust 
-   REAL(wp), PUBLIC  :: nitrfix     !: Nitrogen fixation rate   
-   REAL(wp), PUBLIC  :: diazolight  !: Nitrogen fixation sensitivty to light 
-   REAL(wp), PUBLIC  :: concfediaz  !: Fe half-saturation Cste for diazotrophs 
-   REAL(wp)          :: hratio      !: Fe:3He ratio assumed for vent iron supply
-   REAL(wp), PUBLIC  :: fep_rats    !: Fep/Fer ratio from sed  sources
-   REAL(wp), PUBLIC  :: fep_rath    !: Fep/Fer ratio from hydro sources
-   REAL(wp), PUBLIC  :: lgw_rath    !: Weak ligand ratio from hydro sources
-
-
-   LOGICAL , PUBLIC  :: ll_sbc
-
-   LOGICAL  ::  ll_solub
+   LOGICAL , PUBLIC ::   ln_dust      !: boolean for dust input from the atmosphere
+   LOGICAL , PUBLIC ::   ln_solub     !: boolean for variable solubility of atmospheric iron
+   LOGICAL , PUBLIC ::   ln_river     !: boolean for river input of nutrients
+   LOGICAL , PUBLIC ::   ln_ndepo     !: boolean for atmospheric deposition of N
+   LOGICAL , PUBLIC ::   ln_ironsed   !: boolean for Fe input from sediments
+   LOGICAL , PUBLIC ::   ln_hydrofe   !: boolean for Fe input from hydrothermal vents
+   LOGICAL , PUBLIC ::   ln_ironice   !: boolean for Fe input from sea ice
+   REAL(wp), PUBLIC ::   sedfeinput   !: Coastal release of Iron
+   REAL(wp), PUBLIC ::   dustsolub    !: Solubility of the dust
+   REAL(wp), PUBLIC ::   mfrac        !: Mineral Content of the dust
+   REAL(wp), PUBLIC ::   icefeinput   !: Iron concentration in sea ice
+   REAL(wp), PUBLIC ::   wdust        !: Sinking speed of the dust 
+   REAL(wp), PUBLIC ::   nitrfix      !: Nitrogen fixation rate   
+   REAL(wp), PUBLIC ::   diazolight   !: Nitrogen fixation sensitivty to light 
+   REAL(wp), PUBLIC ::   concfediaz   !: Fe half-saturation Cste for diazotrophs 
+   REAL(wp)         ::   hratio       !: Fe:3He ratio assumed for vent iron supply
+   REAL(wp), PUBLIC ::   fep_rats     !: Fep/Fer ratio from sed  sources
+   REAL(wp), PUBLIC ::   fep_rath     !: Fep/Fer ratio from hydro sources
+   REAL(wp), PUBLIC ::   lgw_rath     !: Weak ligand ratio from hydro sources
+
+   LOGICAL , PUBLIC ::   ll_sbc
+   LOGICAL          ::   ll_solub
 
    INTEGER , PARAMETER  :: jpriv  = 7   !: Maximum number of river input fields
@@ -55,28 +53,27 @@
    INTEGER , PARAMETER  :: jr_dsi = 7   !: index of dissolved silicate
 
-
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_solub      ! structure of input dust
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_river  ! structure of input riverdic
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_solub     ! structure of input dust
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_river     ! structure of input riverdic
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ndepo     ! structure of input nitrogen deposition
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_hydrofe   ! structure of input iron from hydrothermal vents
 
-   INTEGER , PARAMETER :: nbtimes = 365  !: maximum number of times record in a file
-   INTEGER  :: ntimes_dust, ntimes_riv, ntimes_ndep       ! number of time steps in a file
-   INTEGER  :: ntimes_solub, ntimes_hydro                 ! number of time steps in a file
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: dust, solub       !: dust fields
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdic, rivalk    !: river input fields
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdin, rivdip    !: river input fields
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdon, rivdop    !: river input fields
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdoc    !: river input fields
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdsi    !: river input fields
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: nitdep    !: atmospheric N deposition 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ironsed   !: Coastal supply of iron
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hydrofe   !: Hydrothermal vent supply of iron
-
-   REAL(wp), PUBLIC :: sumdepsi, rivalkinput, rivdicinput, nitdepinput
-   REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput
+   INTEGER , PARAMETER ::   nbtimes = 365                          ! maximum number of times record in a file
+   INTEGER             ::   ntimes_dust, ntimes_riv, ntimes_ndep   ! number of time steps in a file
+   INTEGER             ::   ntimes_solub, ntimes_hydro             ! number of time steps in a file
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dust  , solub    !: dust fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdic, rivalk   !: river input fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdin, rivdip   !: river input fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdon, rivdop   !: river input fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdoc           !: river input fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdsi           !: river input fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nitdep           !: atmospheric N deposition 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ironsed          !: Coastal supply of iron
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hydrofe          !: Hydrothermal vent supply of iron
+
+   REAL(wp), PUBLIC ::   rivalkinput, rivdicinput, nitdepinput, sumdepsi
+   REAL(wp), PUBLIC ::   rivdininput, rivdipinput, rivdsiinput
 
    !! * Substitutions
@@ -100,13 +97,11 @@
       !!
       !!----------------------------------------------------------------------
-      !! * arguments
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
-
-      !! * local declarations
-      INTEGER  :: ji,jj 
-      REAL(wp) :: zcoef, zyyss
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !
+      INTEGER  ::   ji, jj 
+      REAL(wp) ::   zcoef, zyyss
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )  CALL timing_start('p4z_sbc')
+      IF( ln_timing )   CALL timing_start('p4z_sbc')
       !
       ! Compute dust at nit000 or only if there is more than 1 time record in dust file
@@ -114,12 +109,10 @@
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
             CALL fld_read( kt, 1, sf_dust )
-            IF( nn_ice_tr == -1 .AND. .NOT. ln_ironice ) THEN
-               dust(:,:) = sf_dust(1)%fnow(:,:,1)
-            ELSE
-               dust(:,:) = sf_dust(1)%fnow(:,:,1) * ( 1.0 - fr_i(:,:) )
+            IF( nn_ice_tr == -1 .AND. .NOT.ln_ironice ) THEN   ;   dust(:,:) = sf_dust(1)%fnow(:,:,1)
+            ELSE                                               ;   dust(:,:) = sf_dust(1)%fnow(:,:,1) * ( 1.-fr_i(:,:) )
             ENDIF
          ENDIF
       ENDIF
-
+      !
       IF( ll_solub ) THEN
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_solub > 1 ) ) THEN
@@ -205,5 +198,4 @@
       !!
       !!----------------------------------------------------------------------
-      !
       INTEGER  :: ji, jj, jk, jm, ifpr
       INTEGER  :: ii0, ii1, ij0, ij1
@@ -224,5 +216,5 @@
       TYPE(FLD_N) ::   sn_riverdoc, sn_riverdic, sn_riverdsi   ! informations about the fields to be read
       TYPE(FLD_N) ::   sn_riverdin, sn_riverdon, sn_riverdip, sn_riverdop
-      !
+      !!
       NAMELIST/nampissbc/cn_dir, sn_dust, sn_solub, sn_riverdic, sn_riverdoc, sn_riverdin, sn_riverdon,     &
         &                sn_riverdip, sn_riverdop, sn_riverdsi, sn_ndepo, sn_ironsed, sn_hydrofe, &
@@ -232,56 +224,58 @@
       !!----------------------------------------------------------------------
       !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_sbc_init : initialization of the external sources of nutrients '
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
       !                            !* set file information
       REWIND( numnatp_ref )              ! Namelist nampissbc in reference namelist : Pisces external sources of nutrients
       READ  ( numnatp_ref, nampissbc, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampissbc in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampissbc in configuration namelist : Pisces external sources of nutrients
       READ  ( numnatp_cfg, nampissbc, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampissbc in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, nampissbc )
 
-      IF ( ( nn_ice_tr >= 0 ) .AND. ln_ironice ) THEN
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist : nampissbc '
+         WRITE(numout,*) '      dust input from the atmosphere           ln_dust     = ', ln_dust
+         WRITE(numout,*) '      Variable solubility of iron input        ln_solub    = ', ln_solub
+         WRITE(numout,*) '      river input of nutrients                 ln_river    = ', ln_river
+         WRITE(numout,*) '      atmospheric deposition of n              ln_ndepo    = ', ln_ndepo
+         WRITE(numout,*) '      Fe input from sediments                  ln_ironsed  = ', ln_ironsed
+         WRITE(numout,*) '      Fe input from seaice                     ln_ironice  = ', ln_ironice
+         WRITE(numout,*) '      fe input from hydrothermal vents         ln_hydrofe  = ', ln_hydrofe
+         WRITE(numout,*) '      coastal release of iron                  sedfeinput  = ', sedfeinput
+         WRITE(numout,*) '      solubility of the dust                   dustsolub   = ', dustsolub
+         WRITE(numout,*) '      Mineral Fe content of the dust           mfrac       = ', mfrac
+         WRITE(numout,*) '      Iron concentration in sea ice            icefeinput  = ', icefeinput
+         WRITE(numout,*) '      sinking speed of the dust                wdust       = ', wdust
+         WRITE(numout,*) '      nitrogen fixation rate                   nitrfix     = ', nitrfix
+         WRITE(numout,*) '      nitrogen fixation sensitivty to light    diazolight  = ', diazolight
+         WRITE(numout,*) '      Fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
+         WRITE(numout,*) '      Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
+         IF( ln_ligand ) THEN
+            WRITE(numout,*) '      Fep/Fer ratio from sed sources            fep_rats   = ', fep_rats
+            WRITE(numout,*) '      Fep/Fer ratio from sed hydro sources      fep_rath   = ', fep_rath
+            WRITE(numout,*) '      Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
+         ENDIF
+      END IF
+
+      IF( nn_ice_tr >= 0 .AND. ln_ironice ) THEN
          IF(lwp) THEN
-            WRITE(numout,*) ' ln_ironice incompatible with nn_ice_tr = ', nn_ice_tr
-            WRITE(numout,*) ' Specify your sea ice iron concentration in nampisice instead '
-            WRITE(numout,*) ' ln_ironice is forced to .FALSE. '
-            ln_ironice = .FALSE.
-         ENDIF
-      ENDIF
-
-      IF(lwp) THEN
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' namelist : nampissbc '
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
-         WRITE(numout,*) '    dust input from the atmosphere           ln_dust     = ', ln_dust
-         WRITE(numout,*) '    Variable solubility of iron input        ln_solub    = ', ln_solub
-         WRITE(numout,*) '    river input of nutrients                 ln_river    = ', ln_river
-         WRITE(numout,*) '    atmospheric deposition of n              ln_ndepo    = ', ln_ndepo
-         WRITE(numout,*) '    Fe input from sediments                  ln_ironsed  = ', ln_ironsed
-         WRITE(numout,*) '    Fe input from seaice                     ln_ironice  = ', ln_ironice
-         WRITE(numout,*) '    fe input from hydrothermal vents         ln_hydrofe  = ', ln_hydrofe
-         WRITE(numout,*) '    coastal release of iron                  sedfeinput  = ', sedfeinput
-         WRITE(numout,*) '    solubility of the dust                   dustsolub   = ', dustsolub
-         WRITE(numout,*) '    Mineral Fe content of the dust           mfrac       = ', mfrac
-         WRITE(numout,*) '    Iron concentration in sea ice            icefeinput  = ', icefeinput
-         WRITE(numout,*) '    sinking speed of the dust                wdust       = ', wdust
-         WRITE(numout,*) '    nitrogen fixation rate                   nitrfix     = ', nitrfix
-         WRITE(numout,*) '    nitrogen fixation sensitivty to light    diazolight  = ', diazolight
-         WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
-         WRITE(numout,*) '    Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
-         IF( ln_ligand ) THEN
-            WRITE(numout,*) '    Fep/Fer ratio from sed sources            fep_rats   = ', fep_rats
-            WRITE(numout,*) '    Fep/Fer ratio from sed hydro sources      fep_rath   = ', fep_rath
-            WRITE(numout,*) '    Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
-         ENDIF
-      END IF
-
-      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN  ;  ll_sbc = .TRUE.
-      ELSE                                            ;  ll_sbc = .FALSE.
-      ENDIF
-
-      IF( ln_dust .AND. ln_solub ) THEN               ;  ll_solub = .TRUE.
-      ELSE                                            ;  ll_solub = .FALSE.
+            WRITE(numout,*) '   ==>>>   ln_ironice incompatible with nn_ice_tr = ', nn_ice_tr
+            WRITE(numout,*) '           Specify your sea ice iron concentration in nampisice instead '
+            WRITE(numout,*) '           ln_ironice is forced to .FALSE. '
+         ENDIF
+         ln_ironice = .FALSE.
+      ENDIF
+
+      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN   ;   ll_sbc = .TRUE.
+      ELSE                                             ;   ll_sbc = .FALSE.
+      ENDIF
+
+      IF( ln_dust .AND. ln_solub ) THEN                ;   ll_solub = .TRUE.
+      ELSE                                             ;   ll_solub = .FALSE.
       ENDIF
 
@@ -322,5 +316,5 @@
             DO jm = 1, ntimes_dust
                sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) * ztimes_dust )
-            ENDDO
+            END DO
             sumdepsi = sumdepsi / ( nyear_len(1) * rday ) * 12. * 8.8 * 0.075 * mfrac / 28.1 
             DEALLOCATE( zdust)
@@ -335,6 +329,6 @@
       IF( ll_solub ) THEN
          !
-         IF(lwp) WRITE(numout,*) '    initialize variable solubility of Fe '
-         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   ll_solub=T , initialize variable solubility of Fe '
          !
          ALLOCATE( solub(jpi,jpj) )    ! allocation
@@ -356,6 +350,6 @@
       IF( ln_river ) THEN
          !
-         slf_river(jr_dic) = sn_riverdic  ;  slf_river(jr_doc) = sn_riverdoc  ;  slf_river(jr_din) = sn_riverdin 
-         slf_river(jr_don) = sn_riverdon  ;  slf_river(jr_dip) = sn_riverdip  ;  slf_river(jr_dop) = sn_riverdop
+         slf_river(jr_dic) = sn_riverdic   ;   slf_river(jr_doc) = sn_riverdoc   ;   slf_river(jr_din) = sn_riverdin 
+         slf_river(jr_don) = sn_riverdon   ;   slf_river(jr_dip) = sn_riverdip   ;   slf_river(jr_dop) = sn_riverdop
          slf_river(jr_dsi) = sn_riverdsi  
          !
@@ -363,6 +357,6 @@
          IF( ln_p5z )  ALLOCATE( rivdon(jpi,jpj), rivdop(jpi,jpj), rivdoc(jpi,jpj) )
          !
-         ALLOCATE( sf_river(jpriv), rivinput(jpriv), STAT=ierr1 )           !* allocate and fill sf_river (forcing structure) with sn_river_
-         rivinput(:) = 0.0
+         ALLOCATE( sf_river(jpriv), rivinput(jpriv), STAT=ierr1 )    !* allocate and fill sf_river (forcing structure) with sn_river_
+         rivinput(:) = 0._wp
 
          IF( ierr1 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_irver structure' )
@@ -409,6 +403,6 @@
       IF( ln_ndepo ) THEN
          !
-         IF(lwp) WRITE(numout,*) '    initialize the nutrient input by dust from ndeposition.orca.nc'
-         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ndepo=T , initialize the nutrient input by dust from NetCDF file'
          !
          ALLOCATE( nitdep(jpi,jpj) )    ! allocation
@@ -446,6 +440,6 @@
       IF( ln_ironsed ) THEN     
          !
-         IF(lwp) WRITE(numout,*) '    computation of an island mask to enhance coastal supply of iron'
-         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
          !
          ALLOCATE( ironsed(jpi,jpj,jpk) )    ! allocation
@@ -458,9 +452,8 @@
          ik50 = 5        !  last level where depth less than 50 m
          DO jk = jpkm1, 1, -1
-            IF( gdept_1d(jk) > 50. )  ik50 = jk - 1
+            IF( gdept_1d(jk) > 50. )   ik50 = jk - 1
          END DO
-         IF (lwp) WRITE(numout,*)
-         IF (lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
-         IF (lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
          DO jk = 1, ik50
             DO jj = 2, jpjm1
@@ -499,6 +492,6 @@
       IF( ln_hydrofe ) THEN
          !
-         IF(lwp) WRITE(numout,*) '    Input of iron from hydrothermal vents '
-         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   ln_hydrofe=T , Input of iron from hydrothermal vents'
          !
          ALLOCATE( hydrofe(jpi,jpj,jpk) )    ! allocation
@@ -521,8 +514,8 @@
          WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '    N Supply   : ', rivdininput*rno3*1E3/1E12*14.,' TgN/yr'
-         WRITE(numout,*) '    Si Supply  : ', rivdsiinput*1E3/1E12*28.1,' TgSi/yr'
+         WRITE(numout,*) '    Si Supply  : ', rivdsiinput*1E3/1E12*28.1    ,' TgSi/yr'
          WRITE(numout,*) '    P Supply   : ', rivdipinput*1E3*po4r/1E12*31.,' TgP/yr'
-         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12,' Teq/yr'
-         WRITE(numout,*) '    DIC Supply : ', rivdicinput*1E3*12./1E12,'TgC/yr'
+         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12         ,' Teq/yr'
+         WRITE(numout,*) '    DIC Supply : ', rivdicinput*1E3*12./1E12     ,' TgC/yr'
          WRITE(numout,*) 
          WRITE(numout,*) '    Total input of elements from atmospheric supply'
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 9169)
@@ -7,23 +7,23 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
    !!----------------------------------------------------------------------
-   !!   p4zsms         :  Time loop of passive tracers sms
+   !!   p4z_sms        : Time loop of passive tracers sms
    !!----------------------------------------------------------------------
-   USE oce_trc         !  shared variables between ocean and passive tracers
-   USE trc             !  passive tracers common variables 
-   USE trcdta
-   USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zbio          !  Biological model
-   USE p4zche          !  Chemical model
-   USE p4zlys          !  Calcite saturation
-   USE p4zflx          !  Gas exchange
-   USE p4zsbc          !  External source of nutrients
-   USE p4zsed          !  Sedimentation
-   USE p4zint          !  time interpolation
-   USE p4zrem          !  remineralisation
-   USE iom             !  I/O manager
-   USE trd_oce         !  Ocean trends variables
-   USE trdtrc          !  TOP trends variables
-   USE sedmodel        !  Sediment model
-   USE prtctl_trc      !  print control for debugging
+   USE oce_trc         ! shared variables between ocean and passive tracers
+   USE trc             ! passive tracers common variables 
+   USE trcdta          ! 
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE p4zbio          ! Biological model
+   USE p4zche          ! Chemical model
+   USE p4zlys          ! Calcite saturation
+   USE p4zflx          ! Gas exchange
+   USE p4zsbc          ! External source of nutrients
+   USE p4zsed          ! Sedimentation
+   USE p4zint          ! time interpolation
+   USE p4zrem          ! remineralisation
+   USE iom             ! I/O manager
+   USE trd_oce         ! Ocean trends variables
+   USE trdtrc          ! TOP trends variables
+   USE sedmodel        ! Sediment model
+   USE prtctl_trc      ! print control for debugging
 
    IMPLICIT NONE
@@ -33,10 +33,9 @@
    PUBLIC   p4z_sms        ! called in p4zsms.F90
 
-   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget
-   REAL(wp) :: xfact1, xfact2, xfact3
-   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget
-
-   !!* Array used to indicate negative tracer values
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ???
+   INTEGER ::    numco2, numnut, numnit      ! logical unit for co2 budget
+   REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget
+   REAL(wp) ::   xfact1, xfact2, xfact3
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     ! Array used to indicate negative tracer values
 
    !!----------------------------------------------------------------------
@@ -45,5 +44,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -197,35 +195,40 @@
       NAMELIST/nampismass/ ln_check_mass
       !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'p4z_sms_init : PISCES initialization'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
 
       REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
       READ  ( numnatp_ref, nampisbio, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisbio in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
       READ  ( numnatp_cfg, nampisbio, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampisbio in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampisbio )
-
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisbio in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampisbio )
+      !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' Namelist : nampisbio'
-         WRITE(numout,*) '    frequence pour la biologie                nrdttrc    =', nrdttrc
-         WRITE(numout,*) '    POC sinking speed                         wsbio      =', wsbio
-         WRITE(numout,*) '    half saturation constant for mortality    xkmort     =', xkmort 
+         WRITE(numout,*) '   Namelist : nampisbio'
+         WRITE(numout,*) '      frequency for the biology                 nrdttrc     =', nrdttrc
+         WRITE(numout,*) '      POC sinking speed                         wsbio       =', wsbio
+         WRITE(numout,*) '      half saturation constant for mortality    xkmort      =', xkmort 
          IF( ln_p5z ) THEN
-            WRITE(numout,*) '    N/C in zooplankton                        no3rat3    =', no3rat3
-            WRITE(numout,*) '    P/C in zooplankton                        po4rat3    =', po4rat3
-         ENDIF
-         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3     =', ferat3
-         WRITE(numout,*) '    Big particles sinking speed               wsbio2     =', wsbio2
-         WRITE(numout,*) '    Big particles maximum sinking speed       wsbio2max  =', wsbio2max
-         WRITE(numout,*) '    Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
-         WRITE(numout,*) '    Maximum number of iterations for POC      niter1max =', niter1max
-         WRITE(numout,*) '    Maximum number of iterations for GOC      niter2max =', niter2max
+            WRITE(numout,*) '      N/C in zooplankton                        no3rat3     =', no3rat3
+            WRITE(numout,*) '      P/C in zooplankton                        po4rat3     =', po4rat3
+         ENDIF
+         WRITE(numout,*) '      Fe/C in zooplankton                       ferat3      =', ferat3
+         WRITE(numout,*) '      Big particles sinking speed               wsbio2      =', wsbio2
+         WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max
+         WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
+         WRITE(numout,*) '      Maximum number of iterations for POC      niter1max   =', niter1max
+         WRITE(numout,*) '      Maximum number of iterations for GOC      niter2max   =', niter2max
          IF( ln_ligand ) THEN
-            WRITE(numout,*) '    FeP sinking speed                             wfep   =', wfep
+            WRITE(numout,*) '      FeP sinking speed                              wfep   =', wfep
             IF( ln_p4z ) THEN
-              WRITE(numout,*) '    Phyto ligand production per unit doc          ldocp  =', ldocp
-              WRITE(numout,*) '    Zoo ligand production per unit doc            ldocz  =', ldocz
-              WRITE(numout,*) '    Proportional loss of ligands due to Fe uptake lthet  =', lthet
+               WRITE(numout,*) '      Phyto ligand production per unit doc           ldocp  =', ldocp
+               WRITE(numout,*) '      Zoo ligand production per unit doc             ldocz  =', ldocz
+               WRITE(numout,*) '      Proportional loss of ligands due to Fe uptake  lthet  =', lthet
             ENDIF
          ENDIF
@@ -235,33 +238,29 @@
       REWIND( numnatp_ref )              ! Namelist nampisdmp in reference namelist : Pisces damping
       READ  ( numnatp_ref, nampisdmp, IOSTAT = ios, ERR = 905)
-905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in reference namelist', lwp )
-
+905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisdmp in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisdmp in configuration namelist : Pisces damping
       READ  ( numnatp_cfg, nampisdmp, IOSTAT = ios, ERR = 906 )
-906   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampisdmp in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampisdmp )
-
+906   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisdmp in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampisdmp )
+      !
       IF(lwp) THEN                         ! control print
          WRITE(numout,*)
-         WRITE(numout,*) ' Namelist : nampisdmp'
-         WRITE(numout,*) '    Relaxation of tracer to glodap mean value             ln_pisdmp      =', ln_pisdmp
-         WRITE(numout,*) '    Frequency of Relaxation                               nn_pisdmp      =', nn_pisdmp
-         WRITE(numout,*) ' '
+         WRITE(numout,*) '   Namelist : nampisdmp --- relaxation to GLODAP'
+         WRITE(numout,*) '      Relaxation of tracer to glodap mean value   ln_pisdmp =', ln_pisdmp
+         WRITE(numout,*) '      Frequency of Relaxation                     nn_pisdmp =', nn_pisdmp
       ENDIF
 
       REWIND( numnatp_ref )              ! Namelist nampismass in reference namelist : Pisces mass conservation check
       READ  ( numnatp_ref, nampismass, IOSTAT = ios, ERR = 907)
-907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in reference namelist', lwp )
-
+907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampismass in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampismass in configuration namelist : Pisces mass conservation check 
       READ  ( numnatp_cfg, nampismass, IOSTAT = ios, ERR = 908 )
-908   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampismass in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampismass )
+908   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampismass in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampismass )
 
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Namelist parameter for mass conservation checking'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    Flag to check mass conservation of NO3/Si/TALK ln_check_mass = ', ln_check_mass
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist : nampismass  --- mass conservation checking'
+         WRITE(numout,*) '      Flag to check mass conservation of NO3/Si/TALK   ln_check_mass = ', ln_check_mass
       ENDIF
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 9169)
@@ -37,5 +37,4 @@
       !! ** Purpose :   Initialisation of the PISCES biochemical model
       !!----------------------------------------------------------------------
-
       !
       CALL trc_nam_pisces
@@ -46,4 +45,5 @@
 
    END SUBROUTINE trc_ini_pisces
+
 
    SUBROUTINE p4z_ini
@@ -53,5 +53,4 @@
       !! ** Purpose :   Initialisation of the PISCES biochemical model
       !!----------------------------------------------------------------------
-      !
       USE p4zsms          ! Main P4Z routine
       USE p4zche          !  Chemical model
@@ -76,13 +75,12 @@
       USE p5zmeso         !  Sources and sinks of mesozooplankton
       USE p5zmort         !  Mortality terms for phytoplankton
-
-      !
-      REAL(wp), SAVE :: sco2   =  2.312e-3_wp
-      REAL(wp), SAVE :: alka0  =  2.426e-3_wp
-      REAL(wp), SAVE :: oxyg0  =  177.6e-6_wp 
-      REAL(wp), SAVE :: po4    =  2.165e-6_wp 
-      REAL(wp), SAVE :: bioma0 =  1.000e-8_wp  
-      REAL(wp), SAVE :: silic1 =  91.51e-6_wp  
-      REAL(wp), SAVE :: no3    =  30.9e-6_wp * 7.625_wp
+      !
+      REAL(wp), SAVE ::   sco2   =  2.312e-3_wp
+      REAL(wp), SAVE ::   alka0  =  2.426e-3_wp
+      REAL(wp), SAVE ::   oxyg0  =  177.6e-6_wp 
+      REAL(wp), SAVE ::   po4    =  2.165e-6_wp 
+      REAL(wp), SAVE ::   bioma0 =  1.000e-8_wp  
+      REAL(wp), SAVE ::   silic1 =  91.51e-6_wp  
+      REAL(wp), SAVE ::   no3    =  30.9e-6_wp * 7.625_wp
       !
       INTEGER  ::  ji, jj, jk, jn, ierr
@@ -90,16 +88,15 @@
       REAL(wp) ::  ztmas, ztmas1
       CHARACTER(len = 20)  ::  cltra
-
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+      !
       IF(lwp) THEN
          WRITE(numout,*)
          IF( ln_p4z ) THEN 
-            WRITE(numout,*) ' p4z_ini :   PISCES biochemical model initialisation'
+            WRITE(numout,*) 'p4z_ini :   PISCES biochemical model initialisation'
+            WRITE(numout,*) '~~~~~~~'
          ELSE
-            WRITE(numout,*) ' p5z_ini :   PISCES biochemical model initialisation'
-            WRITE(numout,*) '             With variable stoichiometry'
+            WRITE(numout,*) 'p5z_ini :   PISCES biochemical model initialisation'
+            WRITE(numout,*) '~~~~~~~     With variable stoichiometry'
          ENDIF
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~'
       ENDIF
       !
@@ -170,8 +167,8 @@
         IF( cltra == 'LGW'      )   jplgw = jn      !: Weak ligands
         IF( cltra == 'LFe'      )   jpfep = jn      !: Fe nanoparticle
-      ENDDO
+      END DO
 
       CALL p4z_sms_init       !  Maint routine
-      !                                            ! Time-step
+      !
 
       ! Set biological ratios
@@ -275,8 +272,9 @@
 
       IF(lwp) WRITE(numout,*) 
-      IF(lwp) WRITE(numout,*) 'Initialization of PISCES tracers done'
+      IF(lwp) WRITE(numout,*) '   ==>>>   Initialization of PISCES tracers done'
       IF(lwp) WRITE(numout,*) 
       !
    END SUBROUTINE p4z_ini
+
 
    SUBROUTINE p2z_ini
@@ -298,5 +296,5 @@
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) ' p2z_ini :   LOBSTER biochemical model initialisation'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~'
 
       ierr =        sms_pisces_alloc()          
@@ -346,5 +344,5 @@
       !
       IF(lwp) WRITE(numout,*) 
-      IF(lwp) WRITE(numout,*) 'Initialization of LOBSTER tracers done'
+      IF(lwp) WRITE(numout,*) '   ==>>>   Initialization of LOBSTER tracers done'
       IF(lwp) WRITE(numout,*) 
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90	(revision 9125)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90	(revision 9169)
@@ -9,5 +9,5 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trcnam.pisces.h90
    !!----------------------------------------------------------------------
-   !! trc_nam_pisces       : PISCES model namelist read
+   !! trc_nam_pisces   : PISCES model namelist read
    !!----------------------------------------------------------------------
    USE oce_trc         ! Ocean variables
@@ -40,6 +40,6 @@
       !!----------------------------------------------------------------------
       INTEGER :: jl, jn
-      INTEGER :: ios, ioptio                 ! Local integer output status for namelist read
-      CHARACTER(LEN=20)   ::   clname
+      INTEGER :: ios, ioptio         ! Local integer
+      CHARACTER(LEN=20)::   clname
       !!
       NAMELIST/nampismod/ln_p2z, ln_p4z, ln_p5z, ln_ligand
@@ -49,36 +49,32 @@
       clname = 'namelist_pisces'
 
-      IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read PISCES namelist'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) 'trc_nam_pisces : read PISCES namelist'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
       CALL ctl_opn( numnatp_ref, TRIM( clname )//'_ref', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
       CALL ctl_opn( numnatp_cfg, TRIM( clname )//'_cfg', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
       IF(lwm) CALL ctl_opn( numonp     , 'output.namelist.pis' , 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
       !
-
       REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
       READ  ( numnatp_ref, nampismod, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismod in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampismod in reference namelist', lwp )
       REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
       READ  ( numnatp_cfg, nampismod, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismod in configuration namelist', lwp )
-      IF(lwm) WRITE ( numonp, nampismod )
-
-     IF(lwp) THEN                  ! control print
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' Flag to use LOBSTER model            ln_p2z    = ', ln_p2z
-         WRITE(numout,*) ' Flag to use PISCES standard  model   ln_p4z    = ', ln_p4z
-         WRITE(numout,*) ' Flag to use PISCES quota     model   ln_p5z    = ', ln_p5z
-         WRITE(numout,*) ' Flag to ligand                       ln_ligand = ', ln_ligand
-         WRITE(numout,*) ' '
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampismod in configuration namelist', lwp )
+      IF(lwm) WRITE( numonp, nampismod )
+      !
+      IF(lwp) THEN                  ! control print
+         WRITE(numout,*) '   Namelist : nampismod '
+         WRITE(numout,*) '      Flag to use LOBSTER model            ln_p2z    = ', ln_p2z
+         WRITE(numout,*) '      Flag to use PISCES standard model    ln_p4z    = ', ln_p4z
+         WRITE(numout,*) '      Flag to use PISCES quota    model    ln_p5z    = ', ln_p5z
+         WRITE(numout,*) '      Flag to ligand                       ln_ligand = ', ln_ligand
       ENDIF
-
+      !
       IF(lwp) THEN                         ! control print
-         WRITE(numout,*) ' '
-         IF( ln_p5z    )  WRITE(numout,*) '  PISCES QUOTA model is used'
-         IF( ln_p4z    )  WRITE(numout,*) '  PISCES STANDARD model is used'
-         IF( ln_p2z    )  WRITE(numout,*) '  LOBSTER model is used'
-         IF( ln_ligand )  WRITE(numout,*) '  Compute remineralization/dissolution of organic ligands'
-         WRITE(numout,*) ' '
+         WRITE(numout,*)
+         IF( ln_p5z    )  WRITE(numout,*) '   ==>>>   PISCES QUOTA model is used'
+         IF( ln_p4z    )  WRITE(numout,*) '   ==>>>   PISCES STANDARD model is used'
+         IF( ln_p2z    )  WRITE(numout,*) '   ==>>>   LOBSTER model is used'
+         IF( ln_ligand )  WRITE(numout,*) '   ==>>>   Compute remineralization/dissolution of organic ligands'
       ENDIF
     
