Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 4990)
@@ -25,5 +25,5 @@
    USE trabbl              ! 
    USE prtctl_trc          ! Print control for debbuging
-   USE trdmod_oce
+   USE trd_oce
    USE trdtra
 
@@ -93,5 +93,5 @@
         DO jn = 1, jptra
            ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
-           CALL trd_tra( kt, 'TRC', jn, jptra_trd_bbl, ztrtrd(:,:,:,jn) )
+           CALL trd_tra( kt, 'TRC', jn, jptra_bbl, ztrtrd(:,:,:,jn) )
         END DO
         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd ) ! temporary save of trends
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 4990)
@@ -23,5 +23,5 @@
    USE prtctl_trc      ! Print control for debbuging
    USE trdtra
-   USE trdmod_oce
+   USE trd_oce
 
    IMPLICIT NONE
@@ -75,5 +75,5 @@
       !! ** Action  : - update the tracer trends tra with the newtonian 
       !!                damping trends.
-      !!              - save the trends ('key_trdmld_trc')
+      !!              - save the trends ('key_trdmxl_trc')
       !!----------------------------------------------------------------------
       !!
@@ -151,5 +151,5 @@
             IF( l_trdtrc ) THEN
                ztrtrd(:,:,:) = tra(:,:,:,jn) -  ztrtrd(:,:,:)
-               CALL trd_tra( kt, 'TRC', jn, jptra_trd_dmp, ztrtrd )
+               CALL trd_tra( kt, 'TRC', jn, jptra_dmp, ztrtrd )
             END IF
             !                                                       ! ===========
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 4990)
@@ -25,5 +25,5 @@
    USE traldf_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine)
    USE traldf_lap      ! lateral mixing            (tra_ldf_lap routine)
-   USE trdmod_oce
+   USE trd_oce
    USE trdtra
    USE prtctl_trc      ! Print control
@@ -73,16 +73,16 @@
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra        )  ! iso-level laplacian
       CASE ( 1 )                                                                                            ! rotated laplacian
                        IF( ln_traldf_grif ) THEN
                           CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
                        ELSE
-                          CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+                          CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra, rn_ahtb_0 )
                        ENDIF
-      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra            )  ! iso-level bilaplacian
       CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, nittrc000, 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
          !
       CASE ( -1 )                                     ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_lap   ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf0 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
@@ -90,9 +90,9 @@
             CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
          ELSE
-            CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+            CALL tra_ldf_iso     ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra, rn_ahtb_0 )
          ENDIF
          WRITE(charout, FMT="('ldf1 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, gtrui, gtrvi, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf2 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
@@ -105,5 +105,5 @@
         DO jn = 1, jptra
            ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
-           CALL trd_tra( kt, 'TRC', jn, jptra_trd_ldf, ztrtrd(:,:,:,jn) )
+           CALL trd_tra( kt, 'TRC', jn, jptra_ldf, ztrtrd(:,:,:,jn) )
         END DO
         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 4990)
@@ -30,5 +30,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl_trc      ! Print control for debbuging
-   USE trdmod_oce
+   USE trd_oce
    USE trdtra
    USE tranxt
@@ -148,5 +148,5 @@
                zfact = 1.e0 / r2dt(jk)  
                ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
-               CALL trd_tra( kt, 'TRC', jn, jptra_trd_atf, ztrdt )
+               CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt )
             END DO
          END DO
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 4990)
@@ -15,5 +15,5 @@
    USE oce_trc             ! ocean dynamics and tracers variables
    USE trc                 ! ocean passive tracers variables
-   USE trdmod_oce
+   USE trd_oce
    USE trdtra
    USE prtctl_trc          ! Print control for debbuging
@@ -156,6 +156,6 @@
                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
-               CALL trd_tra( kt, 'TRC', jn, jptra_trd_radb, ztrtrdb )       ! Asselin-like trend handling
-               CALL trd_tra( kt, 'TRC', jn, jptra_trd_radn, ztrtrdn )       ! standard     trend handling
+               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
+               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
               !
             ENDIF
@@ -187,6 +187,6 @@
                ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
                ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
-               CALL trd_tra( kt, 'TRC', jn, jptra_trd_radb, ztrtrdb )       ! Asselin-like trend handling
-               CALL trd_tra( kt, 'TRC', jn, jptra_trd_radn, ztrtrdn )       ! standard     trend handling
+               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
+               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
               !
             ENDIF
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 4990)
@@ -19,5 +19,5 @@
    USE trc             ! ocean  passive tracers variables
    USE prtctl_trc      ! Print control for debbuging
-   USE trdmod_oce
+   USE trd_oce
    USE trdtra
 
@@ -104,5 +104,5 @@
          IF( l_trdtrc ) THEN
             ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
-            CALL trd_tra( kt, 'TRC', jn, jptra_trd_nsr, ztrtrd )
+            CALL trd_tra( kt, 'TRC', jn, jptra_nsr, ztrtrd )
          END IF
          !                                                       ! ===========
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 4990)
@@ -82,5 +82,5 @@
       IF( .NOT. Agrif_Root())   CALL Agrif_Update_Trc( kstp )   ! Update tracer at AGRIF zoom boundaries : children only
 #endif
-         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
+         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi )  ! Partial steps: now horizontal gradient of passive
                                                                 ! tracers at the bottom ocean level
          !
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 4990)
@@ -19,5 +19,5 @@
    USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine)
    USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine)
-   USE trdmod_oce
+   USE trd_oce
    USE trdtra
    USE prtctl_trc      ! Print control
@@ -106,5 +106,5 @@
                ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dt(jk) ) - ztrtrd(:,:,jk,jn)
             END DO
-            CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:,jn) )
+            CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) )
          END DO
          CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90	(revision 4624)
+++ 	(revision )
@@ -1,1437 +1,0 @@
-MODULE trdmld_trc
-   !!======================================================================
-   !!                       ***  MODULE  trdmld_trc  ***
-   !! Ocean diagnostics:  mixed layer passive tracer trends 
-   !!======================================================================
-   !! History :  9.0  !  06-08  (C. Deltel)  Original code (from trdmld.F90)
-   !!                 !  07-04  (C. Deltel)  Bug fix : add trcrad trends
-   !!                 !  07-06  (C. Deltel)  key_gyre : do not call lbc_lnk
-   !!----------------------------------------------------------------------
-#if   defined key_top && ( defined key_trdmld_trc   ||   defined key_esopa )
-   !!----------------------------------------------------------------------
-   !!   'key_trdmld_trc'                      mixed layer trend diagnostics
-   !!----------------------------------------------------------------------
-   !!   trd_mld_trc      : passive tracer cumulated trends averaged over ML
-   !!   trd_mld_trc_zint : passive tracer trends vertical integration
-   !!   trd_mld_trc_init : initialization step
-   !!----------------------------------------------------------------------
-   USE trc               ! tracer definitions (trn, trb, tra, etc.)
-   USE trc_oce, ONLY :   nn_dttrc  ! frequency of step on passive tracers
-   USE dom_oce           ! domain definition
-   USE zdfmxl  , ONLY : nmln ! number of level in the mixed layer
-   USE zdf_oce , ONLY : avt  ! vert. diffusivity coef. at w-point for temp  
-# if defined key_zdfddm   
-   USE zdfddm  , ONLY : avs  ! salinity vertical diffusivity coeff. at w-point
-# endif
-   USE trcnam_trp        ! passive tracers transport namelist variables
-   USE trdmod_trc_oce    ! definition of main arrays used for trends computations
-   USE in_out_manager    ! I/O manager
-   USE dianam            ! build the name of file (routine)
-   USE ldfslp            ! iso-neutral slopes 
-   USE ioipsl            ! NetCDF library
-   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp           ! MPP library
-   USE trdmld_trc_rst    ! restart for diagnosing the ML trends
-   USE prtctl            ! print control
-   USE sms_pisces        ! PISCES bio-model
-   USE wrk_nemo          ! Memory allocation
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC trd_mld_trc
-   PUBLIC trd_mld_trc_alloc
-   PUBLIC trd_mld_bio
-   PUBLIC trd_mld_trc_init
-   PUBLIC trd_mld_trc_zint
-   PUBLIC trd_mld_bio_zint
-
-   CHARACTER (LEN=40) ::  clhstnam                                ! name of the trends NetCDF file
-   INTEGER ::   nmoymltrd
-   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::   ndextrd1
-   INTEGER, DIMENSION(jptra) ::   nidtrd, nh_t
-   INTEGER ::   ndimtrd1                        
-   INTEGER, SAVE ::  ionce, icount
-#if defined key_pisces_reduced
-   INTEGER ::   nidtrdbio, nh_tb
-   INTEGER, SAVE ::  ioncebio, icountbio
-   INTEGER, SAVE ::   nmoymltrdbio
-#endif
-   LOGICAL :: llwarn  = .TRUE.                                    ! this should always be .TRUE.
-   LOGICAL :: lldebug = .TRUE.
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztmltrd2   !
-#if defined key_pisces_reduced
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ztmltrdbio2  ! only needed for mean diagnostics in trd_mld_bio()
-#endif
-
-   !! * Substitutions
-#  include "top_substitute.h90"
-#  include "zdfddm_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header:  $ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION trd_mld_trc_alloc()
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mld_trc_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) ,      &
-#if defined key_pisces_reduced
-         &      ztmltrdbio2(jpi,jpj,jpdiabio)      ,      &
-#endif
-         &      ndextrd1(jpi*jpj)                  ,  STAT=trd_mld_trc_alloc)
-         !
-      IF( lk_mpp                )   CALL mpp_sum ( trd_mld_trc_alloc )
-      IF( trd_mld_trc_alloc /=0 )   CALL ctl_warn('trd_mld_trc_alloc: failed to allocate arrays')
-      !
-   END FUNCTION trd_mld_trc_alloc
-
-
-   SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mld_trc_zint  ***
-      !! 
-      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments 
-      !!                to the subroutine. This vertical average is performed from ocean
-      !!                surface down to a chosen control surface.
-      !!
-      !! ** Method/usage :
-      !!      The control surface can be either a mixed layer depth (time varying)
-      !!      or a fixed surface (jk level or bowl). 
-      !!      Choose control surface with nctls_trc in namelist NAMTRD :
-      !!        nctls_trc = -2  : use isopycnal surface
-      !!        nctls_trc = -1  : use euphotic layer with light criterion 
-      !!        nctls_trc =  0  : use mixed layer with density criterion 
-      !!        nctls_trc =  1  : read index from file 'ctlsurf_idx'
-      !!        nctls_trc >  1  : use fixed level surface jk = nctls_trc
-      !!      Note: in the remainder of the routine, the volume between the 
-      !!            surface and the control surface is called "mixed-layer"
-      !!----------------------------------------------------------------------
-      !!
-      INTEGER, INTENT( in ) ::   ktrd, kjn                        ! ocean trend index and passive tracer rank
-      CHARACTER(len=2), INTENT( in ) ::  ctype                    ! surface/bottom (2D) or interior (3D) physics
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  ptrc_trdmld ! passive tracer trend
-      !
-      INTEGER ::   ji, jj, jk, isum
-      REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, zvlmsk )
-
-      ! I. Definition of control surface and integration weights
-      ! --------------------------------------------------------
-
-      ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
-         !
-         tmltrd_trc(:,:,:,:) = 0.e0                               ! <<< reset trend arrays to zero
-          
-         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
-         SELECT CASE ( nn_ctls_trc )                                ! choice of the control surface
-            CASE ( -2  )   ;   STOP 'trdmld_trc : not ready '     !     -> isopycnal surface (see ???)
-#if defined key_pisces || defined key_pisces_reduced
-            CASE ( -1  )   ;   nmld_trc(:,:) = neln(:,:)          !     -> euphotic layer with light criterion
-#endif
-            CASE (  0  )   ;   nmld_trc(:,:) = nmln(:,:)          !     -> ML with density criterion (see zdfmxl)
-            CASE (  1  )   ;   nmld_trc(:,:) = nbol_trc(:,:)          !     -> read index from file
-            CASE (  2: )   ;   nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
-                               nmld_trc(:,:) = nn_ctls_trc + 1      !     -> model level
-         END SELECT
-
-         ! ... Compute ndextrd1 and ndimtrd1  ??? role de jpktrd_trc
-         IF( ionce == 1 ) THEN
-            !
-            isum  = 0   ;   zvlmsk(:,:) = 0.e0
-
-            IF( jpktrd_trc < jpk ) THEN                           ! description ???
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
-                        zvlmsk(ji,jj) = tmask(ji,jj,1)
-                     ELSE
-                        isum = isum + 1
-                        zvlmsk(ji,jj) = 0.e0
-                     ENDIF
-                  END DO
-               END DO
-            ENDIF
-
-            IF( isum > 0 ) THEN                                   ! index of ocean points (2D only)
-               WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum 
-               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
-            ELSE 
-               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
-            ENDIF                                
-
-            ionce = 0                                             ! no more pass here
-            !
-         ENDIF   ! ionce == 1
-         
-         ! ... Weights for vertical averaging
-         wkx_trc(:,:,:) = 0.e0
-         DO jk = 1, jpktrd_trc                                    ! initialize wkx_trc with vertical scale factor in mixed-layer
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         
-         rmld_trc(:,:) = 0.e0
-         DO jk = 1, jpktrd_trc                                    ! compute mixed-layer depth : rmld_trc
-            rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
-         END DO
-         
-         DO jk = 1, jpktrd_trc                                    ! compute integration weights
-            wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
-         END DO
-
-         icount = 0                                               ! <<< flag = off : control surface & integr. weights
-         !                                                        !     computed only once per time step
-      ENDIF ONCE_PER_TIME_STEP
-
-      ! II. Vertical integration of trends in the mixed-layer
-      ! -----------------------------------------------------
-
-      SELECT CASE ( ctype )
-         CASE ( '3D' )                                            ! mean passive tracer trends in the mixed-layer
-            DO jk = 1, jpktrd_trc
-               tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,jk) * wkx_trc(:,:,jk)   
-            END DO
-         CASE ( '2D' )                                            ! forcing at upper boundary of the mixed-layer 
-            tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,1) * wkx_trc(:,:,1)  ! non penetrative
-      END SELECT
-      !
-      CALL wrk_dealloc( jpi, jpj, zvlmsk )
-      !
-   END SUBROUTINE trd_mld_trc_zint
-
-
-   SUBROUTINE trd_mld_bio_zint( ptrc_trdmld, ktrd )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mld_bio_zint  ***
-      !!
-      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments
-      !!                to the subroutine. This vertical average is performed from ocean
-      !!                surface down to a chosen control surface.
-      !!
-      !! ** Method/usage :
-      !!      The control surface can be either a mixed layer depth (time varying)
-      !!      or a fixed surface (jk level or bowl).
-      !!      Choose control surface with nctls in namelist NAMTRD :
-      !!        nctls_trc = 0  : use mixed layer with density criterion
-      !!        nctls_trc = 1  : read index from file 'ctlsurf_idx'
-      !!        nctls_trc > 1  : use fixed level surface jk = nctls_trc
-      !!      Note: in the remainder of the routine, the volume between the
-      !!            surface and the control surface is called "mixed-layer"
-      !!----------------------------------------------------------------------
-      !!
-      INTEGER                         , INTENT(in) ::   ktrd          ! bio trend index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   ptrc_trdmld   ! passive trc trend
-#if defined key_pisces_reduced
-      !
-      INTEGER ::   ji, jj, jk, isum
-      REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, zvlmsk )
-
-      ! I. Definition of control surface and integration weights
-      ! --------------------------------------------------------
-      !            ==> only once per time step <==
-
-      IF( icountbio == 1 ) THEN
-         !
-         tmltrd_bio(:,:,:) = 0.e0    ! <<< reset trend arrays to zero
-         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
-         SELECT CASE ( nn_ctls_trc )                                    ! choice of the control surface
-            CASE ( -2  )   ;   STOP 'trdmld_trc : not ready '     !     -> isopycnal surface (see ???)
-            CASE ( -1  )   ;   nmld_trc(:,:) = neln(:,:)          !     -> euphotic layer with light criterion
-            CASE (  0  )   ;   nmld_trc(:,:) = nmln(:,:)          !     -> ML with density criterion (see zdfmxl)
-            CASE (  1  )   ;   nmld_trc(:,:) = nbol_trc(:,:)          !     -> read index from file
-            CASE (  2: )   ;   nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
-                               nmld_trc(:,:) = nn_ctls_trc + 1          !     -> model level
-         END SELECT
-
-         ! ... Compute ndextrd1 and ndimtrd1 only once
-         IF( ioncebio == 1 ) THEN
-            !
-            ! Check of validity : nmld_trc(ji,jj) <= jpktrd_trc
-            isum        = 0
-            zvlmsk(:,:) = 0.e0
-
-            IF( jpktrd_trc < jpk ) THEN
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
-                        zvlmsk(ji,jj) = tmask(ji,jj,1)
-                     ELSE
-                        isum = isum + 1
-                        zvlmsk(ji,jj) = 0.
-                     END IF
-                  END DO
-               END DO
-            END IF
-
-            ! Index of ocean points (2D only)
-            IF( isum > 0 ) THEN
-               WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
-               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
-            ELSE
-               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
-            END IF
-
-            ioncebio = 0                  ! no more pass here
-            !
-         END IF !  ( ioncebio == 1 )
-
-         ! ... Weights for vertical averaging
-         wkx_trc(:,:,:) = 0.e0
-         DO jk = 1, jpktrd_trc         ! initialize wkx_trc with vertical scale factor in mixed-layer
-            DO jj = 1,jpj
-              DO ji = 1,jpi
-                  IF( jk - nmld_trc(ji,jj) < 0. )   wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-
-         rmld_trc(:,:) = 0.
-         DO jk = 1, jpktrd_trc         ! compute mixed-layer depth : rmld_trc
-            rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
-         END DO
-
-         DO jk = 1, jpktrd_trc         ! compute integration weights
-            wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
-         END DO
-
-         icountbio = 0                    ! <<< flag = off : control surface & integr. weights
-         !                             !     computed only once per time step
-      END IF ! ( icountbio == 1 )
-
-      ! II. Vertical integration of trends in the mixed-layer
-      ! -----------------------------------------------------
-
-
-      DO jk = 1, jpktrd_trc
-         tmltrd_bio(:,:,ktrd) = tmltrd_bio(:,:,ktrd) + ptrc_trdmld(:,:,jk) * wkx_trc(:,:,jk)
-      END DO
-
-      CALL wrk_dealloc( jpi, jpj, zvlmsk )
-#endif
-      !
-   END SUBROUTINE trd_mld_bio_zint
-
-
-   SUBROUTINE trd_mld_trc( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mld_trc  ***
-      !! 
-      !! ** Purpose :  Compute and cumulate the mixed layer trends over an analysis
-      !!               period, and write NetCDF (or dimg) outputs.
-      !!
-      !! ** Method/usage :
-      !!          The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant 
-      !!          logical namelist variable) :
-      !!          1) to explain the difference between initial and final 
-      !!             mixed-layer T & S (where initial and final relate to the
-      !!             current analysis window, defined by ntrc_trc in the namelist)
-      !!          2) to explain the difference between the current and previous 
-      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
-      !!             performed over each analysis window).
-      !!
-      !! ** Consistency check : 
-      !!        If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
-      !!        entrainment) should be zero, at machine accuracy. Note that in the case
-      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
-      !!        over the first two analysis windows (except if restart).
-      !!        N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
-      !!             for checking residuals.
-      !!             On a NEC-SX5 computer, this typically leads to:
-      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false.
-      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true.
-      !!
-      !! ** Action :
-      !!       At each time step, mixed-layer averaged trends are stored in the 
-      !!       tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
-      !!       This array is known when trd_mld is called, at the end of the stp subroutine, 
-      !!       except for the purely vertical K_z diffusion term, which is embedded in the
-      !!       lateral diffusion trend.
-      !!
-      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
-      !!       from the lateral diffusion trend.
-      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
-      !!       arrays are updated.
-      !!       In III), called only once per analysis window, we compute the total trends,
-      !!       along with the residuals and the Asselin correction terms.
-      !!       In IV), the appropriate trends are written in the trends NetCDF file.
-      !!
-      !! References :
-      !!       - Vialard & al.
-      !!       - See NEMO documentation (in preparation)
-      !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !
-      INTEGER ::   ji, jj, jk, jl, ik, it, itmod, jn
-      REAL(wp) ::   zavt, zfn, zfn2
-      !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmltot             ! d(trc)/dt over the anlysis window (incl. Asselin)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlres             ! residual = dh/dt entrainment term
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlatf             ! for storage only
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlrad             ! for storage only (for trb<0 corr in trcrad)
-      !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmltot2            ! -+
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlres2            !  | working arrays to diagnose the trends
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmltrdm2           !  | associated with the time meaned ML
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlatf2            !  | passive tracers
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlrad2            !  | (-> for trb<0 corr in trcrad)
-      !
-      CHARACTER (LEN=10) ::   clvar
-#if defined key_dimgout
-      INTEGER ::   iyear,imon,iday
-      CHARACTER(LEN=80) ::   cltext, clmode
-#endif
-      !!----------------------------------------------------------------------
-
-      ! Set-up pointers into sub-arrays of workspaces
-      CALL wrk_alloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad             )
-      CALL wrk_alloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
-
-      IF( nn_dttrc  /= 1  )   CALL ctl_stop( " Be careful, trends diags never validated " )
-
-      ! ======================================================================
-      ! I. Diagnose the purely vertical (K_z) diffusion trend
-      ! ======================================================================
-
-      ! ... These terms can be estimated by flux computation at the lower boundary of the ML 
-      !     (we compute (-1/h) * K_z * d_z( tracer ))
-
-      IF( ln_trcldf_iso ) THEN
-         !
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               ik = nmld_trc(ji,jj)
-               zavt = fsavs(ji,jj,ik)
-               DO jn = 1, jptra
-                  IF( ln_trdtrc(jn) )    &
-                  tmltrd_trc(ji,jj,jpmld_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
-                       &                    * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) )            &
-                       &                    / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
-               END DO
-            END DO
-         END DO
-
-         DO jn = 1, jptra
-            ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
-            IF( ln_trdtrc(jn) ) &
-                 tmltrd_trc(:,:,jpmld_trc_ldf,jn) = tmltrd_trc(:,:,jpmld_trc_ldf,jn) - tmltrd_trc(:,:,jpmld_trc_zdf,jn)
-   
-         END DO
-         !      
-      ENDIF
-
-      IF ( cp_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration
-      ! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm. 
-      ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               DO jl = 1, jpltrd_trc
-                  CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
-               END DO
-            ENDIF
-         END DO
-      ENDIF
-      ! ======================================================================
-      ! II. Cumulate the trends over the analysis window
-      ! ======================================================================
-
-      ztmltrd2(:,:,:,:) = 0.e0   ;   ztmltot2(:,:,:)   = 0.e0     ! <<< reset arrays to zero
-      ztmlres2(:,:,:)   = 0.e0   ;   ztmlatf2(:,:,:)   = 0.e0
-      ztmlrad2(:,:,:)   = 0.e0
-
-      ! II.1 Set before values of vertically averages passive tracers
-      ! -------------------------------------------------------------
-      IF( kt > nittrc000 ) THEN
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               tmlb_trc   (:,:,jn) = tml_trc   (:,:,jn)
-               tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_atf,jn)
-               tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_radb,jn)
-            ENDIF
-         END DO
-      ENDIF
-
-      ! II.2 Vertically averaged passive tracers
-      ! ----------------------------------------
-      tml_trc(:,:,:) = 0.e0
-      DO jk = 1, jpktrd_trc ! - 1 ???
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) &
-               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
-         END DO
-      END DO
-
-      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window    
-      ! ------------------------------------------------------------------------
-      IF( kt == nittrc000 + nn_dttrc ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)    ???
-         !
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
-               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
-               
-               tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0   ;   tmltrd_atf_sumb_trc  (:,:,jn) = 0.e0
-               tmltrd_rad_sumb_trc  (:,:,jn) = 0.e0
-            ENDIF
-         END DO
-         
-         rmldbn_trc(:,:) = rmld_trc(:,:)
-         !
-      ENDIF
-
-      ! II.4 Cumulated trends over the analysis period
-      ! ----------------------------------------------
-      !
-      !         [  1rst analysis window ] [     2nd analysis window     ]                       
-      !
-      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
-      !                            ntrd                             2*ntrd       etc.
-      !     1      2     3     4    =5 e.g.                          =10
-      !
-      IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN                        ! ???
-         !
-         nmoymltrd = nmoymltrd + 1
-
-
-         ! ... Cumulate over BOTH physical contributions AND over time steps
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               DO jl = 1, jpltrd_trc
-                  tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
-               END DO
-            ENDIF
-         END DO
-
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               ! ... Special handling of the Asselin trend 
-               tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
-               tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
-
-               ! ... Trends associated with the time mean of the ML passive tracers
-               tmltrd_sum_trc    (:,:,:,jn) = tmltrd_sum_trc    (:,:,:,jn) + tmltrd_trc    (:,:,:,jn)
-               tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
-               tml_sum_trc       (:,:,jn)   = tml_sum_trc       (:,:,jn)   + tml_trc       (:,:,jn)
-            ENDIF
-         ENDDO
-
-         rmld_sum_trc      (:,:)     = rmld_sum_trc      (:,:)     + rmld_trc      (:,:)
-         !
-      ENDIF
-
-      ! ======================================================================
-      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
-      ! ======================================================================
-
-      ! Convert to appropriate physical units
-      tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
-
-      itmod = kt - nittrc000 + 1
-      it    = kt
-
-      MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN           ! nitend MUST be multiple of nn_trd_trc
-         !
-         ztmltot (:,:,:) = 0.e0                                   ! reset arrays to zero
-         ztmlres (:,:,:) = 0.e0
-         ztmltot2(:,:,:) = 0.e0
-         ztmlres2(:,:,:) = 0.e0
-      
-         zfn  = FLOAT( nmoymltrd )    ;    zfn2 = zfn * zfn
-         
-         ! III.1 Prepare fields for output ("instantaneous" diagnostics) 
-         ! -------------------------------------------------------------
-
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               !-- Compute total trends    (use rdttrc instead of rdt ???)
-               IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN  ! EULER-FORWARD schemes
-                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
-               ELSE                                                                     ! LEAP-FROG schemes
-                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
-               ENDIF
-               
-               !-- Compute residuals
-               ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
-                  &                                                 - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)  )
-               
-               !-- Diagnose Asselin trend over the analysis window 
-               ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
-               ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
-               
-         !-- Lateral boundary conditions
-               IF ( cp_cfg .NE. 'gyre' ) THEN
-                  CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
-                  CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
-               ENDIF
-
-
-#if defined key_diainstant
-               STOP 'tmltrd_trc : key_diainstant was never checked within trdmld. Comment this to proceed.'
-#endif
-            ENDIF
-         END DO
-
-         ! III.2 Prepare fields for output ("mean" diagnostics) 
-         ! ----------------------------------------------------
-         
-         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
-         rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
-
-               !-- Compute passive tracer total trends
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
-               ztmltot2   (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) /  ( 2.*rdt )    ! now tracer unit is /sec
-            ENDIF
-         END DO
-
-         !-- Compute passive tracer residuals
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN
-               !
-               DO jl = 1, jpltrd_trc
-                  ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
-               END DO
-               
-               ztmltrdm2(:,:,jn) = 0.e0
-               DO jl = 1, jpltrd_trc
-                  ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
-               END DO
-               
-               ztmlres2(:,:,jn) =  ztmltot2(:,:,jn)  -       &
-                  & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
-                  &                     - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
-               !
-
-               !-- Diagnose Asselin trend over the analysis window
-               ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) &
-                  &                                               + tmltrd_atf_sumb_trc(:,:,jn)
-               ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) &
-                  &                                               + tmltrd_rad_sumb_trc(:,:,jn)
-
-         !-- Lateral boundary conditions 
-               IF ( cp_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration    
-                  CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
-                  CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
-                  DO jl = 1, jpltrd_trc
-                     CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
-                  END DO
-               ENDIF
-
-            ENDIF
-         END DO
-
-         ! * Debugging information *
-         IF( lldebug ) THEN
-            !
-            WRITE(numout,*) 'trd_mld_trc : write trends in the Mixed Layer for debugging process:'
-            WRITE(numout,*) '~~~~~~~~~~~  '
-            WRITE(numout,*)
-            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
-
-            DO jn = 1, jptra
-
-               IF( ln_trdtrc(jn) ) THEN
-                  WRITE(numout, *)
-                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
-                  
-                  WRITE(numout, *)
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
-                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
-                  
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
-                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
-                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
-                  
-                  WRITE(numout, *)
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
-                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_atf,jn)) 
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_radb,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
-                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
-                  
-                  WRITE(numout, *)
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
-                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
-                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
-                  
-                  WRITE(numout, *)
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
-                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
-                  
-                  WRITE(numout, *)
-                  DO jl = 1, jpltrd_trc
-                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmld_trc_xxx = ', jl, &
-                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
-                  END DO
-               
-                  WRITE(numout,*) 
-                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
-                  WRITE(numout,*)
-                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
-                  WRITE(numout,*)
-                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
-                  
-                  WRITE(numout,*) 
-                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
-                  WRITE(numout,*)
-                  
-                  WRITE(numout,*) '...................................................'
-                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
-                  DO jj = 5, 1, -1
-                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
-                  END DO
-                  
-                  WRITE(numout,*) 
-                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
-                  WRITE(numout,*)
-                  
-                  WRITE(numout,*) '...................................................'
-                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
-                  DO jj = 5, 1, -1
-                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
-                  END DO
-                  !
-               ENDIF
-               !
-            END DO
-
-
-97            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
-98            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
-99            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
-              WRITE(numout,*)
-            !
-         ENDIF
-
-         ! III.3 Time evolution array swap
-         ! -------------------------------
-         ! ML depth
-         rmldbn_trc(:,:)   = rmld_trc(:,:)
-         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
-         DO jn = 1, jptra
-            IF( ln_trdtrc(jn) ) THEN        
-               ! For passive tracer instantaneous diagnostics
-               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
-               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
-               
-               ! For passive tracer mean diagnostics
-               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
-               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
-               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn)
-               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmld_trc_radb,jn)
-               
-               
-               ! III.4 Convert to appropriate physical units
-               ! -------------------------------------------
-               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * rn_ucf_trc/zfn   ! instant diags
-               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * rn_ucf_trc/zfn
-               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * rn_ucf_trc/zfn
-               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * rn_ucf_trc/zfn
-               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
-               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * rn_ucf_trc/zfn2
-               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * rn_ucf_trc/zfn2
-               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * rn_ucf_trc/zfn2
-               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * rn_ucf_trc/zfn2
-               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * rn_ucf_trc/zfn2
-            ENDIF
-         END DO
-         !
-      ENDIF MODULO_NTRD
-
-      ! ======================================================================
-      ! IV. Write trends in the NetCDF file
-      ! ======================================================================
-
-      ! IV.1 Code for dimg mpp output
-      ! -----------------------------
-
-# if defined key_dimgout
-      STOP 'Not implemented'
-# else
-      
-      ! IV.2 Code for IOIPSL/NetCDF output
-      ! ----------------------------------
-
-      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
-         WRITE(numout,*) ' '
-         WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :'
-         WRITE(numout,*) '~~~~~~~~~~~ '
-         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
-         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
-         WRITE(numout,*) ' '
-      ENDIF
-         
-      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
-         !
-
-         DO jn = 1, jptra
-            !
-            IF( ln_trdtrc(jn) ) THEN
-               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
-               !-- Output the fields
-               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
-               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
-               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
-               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
-           
-               DO jl = 1, jpltrd_trc - 2
-                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),             &
-                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
-               END DO
-
-               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
-                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
-
-               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
-                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
-                     
-            ENDIF
-         END DO
-
-         IF( kt == nitend ) THEN 
-            DO jn = 1, jptra
-               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
-            END DO
-         ENDIF
-
-      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
-         
-         DO jn = 1, jptra
-            !
-            IF( ln_trdtrc(jn) ) THEN
-               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
-               !-- Output the fields
-               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
-
-               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
-               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
-               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
-
-               DO jl = 1, jpltrd_trc - 2
-                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),           &
-                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
-               END DO
-            
-               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
-                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
-
-               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
-                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
-
-            ENDIF 
-            !
-         END DO
-         IF( kt == nitend ) THEN 
-            DO jn = 1, jptra
-               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
-            END DO
-         ENDIF
-
-         !
-      ENDIF NETCDF_OUTPUT
-         
-      ! Compute the control surface (for next time step) : flag = on
-      icount = 1
-
-# endif /* key_dimgout */
-
-      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
-         !
-         ! Reset cumulative arrays to zero
-         ! -------------------------------         
-         nmoymltrd = 0
-         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
-         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
-         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
-         rmld_sum_trc       (:,:)     = 0.e0
-         !
-      ENDIF
-
-      ! ======================================================================
-      ! V. Write restart file
-      ! ======================================================================
-
-      IF( lrst_trc )   CALL trd_mld_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
-
-      CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad             )
-      CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
-      !
-   END SUBROUTINE trd_mld_trc
-
-
-   SUBROUTINE trd_mld_bio( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mld  ***
-      !!
-      !! ** Purpose :  Compute and cumulate the mixed layer biological trends over an analysis
-      !!               period, and write NetCDF (or dimg) outputs.
-      !!
-      !! ** Method/usage :
-      !!          The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant
-      !!          logical namelist variable) :
-      !!          1) to explain the difference between initial and final
-      !!             mixed-layer T & S (where initial and final relate to the
-      !!             current analysis window, defined by ntrd in the namelist)
-      !!          2) to explain the difference between the current and previous
-      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
-      !!             performed over each analysis window).
-      !!
-      !! ** Consistency check :
-      !!        If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
-      !!        entrainment) should be zero, at machine accuracy. Note that in the case
-      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
-      !!        over the first two analysis windows (except if restart).
-      !!        N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
-      !!             for checking residuals.
-      !!             On a NEC-SX5 computer, this typically leads to:
-      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false.
-      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true.
-      !!
-      !! ** Action :
-      !!       At each time step, mixed-layer averaged trends are stored in the
-      !!       tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
-      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
-      !!       except for the purely vertical K_z diffusion term, which is embedded in the
-      !!       lateral diffusion trend.
-      !!
-      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
-      !!       from the lateral diffusion trend.
-      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
-      !!       arrays are updated.
-      !!       In III), called only once per analysis window, we compute the total trends,
-      !!       along with the residuals and the Asselin correction terms.
-      !!       In IV), the appropriate trends are written in the trends NetCDF file.
-      !!
-      !! References :
-      !!       - Vialard & al.
-      !!       - See NEMO documentation (in preparation)
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt                       ! ocean time-step index
-#if defined key_pisces_reduced
-      INTEGER  ::  jl, it, itmod
-      LOGICAL  :: llwarn  = .TRUE., lldebug = .TRUE.
-      REAL(wp) :: zfn, zfn2
-#if defined key_dimgout
-      INTEGER ::  iyear,imon,iday
-      CHARACTER(LEN=80) :: cltext, clmode
-#endif
-      !!----------------------------------------------------------------------
-      ! ... Warnings
-      IF( nn_dttrc  /= 1  ) CALL ctl_stop( " Be careful, trends diags never validated " )
-
-      ! ======================================================================
-      ! II. Cumulate the trends over the analysis window
-      ! ======================================================================
-
-      ztmltrdbio2(:,:,:) = 0.e0  ! <<< reset arrays to zero
-
-      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
-      ! ------------------------------------------------------------------------
-      IF( kt == nittrc000 + nn_dttrc ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
-         !
-         tmltrd_csum_ub_bio (:,:,:) = 0.e0
-         !
-      END IF
-
-      ! II.4 Cumulated trends over the analysis period
-      ! ----------------------------------------------
-      !
-      !         [  1rst analysis window ] [     2nd analysis window     ]
-      !
-      !
-      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
-      !                            ntrd                             2*ntrd       etc.
-      !     1      2     3     4    =5 e.g.                          =10
-      !
-      IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
-         !
-         nmoymltrdbio = nmoymltrdbio + 1
-
-         ! ... Trends associated with the time mean of the ML passive tracers
-         tmltrd_sum_bio    (:,:,:) = tmltrd_sum_bio    (:,:,:) + tmltrd_bio    (:,:,:)
-         tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
-         !
-      END IF
-
-      ! ======================================================================
-      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
-      ! ======================================================================
-
-      ! Convert to appropriate physical units
-      tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc
-
-      MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN      ! nitend MUST be multiple of ntrd
-         !
-         zfn  = float(nmoymltrdbio)    ;    zfn2 = zfn * zfn
-
-         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
-         ! -------------------------------------------------------------
-
-#if defined key_diainstant
-         STOP 'tmltrd_bio : key_diainstant was never checked within trdmld. Comment this to proceed.'
-#endif
-         ! III.2 Prepare fields for output ("mean" diagnostics)
-         ! ----------------------------------------------------
-
-         ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
-
-         !-- Lateral boundary conditions
-         IF ( cp_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration 
-            ! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
-            DO jn = 1, jpdiabio
-              CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
-            ENDDO
-         ENDIF
-
-         IF( lldebug ) THEN
-            !
-            WRITE(numout,*) 'trd_mld_bio : write trends in the Mixed Layer for debugging process:'
-            WRITE(numout,*) '~~~~~~~~~~~  '
-            WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
-            WRITE(numout,*)
-
-            DO jl = 1, jpdiabio
-              IF( ln_trdmld_trc_instant ) THEN
-                  WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX  = ', jl, &
-                     & ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
-              ELSE
-                  WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX  = ', jl, &
-                     & ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
-              endif
-            END DO
-
-97          FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
-98          FORMAT(a10, i3, 2x, a30, 2x, g20.10)
-99          FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
-            WRITE(numout,*)
-            !
-         ENDIF
-
-         ! III.3 Time evolution array swap
-         ! -------------------------------
-
-         ! For passive tracer mean diagnostics
-         tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
-
-         ! III.4 Convert to appropriate physical units
-         ! -------------------------------------------
-         ztmltrdbio2    (:,:,:) = ztmltrdbio2    (:,:,:) * rn_ucf_trc/zfn2
-
-      END IF MODULO_NTRD
-
-      ! ======================================================================
-      ! IV. Write trends in the NetCDF file
-      ! ======================================================================
-
-      ! IV.1 Code for dimg mpp output
-      ! -----------------------------
-
-# if defined key_dimgout
-      STOP 'Not implemented'
-# else
-
-      ! IV.2 Code for IOIPSL/NetCDF output
-      ! ----------------------------------
-
-      ! define time axis
-      itmod = kt - nittrc000 + 1
-      it    = kt
-
-      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
-         WRITE(numout,*) ' '
-         WRITE(numout,*) 'trd_mld_bio : write ML bio trends in the NetCDF file :'
-         WRITE(numout,*) '~~~~~~~~~~~ '
-         WRITE(numout,*) '          ', TRIM(clhstnam), ' at kt = ', kt
-         WRITE(numout,*) '          N.B. nmoymltrdbio = ', nmoymltrdbio
-         WRITE(numout,*) ' '
-      END IF
-
-
-      ! 2. Start writing data
-      ! ---------------------
-
-      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN    ! <<< write the trends for passive tracer instant. diags
-         !
-            DO jl = 1, jpdiabio
-               CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) ,            &
-                    &          it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
-            END DO
-
-
-         IF( kt == nitend )   CALL histclo( nidtrdbio )
-
-      ELSE    ! <<< write the trends for passive tracer mean diagnostics
-
-            DO jl = 1, jpdiabio
-               CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) ,            &
-                    &          it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
-            END DO
-
-            IF( kt == nitend )   CALL histclo( nidtrdbio )
-            !
-      END IF NETCDF_OUTPUT
-
-      ! Compute the control surface (for next time step) : flag = on
-      icountbio = 1
-
-
-# endif /* key_dimgout */
-
-      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
-         !
-         ! III.5 Reset cumulative arrays to zero
-         ! -------------------------------------
-         nmoymltrdbio = 0
-         tmltrd_csum_ln_bio (:,:,:) = 0.e0
-         tmltrd_sum_bio     (:,:,:) = 0.e0
-      END IF
-
-      ! ======================================================================
-      ! Write restart file
-      ! ======================================================================
-
-! restart write is done in trd_mld_trc_write which is called by trd_mld_bio (Marina)
-!
-#endif
-   END SUBROUTINE trd_mld_bio
-
-
-   REAL FUNCTION sum2d( ztab )
-      !!----------------------------------------------------------------------
-      !! CD ??? prevoir d'utiliser plutot prtctl
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab      
-      !!----------------------------------------------------------------------
-      sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
-   END FUNCTION sum2d
-
-
-   SUBROUTINE trd_mld_trc_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mld_init  ***
-      !! 
-      !! ** Purpose :   computation of vertically integrated T and S budgets
-      !!      from ocean surface down to control surface (NetCDF output)
-      !!
-      !! ** Method/usage :
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER :: inum   ! logical unit
-      INTEGER :: ilseq, jl, jn, iiter
-      REAL(wp) ::   zjulian, zsto, zout
-      CHARACTER (LEN=40) ::   clop
-      CHARACTER (LEN=15) ::   csuff
-      CHARACTER (LEN=12) ::   clmxl
-      CHARACTER (LEN=16) ::   cltrcu
-      CHARACTER (LEN=10) ::   clvar
-
-      !!----------------------------------------------------------------------
-
-      ! ======================================================================
-      ! I. initialization
-      ! ======================================================================
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' trd_mld_trc_init : Mixed-layer trends for passive tracers                '
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
-         WRITE(numout,*)
-      ENDIF
-
-      
-      ! I.1 Check consistency of user defined preferences
-      ! -------------------------------------------------
-
-      IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
-         WRITE(numout,cform_err)
-         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
-         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
-         WRITE(numout,*) '                          you defined, nn_trd_trc   = ', nn_trd_trc
-         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
-         WRITE(numout,*) '                You should reconsider this choice.                        ' 
-         WRITE(numout,*) 
-         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
-         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
-         nstop = nstop + 1
-      ENDIF
-
-      ! * Debugging information *
-      IF( lldebug ) THEN
-         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
-         WRITE(numout,*) '               ln_trdmld_trc_instant = ', ln_trdmld_trc_instant
-      ENDIF
-
-      IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmld_trc_instant ) THEN
-         WRITE(numout,cform_err)
-         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
-         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
-         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
-         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
-         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
-         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
-         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
-         WRITE(numout,*) 
-         nstop = nstop + 1
-      ENDIF
-
-      ! I.2 Initialize arrays to zero or read a restart file
-      ! ----------------------------------------------------
-      nmoymltrd   = 0
-
-      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
-      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
-      tmlradm_trc        (:,:,:)   = 0.e0
-
-      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
-      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
-
-#if defined key_pisces_reduced
-      nmoymltrdbio   = 0
-      tmltrd_sum_bio     (:,:,:) = 0.e0     ;   tmltrd_csum_ln_bio (:,:,:) = 0.e0
-      DO jl = 1, jp_pisces_trd
-          ctrd_bio(jl,1) = ctrbil(jl)   ! long name
-          ctrd_bio(jl,2) = ctrbio(jl)   ! short name
-       ENDDO
-#endif
-
-      IF( ln_rsttr .AND. ln_trdmld_trc_restart ) THEN
-         CALL trd_mld_trc_rst_read
-      ELSE
-         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
-         tmlbn_trc          (:,:,:)   = 0.e0
-
-         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
-         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
-#if defined key_pisces_reduced
-         tmltrd_csum_ub_bio (:,:,:) = 0.e0
-#endif
-
-       ENDIF
-
-      icount = 1   ;   ionce  = 1  ! open specifier   
-
-#if defined key_pisces_reduced
-      icountbio = 1   ;   ioncebio  = 1  ! open specifier
-#endif
-
-      ! I.3 Read control surface from file ctlsurf_idx
-      ! ----------------------------------------------
-      IF( nn_ctls_trc == 1 ) THEN
-         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
-         READ ( inum ) nbol_trc
-         CLOSE( inum )
-      ENDIF
-
-      ! ======================================================================
-      ! II. netCDF output initialization
-      ! ======================================================================
-
-#if defined key_dimgout 
-      ???
-#else
-      ! clmxl = legend root for netCDF output
-      IF( nn_ctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
-         clmxl = 'Mixed Layer '
-      ELSE IF( nn_ctls_trc == 1 ) THEN                              ! control surface = read index from file 
-         clmxl = '      Bowl '
-      ELSE IF( nn_ctls_trc >= 2 ) THEN                              ! control surface = model level
-         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
-      ENDIF
-
-      ! II.1 Define frequency of output and means
-      ! -----------------------------------------
-      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
-      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cp time)
-      ENDIF
-#  if defined key_diainstant
-      IF( .NOT. ln_trdmld_trc_instant ) THEN
-         STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...'
-      ENDIF
-      zsto = nn_trd_trc * rdt
-      clop = "inst("//TRIM(clop)//")"
-#  else
-      IF( ln_trdmld_trc_instant ) THEN
-         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
-      ELSE
-         zsto = nn_trd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
-      ENDIF
-      clop = "ave("//TRIM(clop)//")"
-#  endif
-      zout = nn_trd_trc * rdt
-      iiter = ( nittrc000 - 1 ) / nn_dttrc
-
-      IF(lwp) WRITE (numout,*) '                netCDF initialization'
-
-      ! II.2 Compute julian date from starting date of the run
-      ! ------------------------------------------------------
-      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
-      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
-      IF(lwp) WRITE(numout,*)' '  
-      IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000               &
-           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
-           &   ,'Julian day : ', zjulian
-
-      ! II.3 Define the T grid trend file (nidtrd)
-      ! ------------------------------------------
-
-      !-- Define long and short names for the NetCDF output variables
-      !       ==> choose them according to trdmld_trc_oce.F90 <==
-
-      ctrd_trc(jpmld_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmld_trc_xad    ,2) = "_xad"
-      ctrd_trc(jpmld_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmld_trc_yad    ,2) = "_yad"
-      ctrd_trc(jpmld_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmld_trc_zad    ,2) = "_zad"
-      ctrd_trc(jpmld_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmld_trc_ldf    ,2) = "_ldf"
-      ctrd_trc(jpmld_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmld_trc_zdf    ,2) = "_zdf"
-      ctrd_trc(jpmld_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmld_trc_bbl    ,2) = "_bbl"
-      ctrd_trc(jpmld_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmld_trc_dmp    ,2) = "_dmp"
-      ctrd_trc(jpmld_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmld_trc_sbc    ,2) = "_sbc"
-      ctrd_trc(jpmld_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmld_trc_sms    ,2) = "_sms"
-      ctrd_trc(jpmld_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radb   ,2) = "_radb"
-      ctrd_trc(jpmld_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radn   ,2) = "_radn"
-      ctrd_trc(jpmld_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmld_trc_atf    ,2) = "_atf"
-
-      DO jn = 1, jptra      
-      !-- Create a NetCDF file and enter the define mode 
-         IF( ln_trdtrc(jn) ) THEN
-            csuff="ML_"//ctrcnm(jn)
-            CALL dia_nam( clhstnam, nn_trd_trc, csuff )
-            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
-               &        1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
-      
-            !-- Define the ML depth variable
-            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
-               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-
-         ENDIF
-      END DO
-
-#if defined key_pisces_reduced
-          !-- Create a NetCDF file and enter the define mode
-          CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
-          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
-             &             1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )
-#endif
-
-      !-- Define physical units
-      IF( rn_ucf_trc == 1. ) THEN
-         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit 
-      ELSEIF ( rn_ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
-         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf 
-      ELSE
-         cltrcu = "unknown?"
-      ENDIF
-
-      !-- Define miscellaneous passive tracer mixed-layer variables 
-      IF( jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb ) THEN
-         STOP 'Error : jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb'  ! see below
-      ENDIF
-
-      DO jn = 1, jptra
-         !
-         IF( ln_trdtrc(jn) ) THEN
-            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
-            CALL histdef(nidtrd(jn), trim(clvar),           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
-              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
-            CALL histdef(nidtrd(jn), trim(clvar)//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
-              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
-            CALL histdef(nidtrd(jn), trim(clvar)//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
-              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
-         
-            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmld_trc_atf
-               CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
-                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
-            END DO                                                                         ! if zsto=rdt above
-         
-            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_radb,1), & 
-              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
-         
-            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1),   & 
-              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
-         !
-         ENDIF
-      END DO
-
-#if defined key_pisces_reduced
-      DO jl = 1, jp_pisces_trd
-         CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1))   ,            &
-             &    cltrcu, jpi, jpj, nh_tb, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
-      END DO                                                                         ! if zsto=rdt above
-#endif
-
-      !-- Leave IOIPSL/NetCDF define mode
-      DO jn = 1, jptra
-         IF( ln_trdtrc(jn) )  CALL histend( nidtrd(jn), snc4set )
-      END DO
-
-#if defined key_pisces_reduced
-      !-- Leave IOIPSL/NetCDF define mode
-      CALL histend( nidtrdbio, snc4set )
-
-      IF(lwp) WRITE(numout,*)
-       IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
-#endif
-
-#endif        /* key_dimgout */
-   END SUBROUTINE trd_mld_trc_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Empty module
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trd_mld_trc( kt )                                   ! Empty routine
-      INTEGER, INTENT( in) ::   kt
-      WRITE(*,*) 'trd_mld_trc: You should not have seen this print! error?', kt
-   END SUBROUTINE trd_mld_trc
-   SUBROUTINE trd_mld_bio( kt )
-      INTEGER, INTENT( in) ::   kt
-      WRITE(*,*) 'trd_mld_bio: You should not have seen this print! error?', kt
-   END SUBROUTINE trd_mld_bio
-   SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
-      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
-      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
-      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmld            ! passive trc trend
-      WRITE(*,*) 'trd_mld_trc_zint: You should not have seen this print! error?', ptrc_trdmld(1,1,1)
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
-   END SUBROUTINE trd_mld_trc_zint
-   SUBROUTINE trd_mld_trc_init                                    ! Empty routine
-      WRITE(*,*) 'trd_mld_trc_init: You should not have seen this print! error?'
-   END SUBROUTINE trd_mld_trc_init
-#endif
-
-   !!======================================================================
-END MODULE trdmld_trc
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc_rst.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc_rst.F90	(revision 4624)
+++ 	(revision )
@@ -1,220 +1,0 @@
-MODULE trdmld_trc_rst
-   !!======================================================================
-   !!                       ***  MODULE  trdmld_rst  ***
-   !! Ocean dynamic :  Input/Output files for restart on mixed-layer diagnostics
-   !!======================================================================
-   !! History :  9.0  ! 07-03 (C. Deltel) Original code
-   !!----------------------------------------------------------------------
-  
-#if defined key_top && defined key_trdmld_trc
-   !!----------------------------------------------------------------------
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE trc             ! for nn_dttrc ctrcnm
-   USE trdmod_trc_oce  ! for lk_trdmld_trc
-
-   IMPLICIT NONE
-   PRIVATE
-  
-   PUBLIC   trd_mld_trc_rst_read    ! routine called by trd_mld_init
-   PUBLIC   trd_mld_trc_rst_write   ! routine called by step.F90
-  
-   INTEGER ::   nummldw_trc               ! logical unit for mld restart
-   !!---------------------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/TRD/trdmld_rst.F90,v 1.6 2006/11/14 09:46:13 opalod Exp $ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!---------------------------------------------------------------------------------
-  
-CONTAINS
-  
-
-    SUBROUTINE trd_mld_trc_rst_write( kt )
-      !!--------------------------------------------------------------------------------
-      !!                  ***  SUBROUTINE trd_mld_rst_wri  ***
-      !!                
-      !! ** Purpose :   Write mixed-layer diagnostics restart fields.
-      !!--------------------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      !
-      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
-      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
-      CHARACTER (len=35) :: charout
-      INTEGER :: jl,  jk, jn               ! loop indice
-      !!--------------------------------------------------------------------------------
-
-      IF( kt == nitrst - nn_dttrc .OR. nitend - nit000 + 1 < 2 * nn_dttrc ) THEN ! idem trcrst.F90
-         IF( nitrst > 1.0e9 ) THEN
-            WRITE(clkt,*) nitrst
-         ELSE
-           WRITE(clkt,'(i8.8)') nitrst
-         ENDIF
-         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trdrst_trc_out)
-         IF(lwp) WRITE(numout,*) '             open ocean restart_mld_trc NetCDF  '//clname
-         CALL iom_open( clname, nummldw_trc, ldwrt = .TRUE., kiolib = jprstlib )
-      ENDIF
-
-      IF( kt == nitend .AND. lk_trdmld_trc ) THEN
-
-         IF( kt == nitend .AND. lwp ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'trdmld_trc_rst: output for ML diags. restart, with trd_mld_trc_rst_write routine'
-            WRITE(numout,*) '~~~~~~~~~~~~~~'
-            WRITE(numout,*)
-         ENDIF
-
-         IF( ln_trdmld_trc_instant ) THEN 
-            !
-            DO jn = 1, jptra
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
-            END DO
-            !
-         ELSE
-            !
-            CALL iom_rstput( kt, nitrst, nummldw_trc, 'rmldbn_trc', rmldbn_trc )  ! 2D x 1
-            
-            !                                                          ! ===========
-            DO jn = 1, jptra                                           ! tracer loop
-               !                                                       ! ===========
-
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc (:,:,jn) ) ! 2D x jptra
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc   (:,:,jn) ) ! 2D x jptra
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc (:,:,jn) ) ! 2D x jptra
-
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! 2D x jptra
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tml_sumb_trc_'//ctrcnm(jn), tml_sumb_trc(:,:,jn) ) ! 2D x jptra
-               
-               DO jk = 1, jpltrd_trc
-                  IF( jk < 10 )   THEN
-                     WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I1)") ctrcnm(jn), jk
-                  ELSE
-                     WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I2)") ctrcnm(jn), jk
-                  ENDIF
-                  CALL iom_rstput( kt, nitrst, nummldw_trc, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
-               END DO
-               
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
-                    &           tmltrd_atf_sumb_trc(:,:,jn) ) ! 2D x jptra
-
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
-                    &           tmltrd_rad_sumb_trc(:,:,jn) ) ! 2D x jptra
-               !                                                       ! ===========
-            END DO                                                     ! tracer loop
-            !                                                          ! ===========
-#if defined key_pisces_reduced
-            DO jl = 1, jp_pisces_trd
-               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmltrd_csum_ub_bio'//ctrd_bio(jl,2), tmltrd_csum_ub_bio(:,:,jl) )
-            ENDDO
-#endif
-
-         ENDIF
-         
-         CALL iom_close( nummldw_trc )
-         lrst_trc = .TRUE.
-
-      ENDIF
-
-    END SUBROUTINE trd_mld_trc_rst_write
-
-
-    SUBROUTINE trd_mld_trc_rst_read
-      !!----------------------------------------------------------------------------
-      !!                   ***  SUBROUTINE trd_mld_rst_lec  ***
-      !!                   
-      !! ** Purpose :   Read file for mixed-layer diagnostics restart.
-      !!----------------------------------------------------------------------------
-      INTEGER  ::  inum       ! temporary logical unit
-      !
-      CHARACTER (len=35) :: charout
-      INTEGER ::  jk, jn, jl     ! loop indice
-      INTEGER ::  jlibalt = jprstlib
-      LOGICAL ::  llok
-      !!-----------------------------------------------------------------------------
-      
-      IF(lwp)  THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' trd_mld_trc_rst_read : read the NetCDF MLD restart file'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
-      ENDIF
-      
-      IF ( jprstlib == jprstdimg ) THEN
-        ! eventually read netcdf file (monobloc)  for restarting on different number of processors
-        ! if {cn_trdrst_trc_in}.nc exists, then set jlibalt to jpnf90
-        INQUIRE( FILE = TRIM(cn_trdrst_trc_in)//'.nc', EXIST = llok )
-        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
-      ENDIF
-
-      CALL iom_open( cn_trdrst_trc_in, inum, kiolib = jlibalt ) 
-      
-      IF( ln_trdmld_trc_instant ) THEN 
-         
-         DO jn = 1, jptra
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
-         END DO
-         
-      ELSE
-         CALL iom_get( inum, jpdom_autoglo, 'rmldbn_trc', rmldbn_trc ) ! needed for rmld_sum
-         
-         !                                                          ! ===========
-         DO jn = 1, jptra                                           ! tracer loop
-            !                                                       ! ===========
-            CALL iom_get( inum, jpdom_autoglo, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc(:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc(:,:,jn) )
-
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! needed for tml_sum
-            CALL iom_get( inum, jpdom_autoglo, 'tml_sumb_trc_'//ctrcnm(jn), tml_sumb_trc(:,:,jn) )
-            
-            DO jk = 1, jpltrd_trc
-               IF( jk < 10 )   THEN
-                  WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I1)") ctrcnm(jn), jk
-               ELSE
-                  WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I2)") ctrcnm(jn), jk
-               ENDIF
-               CALL iom_get( inum, jpdom_autoglo, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
-            END DO
-            
-            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
-                 &        tmltrd_atf_sumb_trc(:,:,jn) )
-
-            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
-                 &        tmltrd_rad_sumb_trc(:,:,jn) )
-            !                                                       ! ===========
-         END DO                                                     ! tracer loop
-         !                                                          ! ===========
-
-#if defined key_pisces_reduced
-         DO jl = 1, jp_pisces_trd
-            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_csum_ub_bio'//ctrd_bio(jl,2), tmltrd_csum_ub_bio(:,:,jl) )
-         ENDDO
-#endif
-         
-         CALL iom_close( inum )
-      ENDIF
-      
-    END SUBROUTINE trd_mld_trc_rst_read
-  
-#else
-  !!=================================================================================
-  !!                       ***  MODULE  trdmld_rst  ***
-  !! Ocean dynamic :  Input/Output files for restart on mixed-layer diagnostics
-  !!=================================================================================
-CONTAINS
-  SUBROUTINE trd_mld_trc_rst_opn( kt )
-    WRITE(*,*) 'trd_mld_trc_rst_opn: You should not have seen this print! error?', kt
-  END SUBROUTINE trd_mld_trc_rst_opn
-  SUBROUTINE trd_mld_trc_rst_write( kt )           !  No ML diags ==> empty routine
-    WRITE(*,*) 'trd_mld_trc_rst_wri: You should not have seen this print! error?', kt
-  END SUBROUTINE trd_mld_trc_rst_write
-  SUBROUTINE trd_mld_trc_rst_read                  !  No ML Diags ==> empty routine
-  END SUBROUTINE trd_mld_trc_rst_read
-#endif
-
-  !!=================================================================================
-END MODULE trdmld_trc_rst
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmod_trc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmod_trc.F90	(revision 4624)
+++ 	(revision )
@@ -1,158 +1,0 @@
-MODULE trdmod_trc
-   !!======================================================================
-   !!                       ***  MODULE  trdmod_trc  ***
-   !! Ocean diagnostics:  mixed layer passive tracer trends 
-   !!======================================================================
-   !! History :  3.0  !  2010-07  (C. Ethe)  Original code (from trdmod.F90)
-   !!----------------------------------------------------------------------
-#if   defined key_top && ( defined key_trdmld_trc   ||   defined key_trdtrc )
-   !!----------------------------------------------------------------------
-   !!   'key_trdmld_trc'                  mixed layer trend diagnostics
-   !!   'key_trdtrc'                      3D trend diagnostics
-   !!----------------------------------------------------------------------
-   !!   trd_mod_trc      : passive tracer trends 
-   !!----------------------------------------------------------------------
-   USE trc               ! tracer definitions (trn, trb, tra, etc.)
-   USE trcnam_trp
-   USE trdmod_oce
-   USE trdmod_trc_oce    ! definition of main arrays used for trends computations
-   USE trdmld_trc        ! Mixed layer trends diag.
-   USE iom               ! I/O library
-
-   IMPLICIT NONE
-   PRIVATE
-
-   INTERFACE trd_mod_trc
-      MODULE PROCEDURE trd_mod_trc_trp, trd_mod_trc_bio
-   END INTERFACE
-
-   PUBLIC trd_mod_trc
-
-   !! * Substitutions
-#  include "top_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header:  $ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trd_mod_trc_trp( ptrtrd, kjn, ktrd, kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mod_trc  ***
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )  ::   kt                                  ! time step
-      INTEGER, INTENT( in )  ::   kjn                                 ! tracer index
-      INTEGER, INTENT( in )  ::   ktrd                                ! tracer trend index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
-      CHARACTER (len=20) :: cltra
-      !!----------------------------------------------------------------------
-
-      IF( kt == nittrc000 ) THEN
-!         IF(lwp)WRITE(numout,*)
-!         IF(lwp)WRITE(numout,*) 'trd_mod_trc:'
-!         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~'
-      ENDIF
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Mixed layer trends for passive tracers
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-#if defined key_trdmld_trc  
-      IF( lk_trdmld_trc .AND. ln_trdtrc( kjn ) ) THEN
-         !
-         SELECT CASE ( ktrd )
-         CASE ( jptra_trd_xad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xad, '3D', kjn )
-         CASE ( jptra_trd_yad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yad, '3D', kjn )
-         CASE ( jptra_trd_zad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zad, '3D', kjn )
-         CASE ( jptra_trd_ldf     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf, '3D', kjn )
-         CASE ( jptra_trd_bbl     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbl, '3D', kjn )
-         CASE ( jptra_trd_zdf     )
-            IF( ln_trcldf_iso ) THEN
-               CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf, '3D', kjn )
-            ELSE
-               CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zdf, '3D', kjn )
-            ENDIF
-         CASE ( jptra_trd_dmp     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_dmp , '3D', kjn )
-         CASE ( jptra_trd_nsr     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sbc , '2D', kjn )
-         CASE ( jptra_trd_sms     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms , '3D', kjn )
-         CASE ( jptra_trd_radb    )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radb, '3D', kjn )
-         CASE ( jptra_trd_radn    )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radn, '3D', kjn )
-         CASE ( jptra_trd_atf     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_atf , '3D', kjn )
-         END SELECT
-         !
-      END IF
-#endif
-
-      IF( lk_trdtrc .AND. ln_trdtrc( kjn ) ) THEN
-         !
-         SELECT CASE( ktrd )
-         CASE( jptra_trd_xad  )       ;    WRITE (cltra,'("XAD_",4a)')
-         CASE( jptra_trd_yad  )       ;    WRITE (cltra,'("YAD_",4a)')
-         CASE( jptra_trd_zad  )       ;    WRITE (cltra,'("ZAD_",4a)')
-         CASE( jptra_trd_ldf  )       ;    WRITE (cltra,'("LDF_",4a)')
-         CASE( jptra_trd_bbl  )       ;    WRITE (cltra,'("BBL_",4a)')
-         CASE( jptra_trd_nsr  )       ;    WRITE (cltra,'("FOR_",4a)')
-         CASE( jptra_trd_zdf  )       ;    WRITE (cltra,'("ZDF_",4a)')
-         CASE( jptra_trd_dmp  )       ;    WRITE (cltra,'("DMP_",4a)')
-         CASE( jptra_trd_sms  )       ;    WRITE (cltra,'("SMS_",4a)')
-         CASE( jptra_trd_atf  )       ;    WRITE (cltra,'("ATF_",4a)')
-         CASE( jptra_trd_radb )       ;    WRITE (cltra,'("RDB_",4a)')
-         CASE( jptra_trd_radn )       ;    WRITE (cltra,'("RDN_",4a)')
-         END SELECT
-                                          cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
-                                          CALL iom_put( cltra,  ptrtrd(:,:,:) )
-         !
-      END IF
-
-   END SUBROUTINE trd_mod_trc_trp
-
-   SUBROUTINE trd_mod_trc_bio( ptrbio, ktrd, kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_mod_bio  ***
-      !!----------------------------------------------------------------------
-
-      INTEGER, INTENT( in )  ::   kt                                  ! time step
-      INTEGER, INTENT( in )  ::   ktrd                                ! bio trend index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout )  ::   ptrbio  ! Bio trend
-      !!----------------------------------------------------------------------
-
-#if defined key_trdmld_trc  
-      CALL trd_mld_bio_zint( ptrbio, ktrd ) ! Verticaly integrated biological trends
-#endif
-
-   END SUBROUTINE trd_mod_trc_bio
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Empty module
-   !!----------------------------------------------------------------------
-
-   INTERFACE trd_mod_trc
-      MODULE PROCEDURE trd_mod_trc_trp, trd_mod_trc_bio
-   END INTERFACE
-
-CONTAINS
-
-   SUBROUTINE trd_mod_trc_trp( ptrtrd, kjn, ktrd, kt )
-      INTEGER               , INTENT( in )     ::   kt      ! time step
-      INTEGER               , INTENT( in )     ::   kjn     ! tracer index
-      INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
-      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
-      WRITE(*,*) 'trd_mod_trc_trp : You should not have seen this print! error?', ptrtrd(1,1,1)
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
-   END SUBROUTINE trd_mod_trc_trp
-
-   SUBROUTINE trd_mod_trc_bio( ptrbio, ktrd, kt )
-      INTEGER               , INTENT( in )     ::   kt      ! time step
-      INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
-      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrbio  ! Temperature or U trend
-      WRITE(*,*) 'trd_mod_trc_trp : You should not have seen this print! error?', ptrbio(1,1,1)
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
-   END SUBROUTINE trd_mod_trc_bio
-
-#endif
-   !!======================================================================
-END MODULE trdmod_trc
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmod_trc_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmod_trc_oce.F90	(revision 4624)
+++ 	(revision )
@@ -1,180 +1,0 @@
-MODULE trdmod_trc_oce
-   !!======================================================================
-   !!                   ***  MODULE trdmod_trc_oce  ***
-   !! Ocean trends :   set tracer and momentum trend variables
-   !!======================================================================
-#if defined key_top   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_top'                                                TOP models
-   !!----------------------------------------------------------------------
-   USE par_oce       ! ocean parameters
-   USE par_trc       ! passive tracers parameters
-
-   IMPLICIT NONE
-   PUBLIC
-
-   !                                         !!* Namelist namtoptrd:  diagnostics on passive tracers trends
-   INTEGER  ::    nn_trd_trc                  !: time step frequency dynamics and tracers trends
-   INTEGER  ::    nn_ctls_trc                 !: control surface type for trends vertical integration
-   REAL(wp) ::    rn_ucf_trc                  !: unit conversion factor (for netCDF trends outputs)
-   LOGICAL  ::    ln_trdmld_trc_instant       !: flag to diagnose inst./mean ML trc trends
-   LOGICAL  ::    ln_trdmld_trc_restart       !: flag to restart mixed-layer trc diagnostics
-   CHARACTER(len=50) ::  cn_trdrst_trc_in     !: suffix of pass. tracer restart name (input)
-   CHARACTER(len=50) ::  cn_trdrst_trc_out    !: suffix of pass. tracer restart name (output)
-   LOGICAL, DIMENSION(jptra) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist)
-
-# if defined key_trdtrc && defined key_iomput
-   LOGICAL, PARAMETER ::   lk_trdtrc = .TRUE. 
-# else
-   LOGICAL, PARAMETER ::   lk_trdtrc = .FALSE.   !: ML trend flag
-# endif
-
-# if defined key_trdmld_trc   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_trdmld_trc'                     mixed layer trends diagnostics
-   !!----------------------------------------------------------------------
-
-   LOGICAL, PARAMETER ::   lk_trdmld_trc = .TRUE.    !: ML trend flag
-
-   INTEGER, PARAMETER ::            & !: mixed layer trends indices
-        jpmld_trc_xad     =  1,     & !:     zonal       advection     
-        jpmld_trc_yad     =  2,     & !:     meridonal   =========
-        jpmld_trc_zad     =  3,     & !:     vertical    =========
-        jpmld_trc_ldf     =  4,     & !:     lateral diffusion (geopot. or iso-neutral)
-        jpmld_trc_zdf     =  5,     & !:     vertical diffusion (TKE)
-        jpmld_trc_bbl     =  6,     & !:     bottom boundary layer (advective/diffusive)
-        jpmld_trc_dmp     =  7,     & !:     internal restoring trend
-        jpmld_trc_sbc     =  8,     & !:     forcing 
-        jpmld_trc_sms     =  9,     & !:     sources minus sinks trend
-  !     jpmld_trc_xxx     = xx,     & !:     add here any additional trend    (** AND UPDATE JPLTRD_TRC BELOW **)
-        jpmld_trc_radn    = 10,     & !:     corr. trn<0 in trcrad
-        jpmld_trc_radb    = 11,     & !:     corr. trb<0 in trcrad (like atf) (** MUST BE BEFORE THE LAST ONE **)
-        jpmld_trc_atf     = 12        !:     asselin trend                    (** MUST BE    THE      LAST ONE**)
-
-   !! Trends diagnostics parameters
-   !!---------------------------------------------------------------------
-   INTEGER, PARAMETER :: jpltrd_trc = 12    !: number of mixed-layer trends arrays
-      
-   INTEGER            :: jpktrd_trc         !: max level for mixed-layer trends diag.
-
-   !! Arrays used for diagnosing mixed-layer trends 
-   !!---------------------------------------------------------------------
-   CHARACTER(LEN=80) :: clname_trc, ctrd_trc(jpltrd_trc+1,2)
-
-   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &
-      nmld_trc       , &                            !: mixed layer depth indexes 
-      nbol_trc                                   !: mixed-layer depth indexes when read from file
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  wkx_trc  !:
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  rmld_trc     !: ML depth (m) corresponding to nmld_trc
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  rmld_sum_trc !: needed to compute the leap-frog time mean of ML depth
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  rmldbn_trc   !: idem
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &
-      tml_trc    ,                        &      !: \ "now" mixed layer temperature/salinity
-      tmlb_trc   ,                        &      !: /  and associated "before" fields
-      tmlbb_trc  ,                        &      !: \  idem, but valid at the 1rst time step of the
-      tmlbn_trc  ,                        &      !: /  current analysis window
-      tml_sum_trc,                        &      !: mixed layer T, summed over the current analysis period
-      tml_sumb_trc,                       &      !: idem, but from the previous analysis period
-      tmltrd_atf_sumb_trc,                &      !: Asselin trends, summed over the previous analysis period
-      tmltrd_rad_sumb_trc                        !: trends due to trb correction in trcrad.F90, summed over the
-                                                 !:     previous analysis period
-                                                 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &      
-      tmlatfb_trc, tmlatfn_trc ,          &      !: "before" Asselin contrib. at beginning of the averaging
-                                                 !:     period (i.e. last contrib. from previous such period)
-                                                 !:     and "now" Asselin contrib. to the ML trc. trends
-      tmlatfm_trc,                        &      !: accumulator for Asselin trends (needed for storage only)
-      tmlradb_trc, tmlradn_trc ,          &      !: similar to Asselin above, but for the trend due to trb
-                                                 !:     correction in trcrad.F90
-      tmlradm_trc                                !: accumulator for the previous trcrad trend
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  &
-      tmltrd_trc,                         &      !: \ physical contributions to the total trend (for T/S),
-                                                 !: / cumulated over the current analysis window
-      tmltrd_sum_trc,                     &      !: sum of these trends over the analysis period
-      tmltrd_csum_ln_trc,                 &      !: now cumulated sum of trends over the "lower triangle"
-      tmltrd_csum_ub_trc                         !: before (prev. analysis period) cumulated sum over the
-                                                 !: upper triangle
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &
-      tmltrdm_trc                                !: total cumulative trends over the analysis window
-
-# else
-   LOGICAL, PARAMETER ::   lk_trdmld_trc = .FALSE.   !: ML trend flag
-# endif
-
-# if defined key_pisces_reduced
-   CHARACTER(LEN=80) :: clname_bio, ctrd_bio(jpdiabio,2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &
-      tmltrd_bio,                         &      !: \ biological contributions to the total trend ,
-                                                 !: / cumulated over the current analysis window
-      tmltrd_sum_bio,                     &      !: sum of these trends over the analysis period
-      tmltrd_csum_ln_bio,                 &      !: now cumulated sum of trends over the "lower triangle"
-      tmltrd_csum_ub_bio                         !: before (prev. analysis period) cumulated sum over the
-                                                 !: upper triangle
-#endif
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/TRD/trdmld_oce.F90,v 1.2 2005/03/27 18:35:23 opalod Exp $ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION trd_mod_trc_oce_alloc()
-      !!----------------------------------------------------------------------
-      !!         *** ROUTINE trd_mod_trc_oce_alloc ***
-      !!----------------------------------------------------------------------
-      USE lib_mpp, ONLY: ctl_warn
-      INTEGER :: ierr(2)
-      !!----------------------------------------------------------------------
-      ierr(:) = 0
-      !
-# if defined key_trdmld_trc
-      ALLOCATE(nmld_trc(jpi,jpj),          nbol_trc(jpi,jpj),           &
-               wkx_trc(jpi,jpj,jpk),       rmld_trc(jpi,jpj),           &
-               rmld_sum_trc(jpi,jpj),      rmldbn_trc(jpi,jpj),         &
-               tml_trc(jpi,jpj,jptra),     tmlb_trc(jpi,jpj,jptra),     &
-               tmlbb_trc(jpi,jpj,jptra),   tmlbn_trc(jpi,jpj,jptra),    &
-               tml_sum_trc(jpi,jpj,jptra), tml_sumb_trc(jpi,jpj,jptra), &
-               tmltrd_atf_sumb_trc(jpi,jpj,jptra),                      &
-               tmltrd_rad_sumb_trc(jpi,jpj,jptra),                      &
-               !
-               tmlatfb_trc(jpi,jpj,jptra), tmlatfn_trc(jpi,jpj,jptra),  &
-               tmlatfm_trc(jpi,jpj,jptra), tmlradb_trc(jpi,jpj,jptra),  &
-               tmlradn_trc(jpi,jpj,jptra), tmlradm_trc(jpi,jpj,jptra),  &
-               !
-               tmltrd_trc(jpi,jpj,jpltrd_trc,jptra)         , &
-               tmltrd_sum_trc(jpi,jpj,jpltrd_trc,jptra)     , &
-               tmltrd_csum_ln_trc(jpi,jpj,jpltrd_trc,jptra) , &
-               tmltrd_csum_ub_trc(jpi,jpj,jpltrd_trc,jptra) , &
-               !
-               tmltrdm_trc(jpi,jpj,jptra)                   , STAT=ierr(1) )
-#endif
-      !
-# if defined key_pisces_reduced
-      ALLOCATE( tmltrd_bio        (jpi,jpj,jpdiabio) ,     &
-         &      tmltrd_sum_bio    (jpi,jpj,jpdiabio) ,     &
-         &      tmltrd_csum_ln_bio(jpi,jpj,jpdiabio) ,     &
-         &      tmltrd_csum_ub_bio(jpi,jpj,jpdiabio) , STAT=ierr(2) )
-# endif
-      !
-      trd_mod_trc_oce_alloc = MAXVAL(ierr)
-      !
-      IF( trd_mod_trc_oce_alloc /= 0 )   CALL ctl_warn('trd_mod_trc_oce_alloc: failed to allocate arrays')
-      !
-# if defined key_trdmld_trc
-      jpktrd_trc = jpk      ! Initialise what used to be a parameter - max level for mixed-layer trends diag.
-# endif
-      !
-   END FUNCTION trd_mod_trc_oce_alloc
-
-#else
-   !!----------------------------------------------------------------------
-   !!  Empty module :                                     No passive tracer
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE trdmod_trc_oce
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90	(revision 4990)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90	(revision 4990)
@@ -0,0 +1,1437 @@
+MODULE trdmxl_trc
+   !!======================================================================
+   !!                       ***  MODULE  trdmxl_trc  ***
+   !! Ocean diagnostics:  mixed layer passive tracer trends 
+   !!======================================================================
+   !! History :  9.0  !  06-08  (C. Deltel)  Original code (from trdmxl.F90)
+   !!                 !  07-04  (C. Deltel)  Bug fix : add trcrad trends
+   !!                 !  07-06  (C. Deltel)  key_gyre : do not call lbc_lnk
+   !!----------------------------------------------------------------------
+#if   defined key_top && ( defined key_trdmxl_trc   ||   defined key_esopa )
+   !!----------------------------------------------------------------------
+   !!   'key_trdmxl_trc'                      mixed layer trend diagnostics
+   !!----------------------------------------------------------------------
+   !!   trd_mxl_trc      : passive tracer cumulated trends averaged over ML
+   !!   trd_mxl_trc_zint : passive tracer trends vertical integration
+   !!   trd_mxl_trc_init : initialization step
+   !!----------------------------------------------------------------------
+   USE trc               ! tracer definitions (trn, trb, tra, etc.)
+   USE trc_oce, ONLY :   nn_dttrc  ! frequency of step on passive tracers
+   USE dom_oce           ! domain definition
+   USE zdfmxl  , ONLY : nmln ! number of level in the mixed layer
+   USE zdf_oce , ONLY : avt  ! vert. diffusivity coef. at w-point for temp  
+# if defined key_zdfddm   
+   USE zdfddm  , ONLY : avs  ! salinity vertical diffusivity coeff. at w-point
+# endif
+   USE trcnam_trp        ! passive tracers transport namelist variables
+   USE trdtrc_oce    ! definition of main arrays used for trends computations
+   USE in_out_manager    ! I/O manager
+   USE dianam            ! build the name of file (routine)
+   USE ldfslp            ! iso-neutral slopes 
+   USE ioipsl            ! NetCDF library
+   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp           ! MPP library
+   USE trdmxl_trc_rst    ! restart for diagnosing the ML trends
+   USE prtctl            ! print control
+   USE sms_pisces        ! PISCES bio-model
+   USE wrk_nemo          ! Memory allocation
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trd_mxl_trc
+   PUBLIC trd_mxl_trc_alloc
+   PUBLIC trd_mxl_bio
+   PUBLIC trd_mxl_trc_init
+   PUBLIC trd_mxl_trc_zint
+   PUBLIC trd_mxl_bio_zint
+
+   CHARACTER (LEN=40) ::  clhstnam                                ! name of the trends NetCDF file
+   INTEGER ::   nmoymltrd
+   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::   ndextrd1
+   INTEGER, DIMENSION(jptra) ::   nidtrd, nh_t
+   INTEGER ::   ndimtrd1                        
+   INTEGER, SAVE ::  ionce, icount
+#if defined key_pisces_reduced
+   INTEGER ::   nidtrdbio, nh_tb
+   INTEGER, SAVE ::  ioncebio, icountbio
+   INTEGER, SAVE ::   nmoymltrdbio
+#endif
+   LOGICAL :: llwarn  = .TRUE.                                    ! this should always be .TRUE.
+   LOGICAL :: lldebug = .TRUE.
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztmltrd2   !
+#if defined key_pisces_reduced
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ztmltrdbio2  ! only needed for mean diagnostics in trd_mxl_bio()
+#endif
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+#  include "zdfddm_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Header:  $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION trd_mxl_trc_alloc()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mxl_trc_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) ,      &
+#if defined key_pisces_reduced
+         &      ztmltrdbio2(jpi,jpj,jpdiabio)      ,      &
+#endif
+         &      ndextrd1(jpi*jpj)                  ,  STAT=trd_mxl_trc_alloc)
+         !
+      IF( lk_mpp                )   CALL mpp_sum ( trd_mxl_trc_alloc )
+      IF( trd_mxl_trc_alloc /=0 )   CALL ctl_warn('trd_mxl_trc_alloc: failed to allocate arrays')
+      !
+   END FUNCTION trd_mxl_trc_alloc
+
+
+   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mxl_trc_zint  ***
+      !! 
+      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments 
+      !!                to the subroutine. This vertical average is performed from ocean
+      !!                surface down to a chosen control surface.
+      !!
+      !! ** Method/usage :
+      !!      The control surface can be either a mixed layer depth (time varying)
+      !!      or a fixed surface (jk level or bowl). 
+      !!      Choose control surface with nctls_trc in namelist NAMTRD :
+      !!        nctls_trc = -2  : use isopycnal surface
+      !!        nctls_trc = -1  : use euphotic layer with light criterion 
+      !!        nctls_trc =  0  : use mixed layer with density criterion 
+      !!        nctls_trc =  1  : read index from file 'ctlsurf_idx'
+      !!        nctls_trc >  1  : use fixed level surface jk = nctls_trc
+      !!      Note: in the remainder of the routine, the volume between the 
+      !!            surface and the control surface is called "mixed-layer"
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, INTENT( in ) ::   ktrd, kjn                        ! ocean trend index and passive tracer rank
+      CHARACTER(len=2), INTENT( in ) ::  ctype                    ! surface/bottom (2D) or interior (3D) physics
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  ptrc_trdmxl ! passive tracer trend
+      !
+      INTEGER ::   ji, jj, jk, isum
+      REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, zvlmsk )
+
+      ! I. Definition of control surface and integration weights
+      ! --------------------------------------------------------
+
+      ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
+         !
+         tmltrd_trc(:,:,:,:) = 0.e0                               ! <<< reset trend arrays to zero
+          
+         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
+         SELECT CASE ( nn_ctls_trc )                                ! choice of the control surface
+            CASE ( -2  )   ;   STOP 'trdmxl_trc : not ready '     !     -> isopycnal surface (see ???)
+#if defined key_pisces || defined key_pisces_reduced
+            CASE ( -1  )   ;   nmld_trc(:,:) = neln(:,:)          !     -> euphotic layer with light criterion
+#endif
+            CASE (  0  )   ;   nmld_trc(:,:) = nmln(:,:)          !     -> ML with density criterion (see zdfmxl)
+            CASE (  1  )   ;   nmld_trc(:,:) = nbol_trc(:,:)          !     -> read index from file
+            CASE (  2: )   ;   nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
+                               nmld_trc(:,:) = nn_ctls_trc + 1      !     -> model level
+         END SELECT
+
+         ! ... Compute ndextrd1 and ndimtrd1  ??? role de jpktrd_trc
+         IF( ionce == 1 ) THEN
+            !
+            isum  = 0   ;   zvlmsk(:,:) = 0.e0
+
+            IF( jpktrd_trc < jpk ) THEN                           ! description ???
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
+                        zvlmsk(ji,jj) = tmask(ji,jj,1)
+                     ELSE
+                        isum = isum + 1
+                        zvlmsk(ji,jj) = 0.e0
+                     ENDIF
+                  END DO
+               END DO
+            ENDIF
+
+            IF( isum > 0 ) THEN                                   ! index of ocean points (2D only)
+               WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum 
+               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
+            ELSE 
+               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
+            ENDIF                                
+
+            ionce = 0                                             ! no more pass here
+            !
+         ENDIF   ! ionce == 1
+         
+         ! ... Weights for vertical averaging
+         wkx_trc(:,:,:) = 0.e0
+         DO jk = 1, jpktrd_trc                                    ! initialize wkx_trc with vertical scale factor in mixed-layer
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         
+         rmld_trc(:,:) = 0.e0
+         DO jk = 1, jpktrd_trc                                    ! compute mixed-layer depth : rmld_trc
+            rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
+         END DO
+         
+         DO jk = 1, jpktrd_trc                                    ! compute integration weights
+            wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
+         END DO
+
+         icount = 0                                               ! <<< flag = off : control surface & integr. weights
+         !                                                        !     computed only once per time step
+      ENDIF ONCE_PER_TIME_STEP
+
+      ! II. Vertical integration of trends in the mixed-layer
+      ! -----------------------------------------------------
+
+      SELECT CASE ( ctype )
+         CASE ( '3D' )                                            ! mean passive tracer trends in the mixed-layer
+            DO jk = 1, jpktrd_trc
+               tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)   
+            END DO
+         CASE ( '2D' )                                            ! forcing at upper boundary of the mixed-layer 
+            tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,1) * wkx_trc(:,:,1)  ! non penetrative
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi, jpj, zvlmsk )
+      !
+   END SUBROUTINE trd_mxl_trc_zint
+
+
+   SUBROUTINE trd_mxl_bio_zint( ptrc_trdmxl, ktrd )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mxl_bio_zint  ***
+      !!
+      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments
+      !!                to the subroutine. This vertical average is performed from ocean
+      !!                surface down to a chosen control surface.
+      !!
+      !! ** Method/usage :
+      !!      The control surface can be either a mixed layer depth (time varying)
+      !!      or a fixed surface (jk level or bowl).
+      !!      Choose control surface with nctls in namelist NAMTRD :
+      !!        nctls_trc = 0  : use mixed layer with density criterion
+      !!        nctls_trc = 1  : read index from file 'ctlsurf_idx'
+      !!        nctls_trc > 1  : use fixed level surface jk = nctls_trc
+      !!      Note: in the remainder of the routine, the volume between the
+      !!            surface and the control surface is called "mixed-layer"
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER                         , INTENT(in) ::   ktrd          ! bio trend index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   ptrc_trdmxl   ! passive trc trend
+#if defined key_pisces_reduced
+      !
+      INTEGER ::   ji, jj, jk, isum
+      REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, zvlmsk )
+
+      ! I. Definition of control surface and integration weights
+      ! --------------------------------------------------------
+      !            ==> only once per time step <==
+
+      IF( icountbio == 1 ) THEN
+         !
+         tmltrd_bio(:,:,:) = 0.e0    ! <<< reset trend arrays to zero
+         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
+         SELECT CASE ( nn_ctls_trc )                                    ! choice of the control surface
+            CASE ( -2  )   ;   STOP 'trdmxl_trc : not ready '     !     -> isopycnal surface (see ???)
+            CASE ( -1  )   ;   nmld_trc(:,:) = neln(:,:)          !     -> euphotic layer with light criterion
+            CASE (  0  )   ;   nmld_trc(:,:) = nmln(:,:)          !     -> ML with density criterion (see zdfmxl)
+            CASE (  1  )   ;   nmld_trc(:,:) = nbol_trc(:,:)          !     -> read index from file
+            CASE (  2: )   ;   nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
+                               nmld_trc(:,:) = nn_ctls_trc + 1          !     -> model level
+         END SELECT
+
+         ! ... Compute ndextrd1 and ndimtrd1 only once
+         IF( ioncebio == 1 ) THEN
+            !
+            ! Check of validity : nmld_trc(ji,jj) <= jpktrd_trc
+            isum        = 0
+            zvlmsk(:,:) = 0.e0
+
+            IF( jpktrd_trc < jpk ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
+                        zvlmsk(ji,jj) = tmask(ji,jj,1)
+                     ELSE
+                        isum = isum + 1
+                        zvlmsk(ji,jj) = 0.
+                     END IF
+                  END DO
+               END DO
+            END IF
+
+            ! Index of ocean points (2D only)
+            IF( isum > 0 ) THEN
+               WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
+               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
+            ELSE
+               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
+            END IF
+
+            ioncebio = 0                  ! no more pass here
+            !
+         END IF !  ( ioncebio == 1 )
+
+         ! ... Weights for vertical averaging
+         wkx_trc(:,:,:) = 0.e0
+         DO jk = 1, jpktrd_trc         ! initialize wkx_trc with vertical scale factor in mixed-layer
+            DO jj = 1,jpj
+              DO ji = 1,jpi
+                  IF( jk - nmld_trc(ji,jj) < 0. )   wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+
+         rmld_trc(:,:) = 0.
+         DO jk = 1, jpktrd_trc         ! compute mixed-layer depth : rmld_trc
+            rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
+         END DO
+
+         DO jk = 1, jpktrd_trc         ! compute integration weights
+            wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
+         END DO
+
+         icountbio = 0                    ! <<< flag = off : control surface & integr. weights
+         !                             !     computed only once per time step
+      END IF ! ( icountbio == 1 )
+
+      ! II. Vertical integration of trends in the mixed-layer
+      ! -----------------------------------------------------
+
+
+      DO jk = 1, jpktrd_trc
+         tmltrd_bio(:,:,ktrd) = tmltrd_bio(:,:,ktrd) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)
+      END DO
+
+      CALL wrk_dealloc( jpi, jpj, zvlmsk )
+#endif
+      !
+   END SUBROUTINE trd_mxl_bio_zint
+
+
+   SUBROUTINE trd_mxl_trc( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mxl_trc  ***
+      !! 
+      !! ** Purpose :  Compute and cumulate the mixed layer trends over an analysis
+      !!               period, and write NetCDF (or dimg) outputs.
+      !!
+      !! ** Method/usage :
+      !!          The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant 
+      !!          logical namelist variable) :
+      !!          1) to explain the difference between initial and final 
+      !!             mixed-layer T & S (where initial and final relate to the
+      !!             current analysis window, defined by ntrc_trc in the namelist)
+      !!          2) to explain the difference between the current and previous 
+      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
+      !!             performed over each analysis window).
+      !!
+      !! ** Consistency check : 
+      !!        If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
+      !!        entrainment) should be zero, at machine accuracy. Note that in the case
+      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
+      !!        over the first two analysis windows (except if restart).
+      !!        N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
+      !!             for checking residuals.
+      !!             On a NEC-SX5 computer, this typically leads to:
+      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
+      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
+      !!
+      !! ** Action :
+      !!       At each time step, mixed-layer averaged trends are stored in the 
+      !!       tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
+      !!       This array is known when trd_mld is called, at the end of the stp subroutine, 
+      !!       except for the purely vertical K_z diffusion term, which is embedded in the
+      !!       lateral diffusion trend.
+      !!
+      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
+      !!       from the lateral diffusion trend.
+      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
+      !!       arrays are updated.
+      !!       In III), called only once per analysis window, we compute the total trends,
+      !!       along with the residuals and the Asselin correction terms.
+      !!       In IV), the appropriate trends are written in the trends NetCDF file.
+      !!
+      !! References :
+      !!       - Vialard & al.
+      !!       - See NEMO documentation (in preparation)
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER ::   ji, jj, jk, jl, ik, it, itmod, jn
+      REAL(wp) ::   zavt, zfn, zfn2
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmltot             ! d(trc)/dt over the anlysis window (incl. Asselin)
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlres             ! residual = dh/dt entrainment term
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlatf             ! for storage only
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlrad             ! for storage only (for trb<0 corr in trcrad)
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmltot2            ! -+
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlres2            !  | working arrays to diagnose the trends
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmltrdm2           !  | associated with the time meaned ML
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlatf2            !  | passive tracers
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmlrad2            !  | (-> for trb<0 corr in trcrad)
+      !
+      CHARACTER (LEN=10) ::   clvar
+#if defined key_dimgout
+      INTEGER ::   iyear,imon,iday
+      CHARACTER(LEN=80) ::   cltext, clmode
+#endif
+      !!----------------------------------------------------------------------
+
+      ! Set-up pointers into sub-arrays of workspaces
+      CALL wrk_alloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad             )
+      CALL wrk_alloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
+
+      IF( nn_dttrc  /= 1  )   CALL ctl_stop( " Be careful, trends diags never validated " )
+
+      ! ======================================================================
+      ! I. Diagnose the purely vertical (K_z) diffusion trend
+      ! ======================================================================
+
+      ! ... These terms can be estimated by flux computation at the lower boundary of the ML 
+      !     (we compute (-1/h) * K_z * d_z( tracer ))
+
+      IF( ln_trcldf_iso ) THEN
+         !
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ik = nmld_trc(ji,jj)
+               zavt = fsavs(ji,jj,ik)
+               DO jn = 1, jptra
+                  IF( ln_trdtrc(jn) )    &
+                  tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
+                       &                    * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) )            &
+                       &                    / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
+               END DO
+            END DO
+         END DO
+
+         DO jn = 1, jptra
+            ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
+            IF( ln_trdtrc(jn) ) &
+                 tmltrd_trc(:,:,jpmxl_trc_ldf,jn) = tmltrd_trc(:,:,jpmxl_trc_ldf,jn) - tmltrd_trc(:,:,jpmxl_trc_zdf,jn)
+   
+         END DO
+         !      
+      ENDIF
+
+      IF ( cp_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration
+      ! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm. 
+      ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               DO jl = 1, jpltrd_trc
+                  CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
+               END DO
+            ENDIF
+         END DO
+      ENDIF
+      ! ======================================================================
+      ! II. Cumulate the trends over the analysis window
+      ! ======================================================================
+
+      ztmltrd2(:,:,:,:) = 0.e0   ;   ztmltot2(:,:,:)   = 0.e0     ! <<< reset arrays to zero
+      ztmlres2(:,:,:)   = 0.e0   ;   ztmlatf2(:,:,:)   = 0.e0
+      ztmlrad2(:,:,:)   = 0.e0
+
+      ! II.1 Set before values of vertically averages passive tracers
+      ! -------------------------------------------------------------
+      IF( kt > nittrc000 ) THEN
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               tmlb_trc   (:,:,jn) = tml_trc   (:,:,jn)
+               tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
+               tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
+            ENDIF
+         END DO
+      ENDIF
+
+      ! II.2 Vertically averaged passive tracers
+      ! ----------------------------------------
+      tml_trc(:,:,:) = 0.e0
+      DO jk = 1, jpktrd_trc ! - 1 ???
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) &
+               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
+         END DO
+      END DO
+
+      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window    
+      ! ------------------------------------------------------------------------
+      IF( kt == nittrc000 + nn_dttrc ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)    ???
+         !
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
+               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
+               
+               tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0   ;   tmltrd_atf_sumb_trc  (:,:,jn) = 0.e0
+               tmltrd_rad_sumb_trc  (:,:,jn) = 0.e0
+            ENDIF
+         END DO
+         
+         rmldbn_trc(:,:) = rmld_trc(:,:)
+         !
+      ENDIF
+
+      ! II.4 Cumulated trends over the analysis period
+      ! ----------------------------------------------
+      !
+      !         [  1rst analysis window ] [     2nd analysis window     ]                       
+      !
+      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
+      !                            ntrd                             2*ntrd       etc.
+      !     1      2     3     4    =5 e.g.                          =10
+      !
+      IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN                        ! ???
+         !
+         nmoymltrd = nmoymltrd + 1
+
+
+         ! ... Cumulate over BOTH physical contributions AND over time steps
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               DO jl = 1, jpltrd_trc
+                  tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
+               END DO
+            ENDIF
+         END DO
+
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               ! ... Special handling of the Asselin trend 
+               tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
+               tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
+
+               ! ... Trends associated with the time mean of the ML passive tracers
+               tmltrd_sum_trc    (:,:,:,jn) = tmltrd_sum_trc    (:,:,:,jn) + tmltrd_trc    (:,:,:,jn)
+               tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
+               tml_sum_trc       (:,:,jn)   = tml_sum_trc       (:,:,jn)   + tml_trc       (:,:,jn)
+            ENDIF
+         ENDDO
+
+         rmld_sum_trc      (:,:)     = rmld_sum_trc      (:,:)     + rmld_trc      (:,:)
+         !
+      ENDIF
+
+      ! ======================================================================
+      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
+      ! ======================================================================
+
+      ! Convert to appropriate physical units
+      tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
+
+      itmod = kt - nittrc000 + 1
+      it    = kt
+
+      MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN           ! nitend MUST be multiple of nn_trd_trc
+         !
+         ztmltot (:,:,:) = 0.e0                                   ! reset arrays to zero
+         ztmlres (:,:,:) = 0.e0
+         ztmltot2(:,:,:) = 0.e0
+         ztmlres2(:,:,:) = 0.e0
+      
+         zfn  = FLOAT( nmoymltrd )    ;    zfn2 = zfn * zfn
+         
+         ! III.1 Prepare fields for output ("instantaneous" diagnostics) 
+         ! -------------------------------------------------------------
+
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               !-- Compute total trends    (use rdttrc instead of rdt ???)
+               IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN  ! EULER-FORWARD schemes
+                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
+               ELSE                                                                     ! LEAP-FROG schemes
+                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
+               ENDIF
+               
+               !-- Compute residuals
+               ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
+                  &                                                 - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)  )
+               
+               !-- Diagnose Asselin trend over the analysis window 
+               ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
+               ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
+               
+         !-- Lateral boundary conditions
+               IF ( cp_cfg .NE. 'gyre' ) THEN
+                  CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
+                  CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
+               ENDIF
+
+
+#if defined key_diainstant
+               STOP 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.'
+#endif
+            ENDIF
+         END DO
+
+         ! III.2 Prepare fields for output ("mean" diagnostics) 
+         ! ----------------------------------------------------
+         
+         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
+         rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
+
+               !-- Compute passive tracer total trends
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
+               ztmltot2   (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) /  ( 2.*rdt )    ! now tracer unit is /sec
+            ENDIF
+         END DO
+
+         !-- Compute passive tracer residuals
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN
+               !
+               DO jl = 1, jpltrd_trc
+                  ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
+               END DO
+               
+               ztmltrdm2(:,:,jn) = 0.e0
+               DO jl = 1, jpltrd_trc
+                  ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
+               END DO
+               
+               ztmlres2(:,:,jn) =  ztmltot2(:,:,jn)  -       &
+                  & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
+                  &                     - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
+               !
+
+               !-- Diagnose Asselin trend over the analysis window
+               ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
+                  &                                               + tmltrd_atf_sumb_trc(:,:,jn)
+               ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
+                  &                                               + tmltrd_rad_sumb_trc(:,:,jn)
+
+         !-- Lateral boundary conditions 
+               IF ( cp_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration    
+                  CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
+                  CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
+                  DO jl = 1, jpltrd_trc
+                     CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
+                  END DO
+               ENDIF
+
+            ENDIF
+         END DO
+
+         ! * Debugging information *
+         IF( lldebug ) THEN
+            !
+            WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
+            WRITE(numout,*) '~~~~~~~~~~~  '
+            WRITE(numout,*)
+            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
+
+            DO jn = 1, jptra
+
+               IF( ln_trdtrc(jn) ) THEN
+                  WRITE(numout, *)
+                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
+                  
+                  WRITE(numout, *)
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
+                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
+                  
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
+                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
+                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
+                  
+                  WRITE(numout, *)
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
+                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn)) 
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
+                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
+                  
+                  WRITE(numout, *)
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
+                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
+                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
+                  
+                  WRITE(numout, *)
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
+                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
+                  
+                  WRITE(numout, *)
+                  DO jl = 1, jpltrd_trc
+                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
+                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
+                  END DO
+               
+                  WRITE(numout,*) 
+                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
+                  WRITE(numout,*)
+                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
+                  WRITE(numout,*)
+                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
+                  
+                  WRITE(numout,*) 
+                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
+                  WRITE(numout,*)
+                  
+                  WRITE(numout,*) '...................................................'
+                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
+                  DO jj = 5, 1, -1
+                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
+                  END DO
+                  
+                  WRITE(numout,*) 
+                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
+                  WRITE(numout,*)
+                  
+                  WRITE(numout,*) '...................................................'
+                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
+                  DO jj = 5, 1, -1
+                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
+                  END DO
+                  !
+               ENDIF
+               !
+            END DO
+
+
+97            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
+98            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
+99            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
+              WRITE(numout,*)
+            !
+         ENDIF
+
+         ! III.3 Time evolution array swap
+         ! -------------------------------
+         ! ML depth
+         rmldbn_trc(:,:)   = rmld_trc(:,:)
+         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
+         DO jn = 1, jptra
+            IF( ln_trdtrc(jn) ) THEN        
+               ! For passive tracer instantaneous diagnostics
+               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
+               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
+               
+               ! For passive tracer mean diagnostics
+               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
+               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
+               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
+               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
+               
+               
+               ! III.4 Convert to appropriate physical units
+               ! -------------------------------------------
+               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * rn_ucf_trc/zfn   ! instant diags
+               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * rn_ucf_trc/zfn
+               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * rn_ucf_trc/zfn
+               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * rn_ucf_trc/zfn
+               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
+               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * rn_ucf_trc/zfn2
+               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * rn_ucf_trc/zfn2
+               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * rn_ucf_trc/zfn2
+               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * rn_ucf_trc/zfn2
+               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * rn_ucf_trc/zfn2
+            ENDIF
+         END DO
+         !
+      ENDIF MODULO_NTRD
+
+      ! ======================================================================
+      ! IV. Write trends in the NetCDF file
+      ! ======================================================================
+
+      ! IV.1 Code for dimg mpp output
+      ! -----------------------------
+
+# if defined key_dimgout
+      STOP 'Not implemented'
+# else
+      
+      ! IV.2 Code for IOIPSL/NetCDF output
+      ! ----------------------------------
+
+      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
+         WRITE(numout,*) ' '
+         WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
+         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
+         WRITE(numout,*) ' '
+      ENDIF
+         
+      NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
+         !
+
+         DO jn = 1, jptra
+            !
+            IF( ln_trdtrc(jn) ) THEN
+               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
+               !-- Output the fields
+               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
+               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
+               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
+               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
+           
+               DO jl = 1, jpltrd_trc - 2
+                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),             &
+                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
+               END DO
+
+               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
+                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
+
+               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
+                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
+                     
+            ENDIF
+         END DO
+
+         IF( kt == nitend ) THEN 
+            DO jn = 1, jptra
+               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
+            END DO
+         ENDIF
+
+      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
+         
+         DO jn = 1, jptra
+            !
+            IF( ln_trdtrc(jn) ) THEN
+               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
+               !-- Output the fields
+               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
+
+               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
+               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
+               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
+
+               DO jl = 1, jpltrd_trc - 2
+                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),           &
+                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
+               END DO
+            
+               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
+                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
+
+               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
+                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
+
+            ENDIF 
+            !
+         END DO
+         IF( kt == nitend ) THEN 
+            DO jn = 1, jptra
+               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
+            END DO
+         ENDIF
+
+         !
+      ENDIF NETCDF_OUTPUT
+         
+      ! Compute the control surface (for next time step) : flag = on
+      icount = 1
+
+# endif /* key_dimgout */
+
+      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
+         !
+         ! Reset cumulative arrays to zero
+         ! -------------------------------         
+         nmoymltrd = 0
+         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
+         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
+         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
+         rmld_sum_trc       (:,:)     = 0.e0
+         !
+      ENDIF
+
+      ! ======================================================================
+      ! V. Write restart file
+      ! ======================================================================
+
+      IF( lrst_trc )   CALL trd_mxl_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
+
+      CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad             )
+      CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
+      !
+   END SUBROUTINE trd_mxl_trc
+
+
+   SUBROUTINE trd_mxl_bio( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mld  ***
+      !!
+      !! ** Purpose :  Compute and cumulate the mixed layer biological trends over an analysis
+      !!               period, and write NetCDF (or dimg) outputs.
+      !!
+      !! ** Method/usage :
+      !!          The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
+      !!          logical namelist variable) :
+      !!          1) to explain the difference between initial and final
+      !!             mixed-layer T & S (where initial and final relate to the
+      !!             current analysis window, defined by ntrd in the namelist)
+      !!          2) to explain the difference between the current and previous
+      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
+      !!             performed over each analysis window).
+      !!
+      !! ** Consistency check :
+      !!        If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
+      !!        entrainment) should be zero, at machine accuracy. Note that in the case
+      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
+      !!        over the first two analysis windows (except if restart).
+      !!        N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
+      !!             for checking residuals.
+      !!             On a NEC-SX5 computer, this typically leads to:
+      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
+      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
+      !!
+      !! ** Action :
+      !!       At each time step, mixed-layer averaged trends are stored in the
+      !!       tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
+      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
+      !!       except for the purely vertical K_z diffusion term, which is embedded in the
+      !!       lateral diffusion trend.
+      !!
+      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
+      !!       from the lateral diffusion trend.
+      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
+      !!       arrays are updated.
+      !!       In III), called only once per analysis window, we compute the total trends,
+      !!       along with the residuals and the Asselin correction terms.
+      !!       In IV), the appropriate trends are written in the trends NetCDF file.
+      !!
+      !! References :
+      !!       - Vialard & al.
+      !!       - See NEMO documentation (in preparation)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt                       ! ocean time-step index
+#if defined key_pisces_reduced
+      INTEGER  ::  jl, it, itmod
+      LOGICAL  :: llwarn  = .TRUE., lldebug = .TRUE.
+      REAL(wp) :: zfn, zfn2
+#if defined key_dimgout
+      INTEGER ::  iyear,imon,iday
+      CHARACTER(LEN=80) :: cltext, clmode
+#endif
+      !!----------------------------------------------------------------------
+      ! ... Warnings
+      IF( nn_dttrc  /= 1  ) CALL ctl_stop( " Be careful, trends diags never validated " )
+
+      ! ======================================================================
+      ! II. Cumulate the trends over the analysis window
+      ! ======================================================================
+
+      ztmltrdbio2(:,:,:) = 0.e0  ! <<< reset arrays to zero
+
+      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
+      ! ------------------------------------------------------------------------
+      IF( kt == nittrc000 + nn_dttrc ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
+         !
+         tmltrd_csum_ub_bio (:,:,:) = 0.e0
+         !
+      END IF
+
+      ! II.4 Cumulated trends over the analysis period
+      ! ----------------------------------------------
+      !
+      !         [  1rst analysis window ] [     2nd analysis window     ]
+      !
+      !
+      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
+      !                            ntrd                             2*ntrd       etc.
+      !     1      2     3     4    =5 e.g.                          =10
+      !
+      IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
+         !
+         nmoymltrdbio = nmoymltrdbio + 1
+
+         ! ... Trends associated with the time mean of the ML passive tracers
+         tmltrd_sum_bio    (:,:,:) = tmltrd_sum_bio    (:,:,:) + tmltrd_bio    (:,:,:)
+         tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
+         !
+      END IF
+
+      ! ======================================================================
+      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
+      ! ======================================================================
+
+      ! Convert to appropriate physical units
+      tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc
+
+      MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN      ! nitend MUST be multiple of ntrd
+         !
+         zfn  = float(nmoymltrdbio)    ;    zfn2 = zfn * zfn
+
+         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
+         ! -------------------------------------------------------------
+
+#if defined key_diainstant
+         STOP 'tmltrd_bio : key_diainstant was never checked within trdmxl. Comment this to proceed.'
+#endif
+         ! III.2 Prepare fields for output ("mean" diagnostics)
+         ! ----------------------------------------------------
+
+         ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
+
+         !-- Lateral boundary conditions
+         IF ( cp_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration 
+            ! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
+            DO jn = 1, jpdiabio
+              CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
+            ENDDO
+         ENDIF
+
+         IF( lldebug ) THEN
+            !
+            WRITE(numout,*) 'trd_mxl_bio : write trends in the Mixed Layer for debugging process:'
+            WRITE(numout,*) '~~~~~~~~~~~  '
+            WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
+            WRITE(numout,*)
+
+            DO jl = 1, jpdiabio
+              IF( ln_trdmxl_trc_instant ) THEN
+                  WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX  = ', jl, &
+                     & ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
+              ELSE
+                  WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX  = ', jl, &
+                     & ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
+              endif
+            END DO
+
+97          FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
+98          FORMAT(a10, i3, 2x, a30, 2x, g20.10)
+99          FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
+            WRITE(numout,*)
+            !
+         ENDIF
+
+         ! III.3 Time evolution array swap
+         ! -------------------------------
+
+         ! For passive tracer mean diagnostics
+         tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
+
+         ! III.4 Convert to appropriate physical units
+         ! -------------------------------------------
+         ztmltrdbio2    (:,:,:) = ztmltrdbio2    (:,:,:) * rn_ucf_trc/zfn2
+
+      END IF MODULO_NTRD
+
+      ! ======================================================================
+      ! IV. Write trends in the NetCDF file
+      ! ======================================================================
+
+      ! IV.1 Code for dimg mpp output
+      ! -----------------------------
+
+# if defined key_dimgout
+      STOP 'Not implemented'
+# else
+
+      ! IV.2 Code for IOIPSL/NetCDF output
+      ! ----------------------------------
+
+      ! define time axis
+      itmod = kt - nittrc000 + 1
+      it    = kt
+
+      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
+         WRITE(numout,*) ' '
+         WRITE(numout,*) 'trd_mxl_bio : write ML bio trends in the NetCDF file :'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '          ', TRIM(clhstnam), ' at kt = ', kt
+         WRITE(numout,*) '          N.B. nmoymltrdbio = ', nmoymltrdbio
+         WRITE(numout,*) ' '
+      END IF
+
+
+      ! 2. Start writing data
+      ! ---------------------
+
+      NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN    ! <<< write the trends for passive tracer instant. diags
+         !
+            DO jl = 1, jpdiabio
+               CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) ,            &
+                    &          it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
+            END DO
+
+
+         IF( kt == nitend )   CALL histclo( nidtrdbio )
+
+      ELSE    ! <<< write the trends for passive tracer mean diagnostics
+
+            DO jl = 1, jpdiabio
+               CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) ,            &
+                    &          it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
+            END DO
+
+            IF( kt == nitend )   CALL histclo( nidtrdbio )
+            !
+      END IF NETCDF_OUTPUT
+
+      ! Compute the control surface (for next time step) : flag = on
+      icountbio = 1
+
+
+# endif /* key_dimgout */
+
+      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
+         !
+         ! III.5 Reset cumulative arrays to zero
+         ! -------------------------------------
+         nmoymltrdbio = 0
+         tmltrd_csum_ln_bio (:,:,:) = 0.e0
+         tmltrd_sum_bio     (:,:,:) = 0.e0
+      END IF
+
+      ! ======================================================================
+      ! Write restart file
+      ! ======================================================================
+
+! restart write is done in trd_mxl_trc_write which is called by trd_mxl_bio (Marina)
+!
+#endif
+   END SUBROUTINE trd_mxl_bio
+
+
+   REAL FUNCTION sum2d( ztab )
+      !!----------------------------------------------------------------------
+      !! CD ??? prevoir d'utiliser plutot prtctl
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab      
+      !!----------------------------------------------------------------------
+      sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
+   END FUNCTION sum2d
+
+
+   SUBROUTINE trd_mxl_trc_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mxl_init  ***
+      !! 
+      !! ** Purpose :   computation of vertically integrated T and S budgets
+      !!      from ocean surface down to control surface (NetCDF output)
+      !!
+      !! ** Method/usage :
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: inum   ! logical unit
+      INTEGER :: ilseq, jl, jn, iiter
+      REAL(wp) ::   zjulian, zsto, zout
+      CHARACTER (LEN=40) ::   clop
+      CHARACTER (LEN=15) ::   csuff
+      CHARACTER (LEN=12) ::   clmxl
+      CHARACTER (LEN=16) ::   cltrcu
+      CHARACTER (LEN=10) ::   clvar
+
+      !!----------------------------------------------------------------------
+
+      ! ======================================================================
+      ! I. initialization
+      ! ======================================================================
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers                '
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
+         WRITE(numout,*)
+      ENDIF
+
+      
+      ! I.1 Check consistency of user defined preferences
+      ! -------------------------------------------------
+
+      IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
+         WRITE(numout,cform_err)
+         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
+         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
+         WRITE(numout,*) '                          you defined, nn_trd_trc   = ', nn_trd_trc
+         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
+         WRITE(numout,*) '                You should reconsider this choice.                        ' 
+         WRITE(numout,*) 
+         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
+         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
+         nstop = nstop + 1
+      ENDIF
+
+      ! * Debugging information *
+      IF( lldebug ) THEN
+         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
+         WRITE(numout,*) '               ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
+      ENDIF
+
+      IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
+         WRITE(numout,cform_err)
+         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
+         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
+         WRITE(numout,*) '                WHY? Everything in trdmxl_trc is coded for leap-frog, and '
+         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
+         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
+         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
+         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
+         WRITE(numout,*) 
+         nstop = nstop + 1
+      ENDIF
+
+      ! I.2 Initialize arrays to zero or read a restart file
+      ! ----------------------------------------------------
+      nmoymltrd   = 0
+
+      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
+      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
+      tmlradm_trc        (:,:,:)   = 0.e0
+
+      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
+      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
+
+#if defined key_pisces_reduced
+      nmoymltrdbio   = 0
+      tmltrd_sum_bio     (:,:,:) = 0.e0     ;   tmltrd_csum_ln_bio (:,:,:) = 0.e0
+      DO jl = 1, jp_pisces_trd
+          ctrd_bio(jl,1) = ctrbil(jl)   ! long name
+          ctrd_bio(jl,2) = ctrbio(jl)   ! short name
+       ENDDO
+#endif
+
+      IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
+         CALL trd_mxl_trc_rst_read
+      ELSE
+         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
+         tmlbn_trc          (:,:,:)   = 0.e0
+
+         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
+         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
+#if defined key_pisces_reduced
+         tmltrd_csum_ub_bio (:,:,:) = 0.e0
+#endif
+
+       ENDIF
+
+      icount = 1   ;   ionce  = 1  ! open specifier   
+
+#if defined key_pisces_reduced
+      icountbio = 1   ;   ioncebio  = 1  ! open specifier
+#endif
+
+      ! I.3 Read control surface from file ctlsurf_idx
+      ! ----------------------------------------------
+      IF( nn_ctls_trc == 1 ) THEN
+         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
+         READ ( inum ) nbol_trc
+         CLOSE( inum )
+      ENDIF
+
+      ! ======================================================================
+      ! II. netCDF output initialization
+      ! ======================================================================
+
+#if defined key_dimgout 
+      ???
+#else
+      ! clmxl = legend root for netCDF output
+      IF( nn_ctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
+         clmxl = 'Mixed Layer '
+      ELSE IF( nn_ctls_trc == 1 ) THEN                              ! control surface = read index from file 
+         clmxl = '      Bowl '
+      ELSE IF( nn_ctls_trc >= 2 ) THEN                              ! control surface = model level
+         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
+      ENDIF
+
+      ! II.1 Define frequency of output and means
+      ! -----------------------------------------
+      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
+      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cp time)
+      ENDIF
+#  if defined key_diainstant
+      IF( .NOT. ln_trdmxl_trc_instant ) THEN
+         STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
+      ENDIF
+      zsto = nn_trd_trc * rdt
+      clop = "inst("//TRIM(clop)//")"
+#  else
+      IF( ln_trdmxl_trc_instant ) THEN
+         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
+      ELSE
+         zsto = nn_trd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
+      ENDIF
+      clop = "ave("//TRIM(clop)//")"
+#  endif
+      zout = nn_trd_trc * rdt
+      iiter = ( nittrc000 - 1 ) / nn_dttrc
+
+      IF(lwp) WRITE (numout,*) '                netCDF initialization'
+
+      ! II.2 Compute julian date from starting date of the run
+      ! ------------------------------------------------------
+      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
+      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
+      IF(lwp) WRITE(numout,*)' '  
+      IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000               &
+           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
+           &   ,'Julian day : ', zjulian
+
+      ! II.3 Define the T grid trend file (nidtrd)
+      ! ------------------------------------------
+
+      !-- Define long and short names for the NetCDF output variables
+      !       ==> choose them according to trdmxl_trc_oce.F90 <==
+
+      ctrd_trc(jpmxl_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmxl_trc_xad    ,2) = "_xad"
+      ctrd_trc(jpmxl_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmxl_trc_yad    ,2) = "_yad"
+      ctrd_trc(jpmxl_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmxl_trc_zad    ,2) = "_zad"
+      ctrd_trc(jpmxl_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmxl_trc_ldf    ,2) = "_ldf"
+      ctrd_trc(jpmxl_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmxl_trc_zdf    ,2) = "_zdf"
+      ctrd_trc(jpmxl_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmxl_trc_bbl    ,2) = "_bbl"
+      ctrd_trc(jpmxl_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmxl_trc_dmp    ,2) = "_dmp"
+      ctrd_trc(jpmxl_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmxl_trc_sbc    ,2) = "_sbc"
+      ctrd_trc(jpmxl_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmxl_trc_sms    ,2) = "_sms"
+      ctrd_trc(jpmxl_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radb   ,2) = "_radb"
+      ctrd_trc(jpmxl_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radn   ,2) = "_radn"
+      ctrd_trc(jpmxl_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmxl_trc_atf    ,2) = "_atf"
+
+      DO jn = 1, jptra      
+      !-- Create a NetCDF file and enter the define mode 
+         IF( ln_trdtrc(jn) ) THEN
+            csuff="ML_"//ctrcnm(jn)
+            CALL dia_nam( clhstnam, nn_trd_trc, csuff )
+            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
+               &        1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
+      
+            !-- Define the ML depth variable
+            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
+               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+
+         ENDIF
+      END DO
+
+#if defined key_pisces_reduced
+          !-- Create a NetCDF file and enter the define mode
+          CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
+          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
+             &             1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )
+#endif
+
+      !-- Define physical units
+      IF( rn_ucf_trc == 1. ) THEN
+         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit 
+      ELSEIF ( rn_ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
+         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf 
+      ELSE
+         cltrcu = "unknown?"
+      ENDIF
+
+      !-- Define miscellaneous passive tracer mixed-layer variables 
+      IF( jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
+         STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb'  ! see below
+      ENDIF
+
+      DO jn = 1, jptra
+         !
+         IF( ln_trdtrc(jn) ) THEN
+            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
+            CALL histdef(nidtrd(jn), trim(clvar),           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
+              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
+            CALL histdef(nidtrd(jn), trim(clvar)//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
+              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
+            CALL histdef(nidtrd(jn), trim(clvar)//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
+              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
+         
+            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmxl_trc_atf
+               CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
+                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
+            END DO                                                                         ! if zsto=rdt above
+         
+            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), & 
+              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
+         
+            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1),   & 
+              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
+         !
+         ENDIF
+      END DO
+
+#if defined key_pisces_reduced
+      DO jl = 1, jp_pisces_trd
+         CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1))   ,            &
+             &    cltrcu, jpi, jpj, nh_tb, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
+      END DO                                                                         ! if zsto=rdt above
+#endif
+
+      !-- Leave IOIPSL/NetCDF define mode
+      DO jn = 1, jptra
+         IF( ln_trdtrc(jn) )  CALL histend( nidtrd(jn), snc4set )
+      END DO
+
+#if defined key_pisces_reduced
+      !-- Leave IOIPSL/NetCDF define mode
+      CALL histend( nidtrdbio, snc4set )
+
+      IF(lwp) WRITE(numout,*)
+       IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
+#endif
+
+#endif        /* key_dimgout */
+   END SUBROUTINE trd_mxl_trc_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                                       Empty module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trd_mxl_trc( kt )                                   ! Empty routine
+      INTEGER, INTENT( in) ::   kt
+      WRITE(*,*) 'trd_mxl_trc: You should not have seen this print! error?', kt
+   END SUBROUTINE trd_mxl_trc
+   SUBROUTINE trd_mxl_bio( kt )
+      INTEGER, INTENT( in) ::   kt
+      WRITE(*,*) 'trd_mxl_bio: You should not have seen this print! error?', kt
+   END SUBROUTINE trd_mxl_bio
+   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
+      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
+      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
+      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmxl            ! passive trc trend
+      WRITE(*,*) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
+   END SUBROUTINE trd_mxl_trc_zint
+   SUBROUTINE trd_mxl_trc_init                                    ! Empty routine
+      WRITE(*,*) 'trd_mxl_trc_init: You should not have seen this print! error?'
+   END SUBROUTINE trd_mxl_trc_init
+#endif
+
+   !!======================================================================
+END MODULE trdmxl_trc
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc_rst.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc_rst.F90	(revision 4990)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc_rst.F90	(revision 4990)
@@ -0,0 +1,220 @@
+MODULE trdmxl_trc_rst
+   !!======================================================================
+   !!                       ***  MODULE  trdmxl_rst  ***
+   !! Ocean dynamic :  Input/Output files for restart on mixed-layer diagnostics
+   !!======================================================================
+   !! History :  9.0  ! 07-03 (C. Deltel) Original code
+   !!----------------------------------------------------------------------
+  
+#if defined key_top && defined key_trdmxl_trc
+   !!----------------------------------------------------------------------
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
+   USE trc             ! for nn_dttrc ctrcnm
+   USE trdmxl_trc_oce  ! for lk_trdmxl_trc
+
+   IMPLICIT NONE
+   PRIVATE
+  
+   PUBLIC   trd_mxl_trc_rst_read    ! routine called by trd_mxl_init
+   PUBLIC   trd_mxl_trc_rst_write   ! routine called by step.F90
+  
+   INTEGER ::   nummldw_trc               ! logical unit for mld restart
+   !!---------------------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/TRD/trdmxl_rst.F90,v 1.6 2006/11/14 09:46:13 opalod Exp $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!---------------------------------------------------------------------------------
+  
+CONTAINS
+  
+
+    SUBROUTINE trd_mxl_trc_rst_write( kt )
+      !!--------------------------------------------------------------------------------
+      !!                  ***  SUBROUTINE trd_mxl_rst_wri  ***
+      !!                
+      !! ** Purpose :   Write mixed-layer diagnostics restart fields.
+      !!--------------------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
+      !
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character
+      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      CHARACTER (len=35) :: charout
+      INTEGER :: jl,  jk, jn               ! loop indice
+      !!--------------------------------------------------------------------------------
+
+      IF( kt == nitrst - nn_dttrc .OR. nitend - nit000 + 1 < 2 * nn_dttrc ) THEN ! idem trcrst.F90
+         IF( nitrst > 1.0e9 ) THEN
+            WRITE(clkt,*) nitrst
+         ELSE
+           WRITE(clkt,'(i8.8)') nitrst
+         ENDIF
+         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trdrst_trc_out)
+         IF(lwp) WRITE(numout,*) '             open ocean restart_mld_trc NetCDF  '//clname
+         CALL iom_open( clname, nummldw_trc, ldwrt = .TRUE., kiolib = jprstlib )
+      ENDIF
+
+      IF( kt == nitend .AND. lk_trdmxl_trc ) THEN
+
+         IF( kt == nitend .AND. lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'trdmxl_trc_rst: output for ML diags. restart, with trd_mxl_trc_rst_write routine'
+            WRITE(numout,*) '~~~~~~~~~~~~~~'
+            WRITE(numout,*)
+         ENDIF
+
+         IF( ln_trdmxl_trc_instant ) THEN 
+            !
+            DO jn = 1, jptra
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
+            END DO
+            !
+         ELSE
+            !
+            CALL iom_rstput( kt, nitrst, nummldw_trc, 'rmldbn_trc', rmldbn_trc )  ! 2D x 1
+            
+            !                                                          ! ===========
+            DO jn = 1, jptra                                           ! tracer loop
+               !                                                       ! ===========
+
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc (:,:,jn) ) ! 2D x jptra
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc   (:,:,jn) ) ! 2D x jptra
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc (:,:,jn) ) ! 2D x jptra
+
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! 2D x jptra
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tml_sumb_trc_'//ctrcnm(jn), tml_sumb_trc(:,:,jn) ) ! 2D x jptra
+               
+               DO jk = 1, jpltrd_trc
+                  IF( jk < 10 )   THEN
+                     WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I1)") ctrcnm(jn), jk
+                  ELSE
+                     WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I2)") ctrcnm(jn), jk
+                  ENDIF
+                  CALL iom_rstput( kt, nitrst, nummldw_trc, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
+               END DO
+               
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
+                    &           tmltrd_atf_sumb_trc(:,:,jn) ) ! 2D x jptra
+
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
+                    &           tmltrd_rad_sumb_trc(:,:,jn) ) ! 2D x jptra
+               !                                                       ! ===========
+            END DO                                                     ! tracer loop
+            !                                                          ! ===========
+#if defined key_pisces_reduced
+            DO jl = 1, jp_pisces_trd
+               CALL iom_rstput( kt, nitrst, nummldw_trc, 'tmltrd_csum_ub_bio'//ctrd_bio(jl,2), tmltrd_csum_ub_bio(:,:,jl) )
+            ENDDO
+#endif
+
+         ENDIF
+         
+         CALL iom_close( nummldw_trc )
+         lrst_trc = .TRUE.
+
+      ENDIF
+
+    END SUBROUTINE trd_mxl_trc_rst_write
+
+
+    SUBROUTINE trd_mxl_trc_rst_read
+      !!----------------------------------------------------------------------------
+      !!                   ***  SUBROUTINE trd_mxl_rst_lec  ***
+      !!                   
+      !! ** Purpose :   Read file for mixed-layer diagnostics restart.
+      !!----------------------------------------------------------------------------
+      INTEGER  ::  inum       ! temporary logical unit
+      !
+      CHARACTER (len=35) :: charout
+      INTEGER ::  jk, jn, jl     ! loop indice
+      INTEGER ::  jlibalt = jprstlib
+      LOGICAL ::  llok
+      !!-----------------------------------------------------------------------------
+      
+      IF(lwp)  THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' trd_mxl_trc_rst_read : read the NetCDF MLD restart file'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
+      ENDIF
+      
+      IF ( jprstlib == jprstdimg ) THEN
+        ! eventually read netcdf file (monobloc)  for restarting on different number of processors
+        ! if {cn_trdrst_trc_in}.nc exists, then set jlibalt to jpnf90
+        INQUIRE( FILE = TRIM(cn_trdrst_trc_in)//'.nc', EXIST = llok )
+        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
+      ENDIF
+
+      CALL iom_open( cn_trdrst_trc_in, inum, kiolib = jlibalt ) 
+      
+      IF( ln_trdmxl_trc_instant ) THEN 
+         
+         DO jn = 1, jptra
+            CALL iom_get( inum, jpdom_autoglo, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_autoglo, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_autoglo, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_autoglo, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
+         END DO
+         
+      ELSE
+         CALL iom_get( inum, jpdom_autoglo, 'rmldbn_trc', rmldbn_trc ) ! needed for rmld_sum
+         
+         !                                                          ! ===========
+         DO jn = 1, jptra                                           ! tracer loop
+            !                                                       ! ===========
+            CALL iom_get( inum, jpdom_autoglo, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_autoglo, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_autoglo, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc(:,:,jn) )
+
+            CALL iom_get( inum, jpdom_autoglo, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! needed for tml_sum
+            CALL iom_get( inum, jpdom_autoglo, 'tml_sumb_trc_'//ctrcnm(jn), tml_sumb_trc(:,:,jn) )
+            
+            DO jk = 1, jpltrd_trc
+               IF( jk < 10 )   THEN
+                  WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I1)") ctrcnm(jn), jk
+               ELSE
+                  WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I2)") ctrcnm(jn), jk
+               ENDIF
+               CALL iom_get( inum, jpdom_autoglo, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
+            END DO
+            
+            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
+                 &        tmltrd_atf_sumb_trc(:,:,jn) )
+
+            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
+                 &        tmltrd_rad_sumb_trc(:,:,jn) )
+            !                                                       ! ===========
+         END DO                                                     ! tracer loop
+         !                                                          ! ===========
+
+#if defined key_pisces_reduced
+         DO jl = 1, jp_pisces_trd
+            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_csum_ub_bio'//ctrd_bio(jl,2), tmltrd_csum_ub_bio(:,:,jl) )
+         ENDDO
+#endif
+         
+         CALL iom_close( inum )
+      ENDIF
+      
+    END SUBROUTINE trd_mxl_trc_rst_read
+  
+#else
+  !!=================================================================================
+  !!                       ***  MODULE  trdmxl_rst  ***
+  !! Ocean dynamic :  Input/Output files for restart on mixed-layer diagnostics
+  !!=================================================================================
+CONTAINS
+  SUBROUTINE trd_mxl_trc_rst_opn( kt )
+    WRITE(*,*) 'trd_mxl_trc_rst_opn: You should not have seen this print! error?', kt
+  END SUBROUTINE trd_mxl_trc_rst_opn
+  SUBROUTINE trd_mxl_trc_rst_write( kt )           !  No ML diags ==> empty routine
+    WRITE(*,*) 'trd_mxl_trc_rst_wri: You should not have seen this print! error?', kt
+  END SUBROUTINE trd_mxl_trc_rst_write
+  SUBROUTINE trd_mxl_trc_rst_read                  !  No ML Diags ==> empty routine
+  END SUBROUTINE trd_mxl_trc_rst_read
+#endif
+
+  !!=================================================================================
+END MODULE trdmxl_trc_rst
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90	(revision 4990)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90	(revision 4990)
@@ -0,0 +1,158 @@
+MODULE trdtrc
+   !!======================================================================
+   !!                       ***  MODULE  trdtrc  ***
+   !! Ocean diagnostics:  mixed layer passive tracer trends 
+   !!======================================================================
+   !! History :  3.0  !  2010-07  (C. Ethe)  Original code (from trdtrc.F90)
+   !!----------------------------------------------------------------------
+#if   defined key_top && ( defined key_trdmxl_trc   ||   defined key_trdtrc )
+   !!----------------------------------------------------------------------
+   !!   'key_trdmxl_trc'                  mixed layer trend diagnostics
+   !!   'key_trdtrc'                      3D trend diagnostics
+   !!----------------------------------------------------------------------
+   !!   trdtrc      : passive tracer trends 
+   !!----------------------------------------------------------------------
+   USE trc               ! tracer definitions (trn, trb, tra, etc.)
+   USE trcnam_trp
+   USE trd_oce
+   USE trdtrc_oce       ! definition of main arrays used for trends computations
+   USE trdmxl_trc        ! Mixed layer trends diag.
+   USE iom               ! I/O library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   INTERFACE trd_trc
+      MODULE PROCEDURE trd_trc_trp, trd_trc_bio
+   END INTERFACE
+
+   PUBLIC trd_trc
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Header:  $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trd_trc_trp( ptrtrd, kjn, ktrd, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_trc  ***
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )  ::   kt                                  ! time step
+      INTEGER, INTENT( in )  ::   kjn                                 ! tracer index
+      INTEGER, INTENT( in )  ::   ktrd                                ! tracer trend index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
+      CHARACTER (len=20) :: cltra
+      !!----------------------------------------------------------------------
+
+      IF( kt == nittrc000 ) THEN
+!         IF(lwp)WRITE(numout,*)
+!         IF(lwp)WRITE(numout,*) 'trd_trc:'
+!         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Mixed layer trends for passive tracers
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+#if defined key_trdmxl_trc  
+      IF( lk_trdmxl_trc .AND. ln_trdtrc( kjn ) ) THEN
+         !
+         SELECT CASE ( ktrd )
+         CASE ( jptra_xad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_xad, '3D', kjn )
+         CASE ( jptra_yad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_yad, '3D', kjn )
+         CASE ( jptra_zad     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_zad, '3D', kjn )
+         CASE ( jptra_ldf     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_ldf, '3D', kjn )
+         CASE ( jptra_bbl     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_bbl, '3D', kjn )
+         CASE ( jptra_zdf     )
+            IF( ln_trcldf_iso ) THEN
+               CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_ldf, '3D', kjn )
+            ELSE
+               CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_zdf, '3D', kjn )
+            ENDIF
+         CASE ( jptra_dmp     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_dmp , '3D', kjn )
+         CASE ( jptra_nsr     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_sbc , '2D', kjn )
+         CASE ( jptra_sms     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_sms , '3D', kjn )
+         CASE ( jptra_radb    )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_radb, '3D', kjn )
+         CASE ( jptra_radn    )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_radn, '3D', kjn )
+         CASE ( jptra_atf     )   ;   CALL trd_mxl_trc_zint( ptrtrd, jpmxl_trc_atf , '3D', kjn )
+         END SELECT
+         !
+      END IF
+#endif
+
+      IF( lk_trdtrc .AND. ln_trdtrc( kjn ) ) THEN
+         !
+         SELECT CASE( ktrd )
+         CASE( jptra_xad  )       ;    WRITE (cltra,'("XAD_",4a)')
+         CASE( jptra_yad  )       ;    WRITE (cltra,'("YAD_",4a)')
+         CASE( jptra_zad  )       ;    WRITE (cltra,'("ZAD_",4a)')
+         CASE( jptra_ldf  )       ;    WRITE (cltra,'("LDF_",4a)')
+         CASE( jptra_bbl  )       ;    WRITE (cltra,'("BBL_",4a)')
+         CASE( jptra_nsr  )       ;    WRITE (cltra,'("FOR_",4a)')
+         CASE( jptra_zdf  )       ;    WRITE (cltra,'("ZDF_",4a)')
+         CASE( jptra_dmp  )       ;    WRITE (cltra,'("DMP_",4a)')
+         CASE( jptra_sms  )       ;    WRITE (cltra,'("SMS_",4a)')
+         CASE( jptra_atf  )       ;    WRITE (cltra,'("ATF_",4a)')
+         CASE( jptra_radb )       ;    WRITE (cltra,'("RDB_",4a)')
+         CASE( jptra_radn )       ;    WRITE (cltra,'("RDN_",4a)')
+         END SELECT
+                                          cltra = TRIM(cltra)//TRIM(ctrcnm(kjn))
+                                          CALL iom_put( cltra,  ptrtrd(:,:,:) )
+         !
+      END IF
+
+   END SUBROUTINE trd_trc_trp
+
+   SUBROUTINE trd_trc_bio( ptrbio, ktrd, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_bio  ***
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT( in )  ::   kt                                  ! time step
+      INTEGER, INTENT( in )  ::   ktrd                                ! bio trend index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout )  ::   ptrbio  ! Bio trend
+      !!----------------------------------------------------------------------
+
+#if defined key_trdmxl_trc  
+      CALL trd_mxl_bio_zint( ptrbio, ktrd ) ! Verticaly integrated biological trends
+#endif
+
+   END SUBROUTINE trd_trc_bio
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                                       Empty module
+   !!----------------------------------------------------------------------
+
+   INTERFACE trd_trc
+      MODULE PROCEDURE trd_trc_trp, trd_trc_bio
+   END INTERFACE
+
+CONTAINS
+
+   SUBROUTINE trd_trc_trp( ptrtrd, kjn, ktrd, kt )
+      INTEGER               , INTENT( in )     ::   kt      ! time step
+      INTEGER               , INTENT( in )     ::   kjn     ! tracer index
+      INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
+      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
+      WRITE(*,*) 'trd_trc_trp : You should not have seen this print! error?', ptrtrd(1,1,1)
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
+   END SUBROUTINE trd_trc_trp
+
+   SUBROUTINE trd_trc_bio( ptrbio, ktrd, kt )
+      INTEGER               , INTENT( in )     ::   kt      ! time step
+      INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
+      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrbio  ! Temperature or U trend
+      WRITE(*,*) 'trd_trc_trp : You should not have seen this print! error?', ptrbio(1,1,1)
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
+      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
+   END SUBROUTINE trd_trc_bio
+
+#endif
+   !!======================================================================
+END MODULE trdtrc
Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90	(revision 4990)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90	(revision 4990)
@@ -0,0 +1,180 @@
+MODULE trdtrc_oce
+   !!======================================================================
+   !!                   ***  MODULE trdtrc_oce  ***
+   !! Ocean trends :   set tracer and momentum trend variables
+   !!======================================================================
+#if defined key_top   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   USE par_oce       ! ocean parameters
+   USE par_trc       ! passive tracers parameters
+
+   IMPLICIT NONE
+   PUBLIC
+
+   !                                         !!* Namelist namtoptrd:  diagnostics on passive tracers trends
+   INTEGER  ::    nn_trd_trc                  !: time step frequency dynamics and tracers trends
+   INTEGER  ::    nn_ctls_trc                 !: control surface type for trends vertical integration
+   REAL(wp) ::    rn_ucf_trc                  !: unit conversion factor (for netCDF trends outputs)
+   LOGICAL  ::    ln_trdmxl_trc_instant       !: flag to diagnose inst./mean ML trc trends
+   LOGICAL  ::    ln_trdmxl_trc_restart       !: flag to restart mixed-layer trc diagnostics
+   CHARACTER(len=50) ::  cn_trdrst_trc_in     !: suffix of pass. tracer restart name (input)
+   CHARACTER(len=50) ::  cn_trdrst_trc_out    !: suffix of pass. tracer restart name (output)
+   LOGICAL, DIMENSION(jptra) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist)
+
+# if defined key_trdtrc && defined key_iomput
+   LOGICAL, PARAMETER ::   lk_trdtrc = .TRUE. 
+# else
+   LOGICAL, PARAMETER ::   lk_trdtrc = .FALSE.   !: ML trend flag
+# endif
+
+# if defined key_trdmxl_trc   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_trdmxl_trc'                     mixed layer trends diagnostics
+   !!----------------------------------------------------------------------
+
+   LOGICAL, PARAMETER ::   lk_trdmxl_trc = .TRUE.    !: ML trend flag
+
+   INTEGER, PARAMETER ::            & !: mixed layer trends indices
+        jpmxl_trc_xad     =  1,     & !:     zonal       advection     
+        jpmxl_trc_yad     =  2,     & !:     meridonal   =========
+        jpmxl_trc_zad     =  3,     & !:     vertical    =========
+        jpmxl_trc_ldf     =  4,     & !:     lateral diffusion (geopot. or iso-neutral)
+        jpmxl_trc_zdf     =  5,     & !:     vertical diffusion (TKE)
+        jpmxl_trc_bbl     =  6,     & !:     bottom boundary layer (advective/diffusive)
+        jpmxl_trc_dmp     =  7,     & !:     internal restoring trend
+        jpmxl_trc_sbc     =  8,     & !:     forcing 
+        jpmxl_trc_sms     =  9,     & !:     sources minus sinks trend
+  !     jpmxl_trc_xxx     = xx,     & !:     add here any additional trend    (** AND UPDATE JPLTRD_TRC BELOW **)
+        jpmxl_trc_radn    = 10,     & !:     corr. trn<0 in trcrad
+        jpmxl_trc_radb    = 11,     & !:     corr. trb<0 in trcrad (like atf) (** MUST BE BEFORE THE LAST ONE **)
+        jpmxl_trc_atf     = 12        !:     asselin trend                    (** MUST BE    THE      LAST ONE**)
+
+   !! Trends diagnostics parameters
+   !!---------------------------------------------------------------------
+   INTEGER, PARAMETER :: jpltrd_trc = 12    !: number of mixed-layer trends arrays
+      
+   INTEGER            :: jpktrd_trc         !: max level for mixed-layer trends diag.
+
+   !! Arrays used for diagnosing mixed-layer trends 
+   !!---------------------------------------------------------------------
+   CHARACTER(LEN=80) :: clname_trc, ctrd_trc(jpltrd_trc+1,2)
+
+   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &
+      nmld_trc       , &                            !: mixed layer depth indexes 
+      nbol_trc                                   !: mixed-layer depth indexes when read from file
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  wkx_trc  !:
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  rmld_trc     !: ML depth (m) corresponding to nmld_trc
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  rmld_sum_trc !: needed to compute the leap-frog time mean of ML depth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  rmldbn_trc   !: idem
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &
+      tml_trc    ,                        &      !: \ "now" mixed layer temperature/salinity
+      tmlb_trc   ,                        &      !: /  and associated "before" fields
+      tmlbb_trc  ,                        &      !: \  idem, but valid at the 1rst time step of the
+      tmlbn_trc  ,                        &      !: /  current analysis window
+      tml_sum_trc,                        &      !: mixed layer T, summed over the current analysis period
+      tml_sumb_trc,                       &      !: idem, but from the previous analysis period
+      tmltrd_atf_sumb_trc,                &      !: Asselin trends, summed over the previous analysis period
+      tmltrd_rad_sumb_trc                        !: trends due to trb correction in trcrad.F90, summed over the
+                                                 !:     previous analysis period
+                                                 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &      
+      tmlatfb_trc, tmlatfn_trc ,          &      !: "before" Asselin contrib. at beginning of the averaging
+                                                 !:     period (i.e. last contrib. from previous such period)
+                                                 !:     and "now" Asselin contrib. to the ML trc. trends
+      tmlatfm_trc,                        &      !: accumulator for Asselin trends (needed for storage only)
+      tmlradb_trc, tmlradn_trc ,          &      !: similar to Asselin above, but for the trend due to trb
+                                                 !:     correction in trcrad.F90
+      tmlradm_trc                                !: accumulator for the previous trcrad trend
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  &
+      tmltrd_trc,                         &      !: \ physical contributions to the total trend (for T/S),
+                                                 !: / cumulated over the current analysis window
+      tmltrd_sum_trc,                     &      !: sum of these trends over the analysis period
+      tmltrd_csum_ln_trc,                 &      !: now cumulated sum of trends over the "lower triangle"
+      tmltrd_csum_ub_trc                         !: before (prev. analysis period) cumulated sum over the
+                                                 !: upper triangle
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &
+      tmltrdm_trc                                !: total cumulative trends over the analysis window
+
+# else
+   LOGICAL, PARAMETER ::   lk_trdmxl_trc = .FALSE.   !: ML trend flag
+# endif
+
+# if defined key_pisces_reduced
+   CHARACTER(LEN=80) :: clname_bio, ctrd_bio(jpdiabio,2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  &
+      tmltrd_bio,                         &      !: \ biological contributions to the total trend ,
+                                                 !: / cumulated over the current analysis window
+      tmltrd_sum_bio,                     &      !: sum of these trends over the analysis period
+      tmltrd_csum_ln_bio,                 &      !: now cumulated sum of trends over the "lower triangle"
+      tmltrd_csum_ub_bio                         !: before (prev. analysis period) cumulated sum over the
+                                                 !: upper triangle
+#endif
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/TRD/trdmxl_oce.F90,v 1.2 2005/03/27 18:35:23 opalod Exp $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION trd_trc_oce_alloc()
+      !!----------------------------------------------------------------------
+      !!         *** ROUTINE trd_trc_oce_alloc ***
+      !!----------------------------------------------------------------------
+      USE lib_mpp, ONLY: ctl_warn
+      INTEGER :: ierr(2)
+      !!----------------------------------------------------------------------
+      ierr(:) = 0
+      !
+# if defined key_trdmxl_trc
+      ALLOCATE(nmld_trc(jpi,jpj),          nbol_trc(jpi,jpj),           &
+               wkx_trc(jpi,jpj,jpk),       rmld_trc(jpi,jpj),           &
+               rmld_sum_trc(jpi,jpj),      rmldbn_trc(jpi,jpj),         &
+               tml_trc(jpi,jpj,jptra),     tmlb_trc(jpi,jpj,jptra),     &
+               tmlbb_trc(jpi,jpj,jptra),   tmlbn_trc(jpi,jpj,jptra),    &
+               tml_sum_trc(jpi,jpj,jptra), tml_sumb_trc(jpi,jpj,jptra), &
+               tmltrd_atf_sumb_trc(jpi,jpj,jptra),                      &
+               tmltrd_rad_sumb_trc(jpi,jpj,jptra),                      &
+               !
+               tmlatfb_trc(jpi,jpj,jptra), tmlatfn_trc(jpi,jpj,jptra),  &
+               tmlatfm_trc(jpi,jpj,jptra), tmlradb_trc(jpi,jpj,jptra),  &
+               tmlradn_trc(jpi,jpj,jptra), tmlradm_trc(jpi,jpj,jptra),  &
+               !
+               tmltrd_trc(jpi,jpj,jpltrd_trc,jptra)         , &
+               tmltrd_sum_trc(jpi,jpj,jpltrd_trc,jptra)     , &
+               tmltrd_csum_ln_trc(jpi,jpj,jpltrd_trc,jptra) , &
+               tmltrd_csum_ub_trc(jpi,jpj,jpltrd_trc,jptra) , &
+               !
+               tmltrdm_trc(jpi,jpj,jptra)                   , STAT=ierr(1) )
+#endif
+      !
+# if defined key_pisces_reduced
+      ALLOCATE( tmltrd_bio        (jpi,jpj,jpdiabio) ,     &
+         &      tmltrd_sum_bio    (jpi,jpj,jpdiabio) ,     &
+         &      tmltrd_csum_ln_bio(jpi,jpj,jpdiabio) ,     &
+         &      tmltrd_csum_ub_bio(jpi,jpj,jpdiabio) , STAT=ierr(2) )
+# endif
+      !
+      trd_trc_oce_alloc = MAXVAL(ierr)
+      !
+      IF( trd_trc_oce_alloc /= 0 )   CALL ctl_warn('trd_trc_oce_alloc: failed to allocate arrays')
+      !
+# if defined key_trdmxl_trc
+      jpktrd_trc = jpk      ! Initialise what used to be a parameter - max level for mixed-layer trends diag.
+# endif
+      !
+   END FUNCTION trd_trc_oce_alloc
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Empty module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE trdtrc_oce
