Index: branches/dev_001_GM/NEMO/TOP_SRC/CFC/trcini.cfc.h90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/CFC/trcini.cfc.h90	(revision 765)
+++ 	(revision )
@@ -1,124 +1,0 @@
-   !!======================================================================
-   !!                         ***  trcini.cfc.h90  ***
-   !! TOP :   Initialisation of CFC chemical model
-   !!======================================================================
-   !! History :    -   !  2004-06  (JC. Dutay) Original code
-   !!             1.0  !  2005-03  (O. Aumont, A. El Moussaoui) F90
-   !!              -   !  2005-10  (C. Ethe) Modularity
-   !!----------------------------------------------------------------------
-
-   CHARACTER (len=34) ::   clname = 'cfc1112.atm' ! ???
-
-   INTEGER  ::   inum                   ! unit number
-   REAL(wp) ::   ylats = -10.           ! 10 degrees south
-   REAL(wp) ::   ylatn =  10.           ! 10 degrees north
-
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_ini
-      !!----------------------------------------------------------------------
-      !!                     ***  trcini.cfc.h90  ***  
-      !!
-      !! ** Purpose : special initialization for cfc model
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jn, jl, jm
-      REAL(wp) ::   zyy,  zyd
-      !!----------------------------------------------------------------------
-
-      ! Initialization of boundaries conditions
-      ! --------------------------------------- 
-      pp_cfc(:,:,:) = 0.e0
-      qtr   (:,:,:) = 0.e0
-      xphem (:,:)   = 0.e0
-      DO jn = 1, jptra
-         DO jm = 1, jphem
-            DO jl = 1, jpyear
-               p_cfc(jl,jm,jn) = 0.0
-            END DO
-         END DO
-      END DO
-      
-      
-      ! Initialization of qint in case of  no restart 
-      !----------------------------------------------
-      IF( .NOT. lrsttr ) THEN    
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'Initialization de qint ; No restart : qint equal zero '
-         ENDIF
-         DO jn = 1, jptra
-            qint(:,:,jn) = 0.e0
-         END DO
-      ENDIF
-
-
-      !   READ CFC partial pressure atmospheric value :
-      !     p11(year,nt) = PCFC11  in northern (1) and southern (2) hemisphere 
-      !     p12(year,nt) = PCFC12  in northern (1) and southern (2) hemisphere 
-      !--------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*) 'read of formatted file cfc1112atm'
-      
-      CALL ctlopn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
-         &           1, numout, .FALSE., 1 )
-      REWIND(inum)
-      
-
-      ! Skip over 1st six descriptor lines
-      !-----------------------------------
-      DO jm = 1, 6
-         READ(inum,'(1x)')
-      END DO
-   
-   
-      !   Read file
-      ! ---------
-      DO jn = 31, 98
-         READ(inum,*) zyy, p_cfc(jn,1,jp11), p_cfc(jn,1,jp12), &
-            &              p_cfc(jn,2,jp11), p_cfc(jn,2,jp12)
-         WRITE(numout,'(f7.2, 4f8.2)' ) &
-            &         zyy, p_cfc(jn,1,jp11), p_cfc(jn,1,jp12), &
-            &              p_cfc(jn,2,jp11), p_cfc(jn,2,jp12)
-      END DO
-
-      p_cfc(32,1:2,jp11) = 5.e-4
-      p_cfc(33,1:2,jp11) = 8.e-4
-      p_cfc(34,1:2,jp11) = 1.e-6
-      p_cfc(35,1:2,jp11) = 2.e-3
-      p_cfc(36,1:2,jp11) = 4.e-3
-      p_cfc(37,1:2,jp11) = 6.e-3
-      p_cfc(38,1:2,jp11) = 8.e-3
-      p_cfc(39,1:2,jp11) = 1.e-2
-      
-      
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' Year   p11HN    p11HS    p12HN    p12HS '
-         DO jn = 30, 100
-            WRITE(numout, '( 1I4, 4F9.2)')   &
-               &         jn, p_cfc(jn,1,jp11), p_cfc(jn,2,jp11), &
-               &             p_cfc(jn,1,jp12), p_cfc(jn,2,jp12)
-         END DO
-      ENDIF
-
-
-      ! Interpolation factor of atmospheric partial pressure
-      ! Linear interpolation between 2 hemispheric function of latitud between ylats and ylatn
-      !---------------------------------------------------------------------------------------
-      zyd = ylatn - ylats      
-      DO jj = 1 , jpj
-         DO ji = 1 , jpi
-            IF(     gphit(ji,jj) >= ylatn ) THEN   ;   xphem(ji,jj) = 1.e0
-            ELSEIF( gphit(ji,jj) <= ylats ) THEN   ;   xphem(ji,jj) = 0.e0
-            ELSE                                   ;   xphem(ji,jj) = ( gphit(ji,jj) - ylats) / zyd
-            ENDIF
-         END DO
-      END DO
-      !
-   END SUBROUTINE trc_ini
Index: branches/dev_001_GM/NEMO/TOP_SRC/CFC/trcini_cfc.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/CFC/trcini_cfc.F90	(revision 766)
+++ branches/dev_001_GM/NEMO/TOP_SRC/CFC/trcini_cfc.F90	(revision 766)
@@ -0,0 +1,154 @@
+MODULE trcini_cfc
+   !!======================================================================
+   !!                         ***  MODULE trcini_cfc  ***
+   !! TOP :   initialisation of the CFC tracers
+   !!======================================================================
+   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) from trcini.cfc.h90
+   !!----------------------------------------------------------------------
+#if defined key_cfc
+   !!----------------------------------------------------------------------
+   !!   'key_cfc'                                               CFC tracers
+   !!----------------------------------------------------------------------
+   !! trc_ini_cfc      : CFC model initialisation
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE trccfc          ! CFC sms trends
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_ini_cfc   ! called by trcini.F90 module
+
+   CHARACTER (len=34) ::   clname = 'cfc1112.atm'   ! ???
+
+   INTEGER  ::   inum                   ! unit number
+   REAL(wp) ::   ylats = -10.           ! 10 degrees south
+   REAL(wp) ::   ylatn =  10.           ! 10 degrees north
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_ini_cfc
+      !!----------------------------------------------------------------------
+      !!                     ***  trc_ini_cfc  ***  
+      !!
+      !! ** Purpose :   initialization for cfc model
+      !!
+      !! ** Method  : - Read the namcfc namelist and check the parameter values
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jn, jl, jm
+      REAL(wp) ::   zyy  ,  zyd
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ini_cfc: initialisation of CFC chemical model'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
+
+      ! Initialization of boundaries conditions
+      ! --------------------------------------- 
+      pp_cfc(:,:,:) = 0.e0
+      qtr   (:,:,:) = 0.e0
+      xphem (:,:)   = 0.e0
+      DO jn = 1, jptra
+         DO jm = 1, jphem
+            DO jl = 1, jpyear
+               p_cfc(jl,jm,jn) = 0.0
+            END DO
+         END DO
+      END DO
+      
+      
+      ! Initialization of qint in case of  no restart 
+      !----------------------------------------------
+      IF( .NOT. lrsttr ) THEN    
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'Initialization de qint ; No restart : qint equal zero '
+         ENDIF
+         DO jn = 1, jptra
+            qint(:,:,jn) = 0.e0
+         END DO
+      ENDIF
+
+
+      !   READ CFC partial pressure atmospheric value :
+      !     p11(year,nt) = PCFC11  in northern (1) and southern (2) hemisphere 
+      !     p12(year,nt) = PCFC12  in northern (1) and southern (2) hemisphere 
+      !--------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*) 'read of formatted file cfc1112atm'
+      
+      CALL ctlopn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
+         &           1, numout, .FALSE., 1 )
+      REWIND(inum)
+      
+
+      ! Skip over 1st six descriptor lines
+      !-----------------------------------
+      DO jm = 1, 6
+         READ(inum,'(1x)')
+      END DO
+   
+   
+      !   Read file
+      ! ---------
+      DO jn = 31, 98
+         READ(inum,*) zyy, p_cfc(jn,1,jp11), p_cfc(jn,1,jp12), &
+            &              p_cfc(jn,2,jp11), p_cfc(jn,2,jp12)
+         WRITE(numout,'(f7.2, 4f8.2)' ) &
+            &         zyy, p_cfc(jn,1,jp11), p_cfc(jn,1,jp12), &
+            &              p_cfc(jn,2,jp11), p_cfc(jn,2,jp12)
+      END DO
+
+      p_cfc(32,1:2,jp11) = 5.e-4
+      p_cfc(33,1:2,jp11) = 8.e-4
+      p_cfc(34,1:2,jp11) = 1.e-6
+      p_cfc(35,1:2,jp11) = 2.e-3
+      p_cfc(36,1:2,jp11) = 4.e-3
+      p_cfc(37,1:2,jp11) = 6.e-3
+      p_cfc(38,1:2,jp11) = 8.e-3
+      p_cfc(39,1:2,jp11) = 1.e-2
+      
+      
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' Year   p11HN    p11HS    p12HN    p12HS '
+         DO jn = 30, 100
+            WRITE(numout, '( 1I4, 4F9.2)')   &
+               &         jn, p_cfc(jn,1,jp11), p_cfc(jn,2,jp11), &
+               &             p_cfc(jn,1,jp12), p_cfc(jn,2,jp12)
+         END DO
+      ENDIF
+
+
+      ! Interpolation factor of atmospheric partial pressure
+      ! Linear interpolation between 2 hemispheric function of latitud between ylats and ylatn
+      !---------------------------------------------------------------------------------------
+      zyd = ylatn - ylats      
+      DO jj = 1 , jpj
+         DO ji = 1 , jpi
+            IF(     gphit(ji,jj) >= ylatn ) THEN   ;   xphem(ji,jj) = 1.e0
+            ELSEIF( gphit(ji,jj) <= ylats ) THEN   ;   xphem(ji,jj) = 0.e0
+            ELSE                                   ;   xphem(ji,jj) = ( gphit(ji,jj) - ylats) / zyd
+            ENDIF
+         END DO
+      END DO
+      !
+   END SUBROUTINE trc_ini_cfc
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                                         No CFC tracers
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_ini_cfc             ! Empty routine
+   END SUBROUTINE trc_ini_cfc
+#endif
+
+   !!======================================================================
+END MODULE trcini_cfc
Index: branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcini.lobster1.h90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcini.lobster1.h90	(revision 765)
+++ 	(revision )
@@ -1,124 +1,0 @@
-   !!======================================================================
-   !!                         ***  trcini.lobster1.h90  ***
-   !! TOP :   Initialisation of LOBSTER 1 biological model
-   !!======================================================================
-   !! History :    -   !  1999-09  (M. Levy) Original code
-   !!              -   !  2000-12  (0. Aumont, E. Kestenare) add sediment 
-   !!             1.0  !  2004-03  (C. Ethe) Modularity
-   !!              -   !  2005-03  (O. Aumont, A. El Moussaoui) F90
-   !!----------------------------------------------------------------------
-
-#  include "domzgr_substitute.h90"
-#  include "passivetrc_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_ini
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE trc_ini  ***
-      !! ** purpose :   specific initialisation for lobster1 model
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zdm0(jpi,jpj,jpk), zrro(jpi,jpj), zfluo, zfluu
-      REAL(wp) ::   ztest, zfluo, zfluu
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrro
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdm0
-      !!----------------------------------------------------------------------
-
-      ! initialization of fields for optical model
-      ! --------------------------------------------
-      xze (:,:)   = 5.e0
-      xpar(:,:,:) = 0.e0
-
-      ! initialization for passive tracer remineralisation-damping  array
-      ! -----------------------------------------------------------------
-
-      DO jn = 1, jptra
-         remdmp(:,jn) = tminr
-      END DO
-
-      IF(lwp) THEN
-         WRITE(numout,*) ' '
-         WRITE(numout,*) ' trcini: compute remineralisation-damping  '
-         WRITE(numout,*) '         arrays for tracers'
-      ENDIF
-
-      ! initialization of biological variables
-      ! ------------------------------------------
-
-      ! Calculate vertical distribution of newly formed biogenic poc
-      ! in the water column in the case of max. possible bottom depth
-      ! ------------------------------------------------------------
-
-      zdm0   = 0.e0
-      zrro = 1.e0
-      DO jk = jpkb,jpkm1
-         DO jj =1, jpj
-            DO ji =1, jpi
-               zfluo = ( fsdepw(ji,jj,jk  ) / fsdepw(ji,jj,jpkb) )**xhr 
-               zfluu = ( fsdepw(ji,jj,jk+1) / fsdepw(ji,jj,jpkb) )**xhr
-               IF( zfluo.GT.1. )   zfluo = 1.e0
-               zdm0(ji,jj,jk) = zfluo - zfluu
-               IF( jk <= jpkb-1 )   zdm0(ji,jj,jk) = 0.e0
-               zrro(ji,jj) = zrro(ji,jj) - zdm0(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      zdm0(:,:,jpk) = zrro(:,:)
-
-      ! Calculate vertical distribution of newly formed biogenic poc
-      ! in the water column with realistic topography (first "dry" layer
-      ! contains total fraction, which has passed to the upper layers)
-      ! ----------------------------------------------------------------------
-      dminl = 0.
-      dmin3 = zdm0
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( tmask(ji,jj,jk) == 0. ) THEN
-                  dminl(ji,jj) = dminl(ji,jj) + dmin3(ji,jj,jk)
-                  dmin3(ji,jj,jk) = 0.e0
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( tmask(ji,jj,1) == 0 )   dmin3(ji,jj,1) = 0.e0
-         END DO
-      END DO
-
-      ! Coastal mask 
-      ! ------------   
-      cmask = 0.e0
-      DO ji = 2, jpi-1
-         DO jj = 2, jpj-1
-            if (tmask(ji,jj,1) == 1) then
-               ztest=tmask(ji+1,jj,1)*tmask(ji-1,jj,1)*tmask(ji,jj+1,1)*tmask(ji,jj-1,1)
-               IF (ztest == 0) cmask(ji,jj) = 1.
-            endif
-         END DO
-      END DO
-
-      cmask( 1 ,:) = cmask(jpi-1,:)
-      cmask(jpi,:) = cmask( 2   ,:)
-
-      !!gm BUG !!!!!   not valid in mpp and also not valid for north fold   !!!!!
-
-      ! Coastal surface
-      ! ---------------
-      areacot = 0.e0
-      DO ji = 2, jpi-1
-         DO jj = 2, jpj-1
-            areacot = areacot + e1t(ji,jj) * e2t(ji,jj) * cmask(ji,jj)
-         END DO
-      END DO
-      !
-   END SUBROUTINE trc_ini
Index: branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcini_lobster.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcini_lobster.F90	(revision 766)
+++ branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcini_lobster.F90	(revision 766)
@@ -0,0 +1,156 @@
+MODULE trcini_lobster
+   !!======================================================================
+   !!                         ***  MODULE trcini_lobster  ***
+   !! TOP :   initialisation of the LOBSTER biological model
+   !!======================================================================
+   !! History :    -   !  1999-09  (M. Levy) Original code
+   !!              -   !  2000-12  (0. Aumont, E. Kestenare) add sediment 
+   !!             1.0  !  2004-03  (C. Ethe) Modularity
+   !!              -   !  2005-03  (O. Aumont, A. El Moussaoui) F90
+   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) from trcini.lobster1.h90
+   !!----------------------------------------------------------------------
+#if defined key_trc_lobster1
+   !!----------------------------------------------------------------------
+   !!   'key_trc_lobster1'                                    LOBSTER model
+   !!----------------------------------------------------------------------
+   !! trc_ini_lobster  : LOBSTER model initialisation
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE trccfc          ! CFC sms trends
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_ini_lobster   ! called by trcini.F90 module
+
+#  include "domzgr_substitute.h90"
+#  include "passivetrc_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_ini_lobster
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE trc_ini_lobster  ***
+      !! ** purpose :   specific initialisation for LOBSTER bio-model
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jn
+      REAL(wp) ::   zdm0(jpi,jpj,jpk), zrro(jpi,jpj), zfluo, zfluu
+      REAL(wp) ::   ztest, zfluo, zfluu
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zrro
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdm0
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ini_lobster :   LOBSTER biochemical model initialisation'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
+
+
+      ! initialization of fields for optical model
+      ! --------------------------------------------
+      xze (:,:)   = 5.e0
+      xpar(:,:,:) = 0.e0
+
+      ! initialization for passive tracer remineralisation-damping  array
+      ! -----------------------------------------------------------------
+
+      DO jn = 1, jptra
+         remdmp(:,jn) = tminr
+      END DO
+
+      IF(lwp) THEN
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' trcini: compute remineralisation-damping  '
+         WRITE(numout,*) '         arrays for tracers'
+      ENDIF
+
+      ! initialization of biological variables
+      ! ------------------------------------------
+
+      ! Calculate vertical distribution of newly formed biogenic poc
+      ! in the water column in the case of max. possible bottom depth
+      ! ------------------------------------------------------------
+
+      zdm0   = 0.e0
+      zrro = 1.e0
+      DO jk = jpkb,jpkm1
+         DO jj =1, jpj
+            DO ji =1, jpi
+               zfluo = ( fsdepw(ji,jj,jk  ) / fsdepw(ji,jj,jpkb) )**xhr 
+               zfluu = ( fsdepw(ji,jj,jk+1) / fsdepw(ji,jj,jpkb) )**xhr
+               IF( zfluo.GT.1. )   zfluo = 1.e0
+               zdm0(ji,jj,jk) = zfluo - zfluu
+               IF( jk <= jpkb-1 )   zdm0(ji,jj,jk) = 0.e0
+               zrro(ji,jj) = zrro(ji,jj) - zdm0(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      zdm0(:,:,jpk) = zrro(:,:)
+
+      ! Calculate vertical distribution of newly formed biogenic poc
+      ! in the water column with realistic topography (first "dry" layer
+      ! contains total fraction, which has passed to the upper layers)
+      ! ----------------------------------------------------------------------
+      dminl = 0.
+      dmin3 = zdm0
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj,jk) == 0. ) THEN
+                  dminl(ji,jj) = dminl(ji,jj) + dmin3(ji,jj,jk)
+                  dmin3(ji,jj,jk) = 0.e0
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( tmask(ji,jj,1) == 0 )   dmin3(ji,jj,1) = 0.e0
+         END DO
+      END DO
+
+      ! Coastal mask 
+      ! ------------   
+      cmask = 0.e0
+      DO ji = 2, jpi-1
+         DO jj = 2, jpj-1
+            if (tmask(ji,jj,1) == 1) then
+               ztest=tmask(ji+1,jj,1)*tmask(ji-1,jj,1)*tmask(ji,jj+1,1)*tmask(ji,jj-1,1)
+               IF (ztest == 0) cmask(ji,jj) = 1.
+            endif
+         END DO
+      END DO
+
+      cmask( 1 ,:) = cmask(jpi-1,:)
+      cmask(jpi,:) = cmask( 2   ,:)
+
+      !!gm BUG !!!!!   not valid in mpp and also not valid for north fold   !!!!!
+
+      ! Coastal surface
+      ! ---------------
+      areacot = 0.e0
+      DO ji = 2, jpi-1
+         DO jj = 2, jpj-1
+            areacot = areacot + e1t(ji,jj) * e2t(ji,jj) * cmask(ji,jj)
+         END DO
+      END DO
+      !
+   END SUBROUTINE trc_ini_lobster
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                                   No LOBSTER bio-model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_ini_lobster             ! Empty routine
+   END SUBROUTINE trc_ini_lobster
+#endif
+
+   !!======================================================================
+END MODULE trcini_lobster
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini.pisces.h90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini.pisces.h90	(revision 765)
+++ 	(revision )
@@ -1,544 +1,0 @@
-   !!======================================================================
-   !!                         ***  trcini.pisces.h90  ***
-   !! TOP :   Initialisation of PISCES biological model
-   !!======================================================================
-   !! History :    -   !  1988-07  (E. Maier-Reiner) Original code
-   !!              -   !  1999-10  (O. Aumont, C. Le Quere)
-   !!              -   !  2002     (O. Aumont)  PISCES
-   !!             1.0  !  2005-03  (O. Aumont, A. El Moussaoui) F90
-   !!----------------------------------------------------------------------
-
-#  include "domzgr_substitute.h90"
-#  include "passivetrc_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trc_ini
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE trc_ini ***
-      !!
-      !! ** Purpose :   Initialisation of PISCES biological and chemical variables
-      !!----------------------------------------------------------------------
-      USE iom
-      !!
-      INTEGER :: ji,jj,jk
-      INTEGER :: ichl,iband,jm
-      INTEGER , PARAMETER :: jpmois = 12, jpan   = 1 
-
-      REAL(wp) :: zcoef
-      REAL(wp) :: ztoto,expide,denitide,ztra,zmaskt
-      REAL(wp) , DIMENSION (jpi,jpj) :: riverdoc,river,ndepo
-      REAL(wp) , DIMENSION (jpi,jpj,jpk) :: cmask
-
-      INTEGER :: numriv, numdust, numbath, numdep
-      INTEGER :: numlight 
-
-#if defined key_trc_kriest
-      INTEGER  ::   jn, kiter
-      REAL(wp) ::   znum, zdiv
-      REAL(wp) ::   zws,zwr, zwl,wmax, xnummax, &
-      REAL(wp) ::   zmin, zmax, zl, zr, xacc
-#endif
-      !!----------------------------------------------------------------------
-
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' trc_ini :   PISCES biological and chemical initialisation'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~'
-
-
-      !                                            ! Time-step
-      rfact   = rdttra(1) * float(ndttrc)          ! ---------
-      rfactr  = 1. / rfact
-      rfact2  = rfact / float(nrdttrc)
-      rfact2r = 1. / rfact2
-
-      IF(lwp) WRITE(numout,*) '    Tracer  time step=', rfact, ' rdt = ', rdt
-      IF(lwp) write(numout,*) '    Biology time step=', rfact2
-
-
-      !                                            ! Dust input from the atmosphere
-      IF( bdustfer ) THEN                          ! ------------------------------
-         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
-         CALL iom_open ( 'dust.orca.nc', numdust )
-         DO jm = 1, jpmois
-            CALL iom_get( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm )
-         END DO
-         CALL iom_close( numdust )
-      ELSE
-         dustmo(:,:,:) = 0.e0
-      ENDIF
-
-
-      !                                            ! Nutrient input from rivers
-      IF( briver ) THEN                            ! --------------------------
-         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
-         CALL iom_open ( 'river.orca.nc', numriv )
-         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
-         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
-         CALL iom_close( numriv )
-      ELSE
-         river   (:,:) = 0.e0
-         riverdoc(:,:) = 0.e0
-      endif
-
-      !                                            ! Nutrient input from dust
-      IF( bndepo ) THEN                            ! ------------------------
-         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
-         CALL iom_open ( 'ndeposition.orca.nc', numdep )
-         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
-         CALL iom_close( numdep )
-      ELSE
-         ndepo(:,:) = 0.e0
-      ENDIF
-
-      !                                            ! Coastal and island masks
-      IF( bsedinput ) THEN                         ! ------------------------
-         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron '
-         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
-         CALL iom_open ( 'bathy.orca.nc', numbath )
-         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
-         CALL iom_close( numbath )
-         !
-         DO jk = 1, 5
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF( tmask(ji,jj,jk) /= 0. ) THEN
-                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
-                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
-                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
-                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
-                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
-               END DO
-            END DO
-         END DO    
-      ELSE
-         cmask(:,:,:) = 0.e0
-      ENDIF
-
-      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
-
-
-      !                                            ! total atmospheric supply of Si
-      !                                            ! ------------------------------
-      sumdepsi = 0.e0
-      DO jm = 1, jpmois
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               sumdepsi = sumdepsi + dustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
-                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
-
-      !                                            ! N/P and Si releases due to coastal rivers
-      !                                            ! -----------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
-            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
-            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
-            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
-      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
-
-      rivpo4input=0.e0
-      rivalkinput=0.e0
-      nitdepinput=0.e0
-      DO jj = 2 , jpjm1
-         DO ji = 2, jpim1
-            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
-            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
-            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
-            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
-         END DO
-      END DO
-      IF( lk_mpp ) THEN
-         CALL mpp_sum( rivpo4input )  ! sum over the global domain
-         CALL mpp_sum( rivalkinput )  ! sum over the global domain
-         CALL mpp_sum( nitdepinput )  ! sum over the global domain
-      ENDIF
-
-
-      !                                            !    Coastal supply of iron
-      !                                            ! -------------------------
-      DO jk = 1, jpkm1
-         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
-      END DO
-      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
-
-
-
-#if defined key_trc_kriest
-      !----------------------------------------------------------------------
-      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
-      !  Search of the maximum number of particles in aggregates for each k-level.  
-      !  Bissection Method
-      !--------------------------------------------------------------------
-      WRITE(numout,*)
-      WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
-            
-      xacc     =  0.001
-      kiter    = 50
-      zmin     =  1.10
-      zmax     = xkr_mass_max / xkr_mass_min
-      xkr_frac = zmax
-
-      DO jk =1,jpk
-         zl = zmin
-         zr = zmax
-         wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
-         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
-         znum = zl - 1.
-         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
-            & - ( xkr_wsbio_max * xkr_eta * znum * &
-            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
-            & - wmax
-  
-         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
-         znum = zr - 1.
-         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
-            & - ( xkr_wsbio_max * xkr_eta * znum * &
-            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
-            & - wmax
-
-iflag:  DO jn = 1, kiter               
-           IF( zwl == 0.e0 ) THEN
-              xnummax = zl
-           ELSE IF ( zwr == 0.e0 ) THEN
-              xnummax = zr
-           ELSE
-              xnummax = ( zr + zl ) / 2.
-              zdiv = xkr_zeta + xkr_eta - xkr_eta * xnummax
-              znum = xnummax - 1.
-              zws =  xkr_wsbio_min * xkr_zeta / zdiv &
-                 & - ( xkr_wsbio_max * xkr_eta * znum * &
-                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
-                 & - wmax
-              IF( zws * zwl < 0. ) THEN
-                 zr = xnummax
-              ELSE
-                 zl = xnummax
-              ENDIF
-              zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
-              znum = zl - 1.
-              zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
-                 & - ( xkr_wsbio_max * xkr_eta * znum * &
-                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
-                 & - wmax
-              
-              zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
-              znum = zr - 1.
-              zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
-                 & - ( xkr_wsbio_max * xkr_eta * znum * &
-                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
-                 & - wmax
-
-              IF ( ABS ( zws )  <= xacc ) EXIT iflag
-
-           ENDIF
-                    
-        END DO iflag
-               
-        xnumm(jk) = xnummax
-        WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
-        
-     END DO
-
-#endif
-
-      !----------------------------------------------------------------------
-      ! Initialize biological variables 
-      !----------------------------------------------------------------------
-      ! Set biological ratios
-      ! ---------------------
-      rno3   = (16.+2.) / 122.
-      po4r   =   1.e0   / 122.
-      o2nit  =  32.     / 122.
-      rdenit =  97.6    /  16.
-      o2ut   = 140.     / 122.
-
-      !----------------------------------------------------------------------
-      ! Initialize chemical variables 
-      !----------------------------------------------------------------------
-
-      ! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
-      ! ----------------------------------------------------------
-      atcox = 0.20946
-
-      ! Set lower/upper limits for temperature and salinity
-      ! ---------------------------------------------------
-      salchl = 1.e0 / 1.80655
-      calcon = 1.03e-2
-
-      ! Set coefficients for apparent solubility equilibrium of calcite
-      ! Millero et al. 1995 from Mucci 1983
-      ! --------------------------------------------------------------
-      akcc1 = -171.9065
-      akcc2 =   -0.077993
-      akcc3 = 2839.319
-      akcc4 =   71.595
-      akcc5 =   -0.77712
-      akcc6 =    0.0028426
-      akcc7 =  178.34
-      akcc8 =   -0.07711
-      akcc9 =    0.0041249
-
-      ! Set coefficients for seawater pressure correction
-      ! -------------------------------------------------
-      devk1(1) = -25.5
-      devk2(1) =   0.1271
-      devk3(1) =   0.e0
-      devk4(1) =  -3.08E-3
-      devk5(1) =   0.0877E-3
-      !
-      devk1(2) = -15.82
-      devk2(2) =  -0.0219
-      devk3(2) =   0.e0
-      devk4(2) =   1.13E-3
-      devk5(2) =  -0.1475E-3
-      !
-      devk1(3) = -29.48
-      devk2(3) =   0.1622
-      devk3(3) =   2.608E-3
-      devk4(3) =  -2.84E-3
-      devk5(3) =   0.e0
-      !
-      devk1(4) = -14.51
-      devk2(4) =   0.1211
-      devk3(4) =  -0.321E-3
-      devk4(4) =  -2.67E-3
-      devk5(4) =   0.0427E-3
-      !
-      devk1(5) = -23.12
-      devk2(5) =   0.1758
-      devk3(5) =  -2.647E-3
-      devk4(5) =  -5.15E-3
-      devk5(5) =   0.09E-3
-      !
-      devk1(6) = -26.57
-      devk2(6) =   0.2020
-      devk3(6) =  -3.042E-3
-      devk4(6) =  -4.08E-3
-      devk5(6) =   0.0714E-3
-      !
-      devk1(7) = -25.60
-      devk2(7) =   0.2324
-      devk3(7) =  -3.6246E-3
-      devk4(7) =  -5.13E-3
-      devk5(7) =   0.0794E-3
-      !
-      ! For calcite with Edmond and Gieske 1970
-      !     devkst = 0.23
-      !     devks  = 35.4
-      ! Millero 95 takes this depth dependance for calcite
-      devk1(8) = -48.76
-      devk2(8) =   0.5304
-      devk3(8) =   0.e0
-      devk4(8) = -11.76E-3
-      devk5(8) =   0.3692E-3
-      !
-      ! Coefficients for sulfate and fluoride
-      devk1(9) = -18.03
-      devk2(9) =   0.0466
-      devk3(9) =   0.316e-3
-      devk4(9) =  -4.53e-3
-      devk5(9) =   0.09e-3
-
-      devk1(10) = -9.78
-      devk2(10) = -0.0090
-      devk3(10) = -0.942e-3
-      devk4(10) = -3.91e-3
-      devk5(10) =  0.054e-3
-
-
-      ! Set universal gas constants
-      ! ---------------------------
-      rgas  = 83.143
-      oxyco =  1.e0 / 22.4144
-
-      ! Set boron constants
-      ! -------------------
-      bor1 = 0.00023
-      bor2 = 1.e0 / 10.82
-
-      ! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
-      ! -----------------------------------------------------------------
-      c00 = -60.2409
-      c01 =  93.4517
-      c02 =  23.3585
-      c03 =   0.023517
-      c04 =  -0.023656
-      c05 =   0.0047036
-      !
-      ca0 = -162.8301
-      ca1 =  218.2968
-      ca2 =   90.9241
-      ca3 =   -1.47696
-      ca4 =    0.025695
-      ca5 =   -0.025225
-      ca6 =    0.0049867
-
-      ! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
-      ! ---------------------------------------------------------------------
-      c10 = -3670.7
-      c11 =    62.008
-      c12 =    -9.7944
-      c13 =     0.0118
-      c14 =    -0.000116
-
-      ! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
-      ! ---------------------------------------------------------------------
-      c20 = -1394.7
-      c21 =    -4.777
-      c22 =     0.0184
-      c23 =    -0.000118
-
-      ! Set constants for calculate concentrations for sulfate and fluoride
-      ! sulfates (Morris & Riley 1966)
-      !----------------------------------------------------------------------
-      st1 = 0.14
-      st2 = 1.e0 / 96.062
-
-      ! fluoride
-      ! --------
-      ft1 = 0.000067
-      ft2 = 1.e0 / 18.9984
-
-      ! sulfates (Dickson 1990 change to mol:kg soln, idem OCMIP)
-      !----------------------------------------------------------
-      ks0  =    141.328
-      ks1  =  -4276.1
-      ks2  =    -23.093
-      ks3  = -13856.
-      ks4  =    324.57
-      ks5  =    -47.986
-      ks6  =  35474.
-      ks7  =   -771.54
-      ks8  =    114.723
-      ks9  =  -2698.
-      ks10 =   1776.
-      ks11 =      1.
-      ks12 =     -0.001005
-
-      ! fluorides (Dickson & Riley 1979 change to mol/kg soln)
-      !-------------------------------------------------------
-      kf0 =  -12.641
-      kf1 = 1590.2
-      kf2 =    1.525
-      kf3 =    1.0
-      kf4 =   -0.001005
-
-      ! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
-      ! ------------------------------------------------------------------
-      cb0  = -8966.90
-      cb1  = -2890.53
-      cb2  =   -77.942
-      cb3  =     1.728
-      cb4  =    -0.0996
-      cb5  =   148.0248
-      cb6  =   137.1942
-      cb7  =     1.62142
-      cb8  =   -24.4344
-      cb9  =   -25.085
-      cb10 =    -0.2474
-      cb11 =     0.053105
-
-      ! Set coeff. for dissoc. of water (Dickson and Riley, 1979, 
-      !   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415 
-      !   after pers. commun. to B. Bacastow, 1988)
-      ! ---------------------------------------------------------
-      cw0 = -13847.26
-      cw1 =    148.9652
-      cw2 =    -23.6521
-      cw3 =    118.67
-      cw4 =     -5.977
-      cw5 =      1.0495
-      cw6 =     -0.01615
-
-      ! Set coeff. for dissoc. of phosphate (Millero (1974)
-      ! ---------------------------------------------------
-      cp10 =   115.54
-      cp11 = -4576.752
-      cp12 =   -18.453
-      cp13 =  -106.736
-      cp14 =     0.69171
-      cp15 =    -0.65643
-      cp16 =    -0.01844
-      !
-      cp20 =   172.1033
-      cp21 = -8814.715
-      cp22 =   -27.927
-      cp23 =  -160.340
-      cp24 =     1.3566
-      cp25 =     0.37335
-      cp26 =    -0.05778
-      !
-      cp30 =   -18.126
-      cp31 = -3070.75
-      cp32 =    17.27039
-      cp33 =     2.81197
-      cp34 =   -44.99486
-      cp35 =    -0.09984
-
-      ! Set coeff. for dissoc. of phosphate (Millero (1974)
-      ! ---------------------------------------------------
-      cs10 =   117.385
-      cs11 = -8904.2
-      cs12 =   -19.334
-      cs13 =  -458.79
-      cs14 =     3.5913
-      cs15 =   188.74
-      cs16 =    -1.5998
-      cs17 =   -12.1652
-      cs18 =     0.07871
-      cs19 =     0.e0
-      cs20 =     1.e0
-      cs21 =    -0.001005
-
-
-      ! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
-      ! ----------------------------------------------------------------
-      ox0 = -58.3877
-      ox1 =  85.8079
-      ox2 =  23.8439
-      ox3 =  -0.034892
-      ox4 =   0.015568
-      ox5 =  -0.0019387
-
-      !  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
-      !  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
-
-      CALL ctlopn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
-         &           1, numout, .TRUE., 1 )
-      DO ichl = 1,61
-         READ(numlight,*) ztoto, ( xkrgb(iband,ichl), iband = 1,3 )
-      END DO
-      CLOSE(numlight)
-
-
-      CALL p4zche        ! initialize the chemical constants
-
-
-      ndayflxtr = 0      !  Initialize a counter for the computation of chemistry
-
-
-      IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done'
-      !
-   END SUBROUTINE trc_ini
Index: branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 766)
+++ branches/dev_001_GM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 766)
@@ -0,0 +1,571 @@
+MODULE trcini_pisces
+   !!======================================================================
+   !!                         ***  MODULE trcini_pisces  ***
+   !! TOP :   initialisation of the PISCES biochemical model
+   !!======================================================================
+   !! History :    -   !  1988-07  (E. Maier-Reiner) Original code
+   !!              -   !  1999-10  (O. Aumont, C. Le Quere)
+   !!              -   !  2002     (O. Aumont)  PISCES
+   !!             1.0  !  2005-03  (O. Aumont, A. El Moussaoui) F90
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trcini.pisces.h90
+   !!----------------------------------------------------------------------
+#if defined key_trc_pisces  &&  defined key_trc_kriest
+   !!---------------------------------------------------------------------
+   !!   'key_trc_pisces' & 'key_trc_kriest'         PISCES bio-model + ???
+   !!---------------------------------------------------------------------
+   !! trc_ini_pisces   : PISCES biochemical model initialisation
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE trccfc          ! CFC sms trends
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_ini_pisces   ! called by trcini.F90 module
+
+#  include "domzgr_substitute.h90"
+#  include "passivetrc_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_ini_pisces
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE trc_ini_pisces ***
+      !!
+      !! ** Purpose :   Initialisation of the PISCES biochemical model
+      !!----------------------------------------------------------------------
+      INTEGER :: ji, jj, jk, jm
+      INTEGER :: ichl, iband
+      INTEGER , PARAMETER ::   jpmois = 12, jpan   = 1 
+
+      REAL(wp) :: zcoef
+      REAL(wp) :: ztoto,expide,denitide,ztra,zmaskt
+      REAL(wp) , DIMENSION (jpi,jpj) :: riverdoc,river,ndepo
+      REAL(wp) , DIMENSION (jpi,jpj,jpk) :: cmask
+
+      INTEGER :: numriv, numdust, numbath, numdep
+      INTEGER :: numlight 
+
+#if defined key_trc_kriest
+      INTEGER  ::   jn, kiter
+      REAL(wp) ::   znum, zdiv
+      REAL(wp) ::   zws,zwr, zwl,wmax, xnummax, &
+      REAL(wp) ::   zmin, zmax, zl, zr, xacc
+#endif
+      !!----------------------------------------------------------------------
+
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ini_pisces :   PISCES biochemical model initialisation'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+
+      !                                            ! Time-step
+      rfact   = rdttra(1) * float(ndttrc)          ! ---------
+      rfactr  = 1. / rfact
+      rfact2  = rfact / float(nrdttrc)
+      rfact2r = 1. / rfact2
+
+      IF(lwp) WRITE(numout,*) '    Tracer  time step    rfact  = ', rfact, ' rdt = ', rdt
+      IF(lwp) write(numout,*) '    Biology time step    rfact2 = ', rfact2
+
+
+      !                                            ! Dust input from the atmosphere
+      IF( bdustfer ) THEN                          ! ------------------------------
+         IF(lwp) WRITE(numout,*) '    Initialize dust input from atmosphere '
+         CALL iom_open ( 'dust.orca.nc', numdust )
+         DO jm = 1, jpmois
+            CALL iom_get( numdust, jpdom_data, 'dust', dustmo(:,:,jm), jm )
+         END DO
+         CALL iom_close( numdust )
+      ELSE
+         dustmo(:,:,:) = 0.e0
+      ENDIF
+
+
+      !                                            ! Nutrient input from rivers
+      IF( briver ) THEN                            ! --------------------------
+         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by rivers from river.orca.nc file'
+         CALL iom_open ( 'river.orca.nc', numriv )
+         CALL iom_get  ( numriv, jpdom_data, 'riverdic', river   (:,:), jpan )
+         CALL iom_get  ( numriv, jpdom_data, 'riverdoc', riverdoc(:,:), jpan )
+         CALL iom_close( numriv )
+      ELSE
+         river   (:,:) = 0.e0
+         riverdoc(:,:) = 0.e0
+      endif
+
+      !                                            ! Nutrient input from dust
+      IF( bndepo ) THEN                            ! ------------------------
+         IF(lwp) WRITE(numout,*) '    Initialize the nutrient input by dust from ndeposition.orca.nc'
+         CALL iom_open ( 'ndeposition.orca.nc', numdep )
+         CALL iom_get  ( numdep, jpdom_data, 'ndep', ndepo(:,:), jpan )
+         CALL iom_close( numdep )
+      ELSE
+         ndepo(:,:) = 0.e0
+      ENDIF
+
+      !                                            ! Coastal and island masks
+      IF( bsedinput ) THEN                         ! ------------------------
+         IF(lwp) WRITE(numout,*) '    Computation of an island mask to enhance coastal supply of iron '
+         IF(lwp) WRITE(numout,*) '       from bathy.orca.nc file '
+         CALL iom_open ( 'bathy.orca.nc', numbath )
+         CALL iom_get  ( numbath, jpdom_data, 'bathy', cmask(:,:,:), jpan )
+         CALL iom_close( numbath )
+         !
+         DO jk = 1, 5
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF( tmask(ji,jj,jk) /= 0. ) THEN
+                     zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk)    &
+                        &                       * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
+                     IF( zmaskt == 0. )   cmask(ji,jj,jk ) = 0.1
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  expide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
+                  denitide = -0.9543 + 0.7662 * LOG( expide ) - 0.235 * LOG( expide )**2
+                  cmask(ji,jj,jk) = cmask(ji,jj,jk) * MIN( 1., EXP( denitide ) / 0.5 )
+               END DO
+            END DO
+         END DO    
+      ELSE
+         cmask(:,:,:) = 0.e0
+      ENDIF
+
+      CALL lbc_lnk( cmask , 'T', 1. )      ! Lateral boundary conditions on cmask   (sign unchanged)
+
+
+      !                                            ! total atmospheric supply of Si
+      !                                            ! ------------------------------
+      sumdepsi = 0.e0
+      DO jm = 1, jpmois
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               sumdepsi = sumdepsi + dustmo(ji,jj,jm) / (12.*rmoss) * 8.8        &
+                  &     * 0.075/28.1 * e1t(ji,jj) * e2t(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
+            END DO
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( sumdepsi )  ! sum over the global domain
+
+      !                                            ! N/P and Si releases due to coastal rivers
+      !                                            ! -----------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zcoef = raass * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
+            cotdep(ji,jj) =  river(ji,jj)                  *1E9 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
+            rivinp(ji,jj) = (river(ji,jj)+riverdoc(ji,jj)) *1E9 / ( 31.6* zcoef + rtrn ) * tmask(ji,jj,1)
+            nitdep(ji,jj) = 7.6 * ndepo(ji,jj)                  / ( 14E6*raass*fse3t(ji,jj,1) + rtrn ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      ! Lateral boundary conditions on ( cotdep, rivinp, nitdep )   (sign unchanged)
+      CALL lbc_lnk( cotdep , 'T', 1. )  ;  CALL lbc_lnk( rivinp , 'T', 1. )  ;  CALL lbc_lnk( nitdep , 'T', 1. )
+
+      rivpo4input=0.e0
+      rivalkinput=0.e0
+      nitdepinput=0.e0
+      DO jj = 2 , jpjm1
+         DO ji = 2, jpim1
+            zcoef = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)) * tmask(ji,jj,1) * tmask_i(ji,jj) * raass
+            rivpo4input = rivpo4input + rivinp(ji,jj) * zcoef
+            rivalkinput = rivalkinput + cotdep(ji,jj) * zcoef
+            nitdepinput = nitdepinput + nitdep(ji,jj) * zcoef
+         END DO
+      END DO
+      IF( lk_mpp ) THEN
+         CALL mpp_sum( rivpo4input )  ! sum over the global domain
+         CALL mpp_sum( rivalkinput )  ! sum over the global domain
+         CALL mpp_sum( nitdepinput )  ! sum over the global domain
+      ENDIF
+
+
+      !                                            !    Coastal supply of iron
+      !                                            ! -------------------------
+      DO jk = 1, jpkm1
+         ironsed(:,:,jk) = sedfeinput * cmask(:,:,jk) / ( fse3t(:,:,jk) * rjjss )
+      END DO
+      CALL lbc_lnk( ironsed , 'T', 1. )      ! Lateral boundary conditions on ( ironsed )   (sign unchanged)
+
+
+
+#if defined key_trc_kriest
+      !----------------------------------------------------------------------
+      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
+      !  Search of the maximum number of particles in aggregates for each k-level.  
+      !  Bissection Method
+      !--------------------------------------------------------------------
+      WRITE(numout,*)
+      WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
+            
+      xacc     =  0.001
+      kiter    = 50
+      zmin     =  1.10
+      zmax     = xkr_mass_max / xkr_mass_min
+      xkr_frac = zmax
+
+      DO jk =1,jpk
+         zl = zmin
+         zr = zmax
+         wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
+         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
+         znum = zl - 1.
+         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
+            & - ( xkr_wsbio_max * xkr_eta * znum * &
+            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+            & - wmax
+  
+         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
+         znum = zr - 1.
+         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
+            & - ( xkr_wsbio_max * xkr_eta * znum * &
+            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+            & - wmax
+
+iflag:  DO jn = 1, kiter               
+           IF( zwl == 0.e0 ) THEN
+              xnummax = zl
+           ELSE IF ( zwr == 0.e0 ) THEN
+              xnummax = zr
+           ELSE
+              xnummax = ( zr + zl ) / 2.
+              zdiv = xkr_zeta + xkr_eta - xkr_eta * xnummax
+              znum = xnummax - 1.
+              zws =  xkr_wsbio_min * xkr_zeta / zdiv &
+                 & - ( xkr_wsbio_max * xkr_eta * znum * &
+                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+                 & - wmax
+              IF( zws * zwl < 0. ) THEN
+                 zr = xnummax
+              ELSE
+                 zl = xnummax
+              ENDIF
+              zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
+              znum = zl - 1.
+              zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
+                 & - ( xkr_wsbio_max * xkr_eta * znum * &
+                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+                 & - wmax
+              
+              zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
+              znum = zr - 1.
+              zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
+                 & - ( xkr_wsbio_max * xkr_eta * znum * &
+                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
+                 & - wmax
+
+              IF ( ABS ( zws )  <= xacc ) EXIT iflag
+
+           ENDIF
+                    
+        END DO iflag
+               
+        xnumm(jk) = xnummax
+        WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
+        
+     END DO
+
+#endif
+
+      !----------------------------------------------------------------------
+      ! Initialize biological variables 
+      !----------------------------------------------------------------------
+      ! Set biological ratios
+      ! ---------------------
+      rno3   = (16.+2.) / 122.
+      po4r   =   1.e0   / 122.
+      o2nit  =  32.     / 122.
+      rdenit =  97.6    /  16.
+      o2ut   = 140.     / 122.
+
+      !----------------------------------------------------------------------
+      ! Initialize chemical variables 
+      !----------------------------------------------------------------------
+
+      ! set pre-industrial atmospheric [co2] (ppm) and o2/n2 ratio
+      ! ----------------------------------------------------------
+      atcox = 0.20946
+
+      ! Set lower/upper limits for temperature and salinity
+      ! ---------------------------------------------------
+      salchl = 1.e0 / 1.80655
+      calcon = 1.03e-2
+
+      ! Set coefficients for apparent solubility equilibrium of calcite
+      ! Millero et al. 1995 from Mucci 1983
+      ! --------------------------------------------------------------
+      akcc1 = -171.9065
+      akcc2 =   -0.077993
+      akcc3 = 2839.319
+      akcc4 =   71.595
+      akcc5 =   -0.77712
+      akcc6 =    0.0028426
+      akcc7 =  178.34
+      akcc8 =   -0.07711
+      akcc9 =    0.0041249
+
+      ! Set coefficients for seawater pressure correction
+      ! -------------------------------------------------
+      devk1(1) = -25.5
+      devk2(1) =   0.1271
+      devk3(1) =   0.e0
+      devk4(1) =  -3.08E-3
+      devk5(1) =   0.0877E-3
+      !
+      devk1(2) = -15.82
+      devk2(2) =  -0.0219
+      devk3(2) =   0.e0
+      devk4(2) =   1.13E-3
+      devk5(2) =  -0.1475E-3
+      !
+      devk1(3) = -29.48
+      devk2(3) =   0.1622
+      devk3(3) =   2.608E-3
+      devk4(3) =  -2.84E-3
+      devk5(3) =   0.e0
+      !
+      devk1(4) = -14.51
+      devk2(4) =   0.1211
+      devk3(4) =  -0.321E-3
+      devk4(4) =  -2.67E-3
+      devk5(4) =   0.0427E-3
+      !
+      devk1(5) = -23.12
+      devk2(5) =   0.1758
+      devk3(5) =  -2.647E-3
+      devk4(5) =  -5.15E-3
+      devk5(5) =   0.09E-3
+      !
+      devk1(6) = -26.57
+      devk2(6) =   0.2020
+      devk3(6) =  -3.042E-3
+      devk4(6) =  -4.08E-3
+      devk5(6) =   0.0714E-3
+      !
+      devk1(7) = -25.60
+      devk2(7) =   0.2324
+      devk3(7) =  -3.6246E-3
+      devk4(7) =  -5.13E-3
+      devk5(7) =   0.0794E-3
+      !
+      ! For calcite with Edmond and Gieske 1970
+      !     devkst = 0.23
+      !     devks  = 35.4
+      ! Millero 95 takes this depth dependance for calcite
+      devk1(8) = -48.76
+      devk2(8) =   0.5304
+      devk3(8) =   0.e0
+      devk4(8) = -11.76E-3
+      devk5(8) =   0.3692E-3
+      !
+      ! Coefficients for sulfate and fluoride
+      devk1(9) = -18.03
+      devk2(9) =   0.0466
+      devk3(9) =   0.316e-3
+      devk4(9) =  -4.53e-3
+      devk5(9) =   0.09e-3
+
+      devk1(10) = -9.78
+      devk2(10) = -0.0090
+      devk3(10) = -0.942e-3
+      devk4(10) = -3.91e-3
+      devk5(10) =  0.054e-3
+
+
+      ! Set universal gas constants
+      ! ---------------------------
+      rgas  = 83.143
+      oxyco =  1.e0 / 22.4144
+
+      ! Set boron constants
+      ! -------------------
+      bor1 = 0.00023
+      bor2 = 1.e0 / 10.82
+
+      ! Set volumetric solubility constants for co2 in ml/l (Weiss, 1974)
+      ! -----------------------------------------------------------------
+      c00 = -60.2409
+      c01 =  93.4517
+      c02 =  23.3585
+      c03 =   0.023517
+      c04 =  -0.023656
+      c05 =   0.0047036
+      !
+      ca0 = -162.8301
+      ca1 =  218.2968
+      ca2 =   90.9241
+      ca3 =   -1.47696
+      ca4 =    0.025695
+      ca5 =   -0.025225
+      ca6 =    0.0049867
+
+      ! Set coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
+      ! ---------------------------------------------------------------------
+      c10 = -3670.7
+      c11 =    62.008
+      c12 =    -9.7944
+      c13 =     0.0118
+      c14 =    -0.000116
+
+      ! Set coeff. for 2. dissoc. of carbonic acid (Edmond and Gieskes, 1970)
+      ! ---------------------------------------------------------------------
+      c20 = -1394.7
+      c21 =    -4.777
+      c22 =     0.0184
+      c23 =    -0.000118
+
+      ! Set constants for calculate concentrations for sulfate and fluoride
+      ! sulfates (Morris & Riley 1966)
+      !----------------------------------------------------------------------
+      st1 = 0.14
+      st2 = 1.e0 / 96.062
+
+      ! fluoride
+      ! --------
+      ft1 = 0.000067
+      ft2 = 1.e0 / 18.9984
+
+      ! sulfates (Dickson 1990 change to mol:kg soln, idem OCMIP)
+      !----------------------------------------------------------
+      ks0  =    141.328
+      ks1  =  -4276.1
+      ks2  =    -23.093
+      ks3  = -13856.
+      ks4  =    324.57
+      ks5  =    -47.986
+      ks6  =  35474.
+      ks7  =   -771.54
+      ks8  =    114.723
+      ks9  =  -2698.
+      ks10 =   1776.
+      ks11 =      1.
+      ks12 =     -0.001005
+
+      ! fluorides (Dickson & Riley 1979 change to mol/kg soln)
+      !-------------------------------------------------------
+      kf0 =  -12.641
+      kf1 = 1590.2
+      kf2 =    1.525
+      kf3 =    1.0
+      kf4 =   -0.001005
+
+      ! Set coeff. for 1. dissoc. of boric acid (Edmond and Gieskes, 1970)
+      ! ------------------------------------------------------------------
+      cb0  = -8966.90
+      cb1  = -2890.53
+      cb2  =   -77.942
+      cb3  =     1.728
+      cb4  =    -0.0996
+      cb5  =   148.0248
+      cb6  =   137.1942
+      cb7  =     1.62142
+      cb8  =   -24.4344
+      cb9  =   -25.085
+      cb10 =    -0.2474
+      cb11 =     0.053105
+
+      ! Set coeff. for dissoc. of water (Dickson and Riley, 1979, 
+      !   eq. 7, coefficient cw2 corrected from 0.9415 to 0.09415 
+      !   after pers. commun. to B. Bacastow, 1988)
+      ! ---------------------------------------------------------
+      cw0 = -13847.26
+      cw1 =    148.9652
+      cw2 =    -23.6521
+      cw3 =    118.67
+      cw4 =     -5.977
+      cw5 =      1.0495
+      cw6 =     -0.01615
+
+      ! Set coeff. for dissoc. of phosphate (Millero (1974)
+      ! ---------------------------------------------------
+      cp10 =   115.54
+      cp11 = -4576.752
+      cp12 =   -18.453
+      cp13 =  -106.736
+      cp14 =     0.69171
+      cp15 =    -0.65643
+      cp16 =    -0.01844
+      !
+      cp20 =   172.1033
+      cp21 = -8814.715
+      cp22 =   -27.927
+      cp23 =  -160.340
+      cp24 =     1.3566
+      cp25 =     0.37335
+      cp26 =    -0.05778
+      !
+      cp30 =   -18.126
+      cp31 = -3070.75
+      cp32 =    17.27039
+      cp33 =     2.81197
+      cp34 =   -44.99486
+      cp35 =    -0.09984
+
+      ! Set coeff. for dissoc. of phosphate (Millero (1974)
+      ! ---------------------------------------------------
+      cs10 =   117.385
+      cs11 = -8904.2
+      cs12 =   -19.334
+      cs13 =  -458.79
+      cs14 =     3.5913
+      cs15 =   188.74
+      cs16 =    -1.5998
+      cs17 =   -12.1652
+      cs18 =     0.07871
+      cs19 =     0.e0
+      cs20 =     1.e0
+      cs21 =    -0.001005
+
+
+      ! Set volumetric solubility constants for o2 in ml/l (Weiss, 1970)
+      ! ----------------------------------------------------------------
+      ox0 = -58.3877
+      ox1 =  85.8079
+      ox2 =  23.8439
+      ox3 =  -0.034892
+      ox4 =   0.015568
+      ox5 =  -0.0019387
+
+      !  FROM THE NEW BIOOPTIC MODEL PROPOSED JM ANDRE, WE READ HERE
+      !  A PRECOMPUTED ARRAY CORRESPONDING TO THE ATTENUATION COEFFICIENT
+
+      CALL ctlopn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
+         &           1, numout, .TRUE., 1 )
+      DO ichl = 1,61
+         READ(numlight,*) ztoto, ( xkrgb(iband,ichl), iband = 1,3 )
+      END DO
+      CLOSE(numlight)
+
+
+      CALL p4zche        ! initialize the chemical constants
+
+
+      ndayflxtr = 0      !  Initialize a counter for the computation of chemistry
+
+
+      IF(lwp) WRITE(numout,*) ' Initialisation of PISCES done'
+      !
+   END SUBROUTINE trc_ini_pisces
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                            No PISCES biochemical model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_ini_pisces             ! Empty routine
+   END SUBROUTINE trc_ini_pisces
+#endif
+
+   !!======================================================================
+END MODULE trcini_pisces
+
Index: branches/dev_001_GM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/trcini.F90	(revision 765)
+++ branches/dev_001_GM/NEMO/TOP_SRC/trcini.F90	(revision 766)
@@ -8,15 +8,11 @@
    !!              -   !  2005-10  (C. Ethe) Module
    !!----------------------------------------------------------------------
-#if defined key_passivetrc
-   !!----------------------------------------------------------------------
-   !!   'key_passivetrc'                                    Passive tracers
-   !!----------------------------------------------------------------------
    !!   trc_ini :   Initialization for passive tracer
    !!----------------------------------------------------------------------
-   USE oce_trc
-   USE trc
-   USE sms
-   USE lib_mpp
-   USE lbclnk
+   USE par_trc         ! TOP parameters
+   USE trcini_cfc      ! CFC     initialisation
+   USE trcini_lobster  ! LOBSTER initialisation
+   USE trcini_pisces   ! PISCES  initialisation
+   USE trcini_my_trc   ! MY_TRC  initialisation
 
    IMPLICIT NONE
@@ -25,26 +21,4 @@
    PUBLIC   trc_ini   ! called by ???
 
-# if defined key_trc_lobster1
-   !!----------------------------------------------------------------------
-   !!   'key_trc_lobster1'                        LOBSTER1 biological model  
-   !!----------------------------------------------------------------------
-#  include "trcini.lobster1.h90"
-
-# elif defined key_trc_pisces
-   !!----------------------------------------------------------------------
-   !!   'key_trc_pisces'                            PISCES biological model                  
-   !!----------------------------------------------------------------------
-#  include "trcini.pisces.h90"
-
-# elif defined key_cfc
-   !!----------------------------------------------------------------------
-   !!   'key_cfc  '                                          CFC model                  
-   !!----------------------------------------------------------------------
-#  include "trcini.cfc.h90"
-
-# else
-   !!----------------------------------------------------------------------
-   !!   Default option                               
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
@@ -59,60 +33,32 @@
       !!                    ***  ROUTINE trc_ini  ***
       !!              
-      !! ** Purpose :   Initialization of all passive tracer to zero
-      !!              (default case)
-      !!-------------------------------------------------------------------
-      INTEGER ::   ji ,jj ,jk ,jn, jl        ! dummy loop indices  
+      !! ** Purpose :   Initialization of passive tracer to zero
+      !!
+      !! ** Method  :   call the initialisation of each defined tracer
+      !!             model (LOBSTER, PISCES, CFC, MY_TRC)
       !!-------------------------------------------------------------------
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'trc_ini : initial set up of the passive tracers'
+      IF(lwp) WRITE(numout,*) 'trc_ini : initial set up of the tracer models'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
+      !
+      IF( lk_trc_lobster ) THEN   ;   CALL trc_ini_lobster      ! LOBSTER bio-model
+      ELSE                        ;   IF(lwp) WRITE(numout,*) '          LOBSTER not used'
+      ENDIF
 
-      ! initialization of passives tracers field
-      ! ----------------------------------------
-      DO jn = 1, jptra
-         trn(:,:,:,jn) = 0.e0
-         tra(:,:,:,jn) = 0.e0
-      END DO
+      IF( lk_trc_pisces  ) THEN   ;   CALL trc_init_pisces      ! PISCES  bio-model
+      ELSE                        ;   IF(lwp) WRITE(numout,*) '          PISCES not used'
+      ENDIF
 
-#  if defined key_trc_diaadd
-      ! initialization of output 2d and 3d arrays
-      DO jn = 1, jpdia2d
-         trc2d(:,:,jn) = 0.e0
-      END DO
-      DO jn = 1, jpdia3d
-         trc3d(:,:,:,jn) = 0.e0
-      END DO
-#  endif
+      IF( lk_trc_cfc     ) THEN   ;   CALL trc_init_cfc         ! CFC     tracers
+      ELSE                        ;   IF(lwp) WRITE(numout,*) '          CFC not used'
+      ENDIF
 
-#  if defined key_trc_diabio
-      ! initialization of biological trends
-      DO jn = 1, jpdiabio
-         trbio(:,:,:,jn) = 0.e0
-      END DO
-#  endif
-
-#  if defined key_trc_diatrd
-      ! initialization of tracer trends
-      DO jl = 1, jpdiatrc
-         DO jn = 1, jptra
-            IF( luttrd(jn) )   trtrd(:,:,:,ikeep(jn),jl) = 0.e0
-         END DO
-      END DO
-#  endif      
+      IF( lk_trc_my_trc  ) THEN   ;   CALL trc_init_my_trc      ! MY_TRC  tracers
+      ELSE                        ;   IF(lwp) WRITE(numout,*) '          MY_TRC not used'
+      ENDIF
       !
    END SUBROUTINE trc_ini
 
-# endif
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                      NO passive tracer
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_ini              ! Empty routine
-   END SUBROUTINE trc_ini
-#endif
-
    !!======================================================================
 END MODULE trcini
