Index: /trunk/NEMO/TOP_SRC/TRP/trdmld_trc.F90
===================================================================
--- /trunk/NEMO/TOP_SRC/TRP/trdmld_trc.F90	(revision 1174)
+++ /trunk/NEMO/TOP_SRC/TRP/trdmld_trc.F90	(revision 1174)
@@ -0,0 +1,1172 @@
+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 trp_trc               ! tracer definitions (trn, trb, tra, etc.)
+   USE oce_trc           ! needed for namelist logicals, and euphotic layer arrays
+   USE trctrp_lec
+   USE trdmld_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 trdmld_trc_rst    ! restart for diagnosing the ML trends
+   USE prtctl            ! print control
+   USE sms_pisces        
+   USE sms_lobster
+   USE trcsms_cfc
+   USE trc
+   USE trcrst            ! for lrst_trc -> circ. dep. ??? we put lrst_trc in trc_oce
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trd_mod_trc                                             ! routine called by step.F90
+   PUBLIC trd_mld_trc
+   PUBLIC trd_mld_trc_init
+
+   CHARACTER (LEN=40) ::  clhstnam                                ! name of the trends NetCDF file
+   INTEGER ::   nmoymltrd
+   INTEGER ::   ndextrd1(jpi*jpj)
+   INTEGER, DIMENSION(jptra) ::   nidtrd, nh_t
+   INTEGER ::   ndimtrd1                        
+   INTEGER, SAVE ::  ionce, icount
+   LOGICAL :: llwarn  = .TRUE.                                    ! this should always be .TRUE.
+   LOGICAL :: lldebug = .TRUE.
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   TOP 1.0 , LOCEAN-IPSL (2007) 
+   !! $Header:  $ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trd_mod_trc( ptrtrd, kjn, ktrd, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mod_trc  ***
+      !!----------------------------------------------------------------------
+#if defined key_trcbbl_adv
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zun, zvn                   ! temporary arrays
+#else
+      USE oce_trc,   zun => un                                            ! When no bbl, zun == un
+      USE oce_trc,   zvn => vn                                            ! When no bbl, zvn == vn
+#endif
+      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
+      !!----------------------------------------------------------------------
+
+      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
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      SELECT CASE ( ktrd )
+         CASE ( jptrc_trd_xad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xad    , '3D', kjn )
+         CASE ( jptrc_trd_yad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yad    , '3D', kjn )
+         CASE ( jptrc_trd_zad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zad    , '3D', kjn )
+         CASE ( jptrc_trd_ldf     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf    , '3D', kjn )
+         CASE ( jptrc_trd_xei     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xei    , '3D', kjn )
+         CASE ( jptrc_trd_yei     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yei    , '3D', kjn )
+         CASE ( jptrc_trd_bbl     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbl    , '3D', kjn )
+         CASE ( jptrc_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 ( jptrc_trd_zei     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zei    , '3D', kjn )
+         CASE ( jptrc_trd_dmp     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_dmp    , '3D', kjn )
+         CASE ( jptrc_trd_sbc     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sbc    , '2D', kjn )
+#if defined key_lobster
+         CASE ( jptrc_trd_sms_sed )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_sed, '3D', kjn )
+         CASE ( jptrc_trd_sms_bio )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_bio, '3D', kjn )
+         CASE ( jptrc_trd_sms_exp )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_exp, '3D', kjn )
+#else
+         CASE ( jptrc_trd_sms     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms    , '3D', kjn )
+#endif
+         CASE ( jptrc_trd_bbc     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbc    , '3D', kjn )
+         CASE ( jptrc_trd_radb    )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radb   , '3D', kjn )
+         CASE ( jptrc_trd_radn    )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radn   , '3D', kjn )
+         CASE ( jptrc_trd_atf     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_atf    , '3D', kjn )
+      END SELECT
+
+
+   END SUBROUTINE trd_mod_trc
+
+   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), DIMENSION(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 ( nctls_trc )                                ! choice of the control surface
+            CASE ( -2  )   ;   STOP 'trdmld_trc : not ready '     !     -> isopycnal surface (see ???)
+#if defined key_pisces || defined key_lobster
+            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: )   ;   nctls_trc = MIN( nctls_trc, jpktrd_trc - 1 )
+                               nmld_trc(:,:) = nctls_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
+
+    END SUBROUTINE trd_mld_trc_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, 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, jn
+      REAL(wp) ::   zavt, zfn, zfn2
+      !!
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltot             ! d(trc)/dt over the anlysis window (incl. Asselin)
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlres             ! residual = dh/dt entrainment term
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlatf             ! for storage only
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlrad             ! for storage only (for trb<0 corr in trcrad)
+      !!
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltot2            ! -+
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlres2            !  | working arrays to diagnose the trends
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltrdm2           !  | associated with the time meaned ML
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlatf2            !  | passive tracers
+      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlrad2            !  | (-> for trb<0 corr in trcrad)
+      REAL(wp), DIMENSION(jpi,jpj,jpltrd_trc,jptra) ::  ztmltrd2  ! -+
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                       ! temporary array, used for eiv arrays
+      CHARACTER (LEN= 5) ::   clvar
+#if defined key_dimgout
+      INTEGER ::   iyear,imon,iday
+      CHARACTER(LEN=80) ::   cltext, clmode
+#endif
+      !!----------------------------------------------------------------------
+
+      IF( llwarn ) THEN                                           ! warnings
+         IF(      ( nittrc000 /= nit000   ) &
+              .OR.( ndttrc    /= 1        )    ) THEN
+
+            WRITE(numout,*) 'Be careful, trends diags never validated'
+            STOP 'Uncomment this line to proceed'
+         ENDIF
+      ENDIF
+
+      ! ======================================================================
+      ! 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 = avt(ji,jj,ik)
+               DO jn = 1, jptra
+                  IF( luttrd(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( luttrd(jn) ) &
+                 tmltrd_trc(:,:,jpmld_trc_ldf,jn) = tmltrd_trc(:,:,jpmld_trc_ldf,jn) - tmltrd_trc(:,:,jpmld_trc_zdf,jn)
+   
+         END DO
+         !      
+      ENDIF
+
+#if ! defined key_gyre
+      ! 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( luttrd(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( luttrd(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( luttrd(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 == 2 ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)    ???
+         !
+         DO jn = 1, jptra
+            IF( luttrd(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.( lrsttr ) ) THEN                        ! ???
+         !
+         nmoymltrd = nmoymltrd + 1
+
+
+         ! ... Cumulate over BOTH physical contributions AND over time steps
+         DO jn = 1, jptra
+            IF( luttrd(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( luttrd(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(:,:,:,:) * ucf_trc
+
+      MODULO_NTRD : IF( MOD( kt, ntrd_trc ) == 0 ) THEN           ! nitend MUST be multiple of ntrd_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( luttrd(jn) ) THEN
+               !-- Compute total trends    (use rdttrc instead of rdt ???)
+               IF ( ln_trcadv_smolar .OR. 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 ! defined key_gyre
+
+               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( luttrd(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( luttrd(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 ! defined key_gyre         
+               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( luttrd(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( luttrd(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)   * ucf_trc/zfn   ! instant diags
+               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * ucf_trc/zfn
+               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * ucf_trc/zfn
+               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * ucf_trc/zfn
+               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
+               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * ucf_trc/zfn2
+               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * ucf_trc/zfn2
+               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * ucf_trc/zfn2
+               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * ucf_trc/zfn2
+               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * 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( kt , ntrd_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
+         
+      it = kt - nit000 + 1
+
+      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
+         !
+
+         DO jn = 1, jptra
+            !
+            IF( luttrd(jn) ) THEN
+               !-- Specific treatment for EIV trends
+               !   WARNING : When eiv is switched on but key_diaeiv is not, we do NOT diagnose
+               !   u_eiv, v_eiv, and w_eiv : the exact eiv advective trends thus cannot be computed,
+               !   only their sum makes sense => mask directional contrib. to avoid confusion
+               z2d(:,:) = tmltrd_trc(:,:,jpmld_trc_xei,jn) + tmltrd_trc(:,:,jpmld_trc_yei,jn) &
+                    &   + tmltrd_trc(:,:,jpmld_trc_zei,jn)
+#if ( defined key_trcldf_eiv && defined key_diaeiv )
+               tmltrd_trc(:,:,jpmld_trc_xei,jn) = -999.
+               tmltrd_trc(:,:,jpmld_trc_yei,jn) = -999.
+               tmltrd_trc(:,:,jpmld_trc_zei,jn) = -999.
+#endif   
+               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), clvar         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
+               CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
+               CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
+           
+               DO jl = 1, jpltrd_trc - 2
+                  CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)),             &
+                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
+               END DO
+
+               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
+                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
+
+               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
+                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
+                     
+               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)),     &  ! now total EIV : jpltrd_trc + 1
+                    &          it, z2d(:,:), ndimtrd1, ndextrd1 )                     
+            !
+            ENDIF
+         END DO
+
+         IF( kt == nitend ) THEN 
+            DO jn = 1, jptra
+               IF( luttrd(jn) )  CALL histclo( nidtrd(jn) )
+            END DO
+         ENDIF
+
+      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
+         
+                 
+         DO jn = 1, jptra
+            !
+            IF( luttrd(jn) ) THEN
+               !-- Specific treatment for EIV trends
+               !   WARNING : see above
+               z2d(:,:) = ztmltrd2(:,:,jpmld_trc_xei,jn) + ztmltrd2(:,:,jpmld_trc_yei,jn) &
+                   &   + ztmltrd2(:,:,jpmld_trc_zei,jn)
+
+#if ( defined key_trcldf_eiv && defined key_diaeiv )
+               ztmltrd2(:,:,jpmld_trc_xei,jn) = -999.
+               ztmltrd2(:,:,jpmld_trc_yei,jn) = -999.
+               ztmltrd2(:,:,jpmld_trc_zei,jn) = -999.
+#endif
+               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), clvar         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
+               CALL histwrite( nidtrd(jn), clvar//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
+               CALL histwrite( nidtrd(jn), clvar//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
+
+               DO jl = 1, jpltrd_trc - 2
+                  CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)),           &
+                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
+               END DO
+            
+               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
+                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
+
+               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
+                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
+
+               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)),    &  ! now total EIV : jpltrd_trc + 1
+                 &          it, z2d(:,:), ndimtrd1, ndextrd1 )
+
+            ENDIF 
+            !
+         END DO
+         IF( kt == nitend ) THEN 
+            DO jn = 1, jptra
+               IF( luttrd(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( kt, ntrd_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)
+
+   END SUBROUTINE trd_mld_trc
+
+
+   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 :: ilseq, jl, jn
+      REAL(wp) ::   zjulian, zsto, zout
+      CHARACTER (LEN=21) ::   clold ='OLD'         ! open specifier (direct access files)
+      CHARACTER (LEN=21) ::   clunf ='UNFORMATTED' ! open specifier (direct access files)
+      CHARACTER (LEN=21) ::   clseq ='SEQUENTIAL'  ! open specifier (direct access files)
+      CHARACTER (LEN=40) ::   clop, cleiv
+      CHARACTER (LEN=15) ::   csuff
+      CHARACTER (LEN=12) ::   clmxl
+      CHARACTER (LEN=16) ::   cltrcu
+      CHARACTER (LEN= 5) ::   clvar
+      CHARACTER (LEN=80) ::   clname
+
+      NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, &
+                          ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd
+
+      !!----------------------------------------------------------------------
+
+      ! ======================================================================
+      ! 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 defined key_trcldf_eiv
+      IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN
+         WRITE(numout,cform_war)
+         WRITE(numout,*) '                You asked for ML diagnostics with iso-neutral diffusion   '
+         WRITE(numout,*) '                and eiv physics.                                          '
+         WRITE(numout,*) '                Yet, key_diaeiv is NOT switched on, so the eddy induced   '
+         WRITE(numout,*) '                velocity is not diagnosed.                                '
+         WRITE(numout,*) '                Therefore, we cannot deduce the eiv advective trends.     '
+         WRITE(numout,*) '                Only THE SUM of the i,j,k directional contributions then  '
+         WRITE(numout,*) '                makes sense => To avoid any confusion, we choosed to mask '
+         WRITE(numout,*) '                these i,j,k directional contributions (with -999.)        '
+         nwarn = nwarn + 1
+      ENDIF
+#  endif
+
+      IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_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, ntrd_trc   = ', ntrd_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
+
+      IF( ( lk_trdmld_trc ) .AND. ( n_cla == 1 ) ) THEN
+         WRITE(numout,cform_war)
+         WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
+         WRITE(numout,*) '                are not exact along the corresponding straits.            '
+         nwarn = nwarn + 1
+      ENDIF
+
+
+      ! * Debugging information *
+      IF( lldebug ) THEN
+         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
+         WRITE(numout,*) '               ln_trcadv_smolar = '     , ln_trcadv_smolar
+         WRITE(numout,*) '               ln_trdmld_trc_instant = ', ln_trdmld_trc_instant
+      ENDIF
+
+      IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN
+         WRITE(numout,cform_err)
+         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer Smolark. '
+         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,*) '                Smolarkiewicz scheme is Euler forward.                    '
+         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
+
+      IF( ln_trcadv_muscl .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
+
+      IF( ln_trcadv_muscl2 .AND. .NOT. ln_trdmld_trc_instant ) THEN
+         WRITE(numout,cform_err)
+         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL2    '
+         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( lrsttr .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 
+      ENDIF
+
+      ilseq  = 1   ;   icount = 1   ;   ionce  = 1  ! open specifier   
+
+      ! I.3 Read control surface from file ctlsurf_idx
+      ! ----------------------------------------------
+      IF( nctls_trc == 1 ) THEN
+         clname = 'ctlsurf_idx'
+         CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 )
+         REWIND( numbol )
+         READ  ( numbol ) nbol_trc
+      ENDIF
+
+      ! ======================================================================
+      ! II. netCDF output initialization
+      ! ======================================================================
+
+#if defined key_dimgout 
+      ???
+#else
+      ! clmxl = legend root for netCDF output
+      IF( nctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
+         clmxl = 'Mixed Layer '
+      ELSE IF( nctls_trc == 1 ) THEN                              ! control surface = read index from file 
+         clmxl = '      Bowl '
+      ELSE IF( nctls_trc >= 2 ) THEN                              ! control surface = model level
+         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls_trc
+      ENDIF
+
+      ! II.1 Define frequency of output and means
+      ! -----------------------------------------
+#  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 = ntrd_trc * rdt
+      clop ="inst(only(x))"
+#  else
+      IF( ln_trdmld_trc_instant ) THEN
+         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
+      ELSE
+         zsto = ntrd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
+      ENDIF
+      clop ="ave(only(x))"
+#  endif
+      zout = ntrd_trc * rdt
+
+      IF(lwp) WRITE (numout,*) '                netCDF initialization'
+
+      ! II.2 Compute julian date from starting date of the run
+      ! ------------------------------------------------------
+      CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
+      IF(lwp) WRITE(numout,*)' '  
+      IF(lwp) WRITE(numout,*)' Date 0 used :', nit000                  &
+           &   ,' 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 <==
+
+#if defined key_diaeiv
+      cleiv = " (*** only total EIV is meaningful ***)"           ! eiv advec. trends require u_eiv, v_eiv
+#else
+      cleiv = " "
+#endif
+      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_xei    ,1) = " Zonal EIV advection"//cleiv      ;   ctrd_trc(jpmld_trc_xei    ,2) = "_xei"
+      ctrd_trc(jpmld_trc_yei    ,1) = " Merid. EIV advection"//cleiv     ;   ctrd_trc(jpmld_trc_yei    ,2) = "_yei"
+      ctrd_trc(jpmld_trc_zei    ,1) = " Vertical EIV advection"//cleiv   ;   ctrd_trc(jpmld_trc_zei    ,2) = "_zei"
+      ctrd_trc(jpmld_trc_bbc    ,1) = " Geothermal flux"                 ;   ctrd_trc(jpmld_trc_bbc    ,2) = "_bbc"
+      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"
+#if defined key_lobster
+      ctrd_trc(jpmld_trc_sms_sed,1) = " Sources minus sinks : sed"       ;   ctrd_trc(jpmld_trc_sms_sed,2) = "_sms_sed"
+      ctrd_trc(jpmld_trc_sms_bio,1) = " Sources minus sinks : bio"       ;   ctrd_trc(jpmld_trc_sms_bio,2) = "_sms_bio"
+      ctrd_trc(jpmld_trc_sms_exp,1) = " Sources minus sinks : exp"       ;   ctrd_trc(jpmld_trc_sms_exp,2) = "_sms_exp"
+#else
+      ctrd_trc(jpmld_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmld_trc_sms    ,2) = "_sms"
+#endif                                                                                                 
+      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"
+      ctrd_trc(jpltrd_trc+1     ,1) = " Total EIV"//cleiv                ;   ctrd_trc(jpltrd_trc+1     ,2) = "_tei"
+
+      DO jn = 1, jptra      
+      !-- Create a NetCDF file and enter the define mode 
+         IF( luttrd(jn) ) THEN
+            csuff="TD_"//ctrcnm(jn)
+            CALL dia_nam( clhstnam, ntrd_trc, csuff )
+            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
+               &        1, jpi, 1, jpj, 0, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )
+      
+            !-- 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
+
+      !-- Define physical units
+      IF( ucf_trc == 1. ) THEN
+         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit 
+      ELSEIF ( 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
+#if defined key_lobster
+      IF( lldebug ) THEN
+         DO jn = 1, jptra
+            WRITE(numout, *) 'TRC jpdet=', jpdet, ' jpnh4=', jpnh4
+            WRITE(numout, *) 'TRC short title  ctrcnm  jn=", jn, " : ', ctrcnm(jn)
+            WRITE(numout, *) 'TRC trim(ctrcnm(jn))//"_tot" = ', trim(ctrcnm(jn))//"ml_tot"  ! tml_tot -> detml_tot
+         END DO
+         CALL flush(numout)
+      ENDIF
+#else
+!!      Error : this is not ready (PISCES)
+#endif      
+
+      DO jn = 1, jptra
+         !
+         IF( luttrd(jn) ) THEN
+            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
+            CALL histdef(nidtrd(jn), 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), 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), 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//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//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//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
+         
+            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpltrd_trc+1,2)),  clmxl//" "//clvar//ctrd_trc(jpltrd_trc+1 ,1),   & 
+              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! Total EIV 
+         !
+         ENDIF
+      END DO
+
+      !-- Leave IOIPSL/NetCDF define mode
+      DO jn = 1, jptra
+         IF( luttrd(jn) )  CALL histend( nidtrd(jn) )
+      END DO
+
+#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_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
+   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/NEMO/TOP_SRC/TRP/trdmld_trc_oce.F90
===================================================================
--- /trunk/NEMO/TOP_SRC/TRP/trdmld_trc_oce.F90	(revision 1174)
+++ /trunk/NEMO/TOP_SRC/TRP/trdmld_trc_oce.F90	(revision 1174)
@@ -0,0 +1,161 @@
+MODULE trdmld_trc_oce
+   !!======================================================================
+   !!                   ***  MODULE trdmld_trc_oce  ***
+   !! Ocean trends :   set tracer and momentum trend variables
+   !!======================================================================
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $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 (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   '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  ::    ntrd_trc                !: time step frequency dynamics and tracers trends
+   INTEGER  ::    nctls_trc               !: control surface type for trends vertical integration
+   REAL(wp) ::    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
+   LOGICAL, DIMENSION (jptra) ::   luttrd   !: large trends diagnostic to write or not (namelist)
+
+   !!* Passive tracers trends indexes
+   INTEGER, PARAMETER ::   jptrc_trd_xad     =  1   !: x- horizontal advection
+   INTEGER, PARAMETER ::   jptrc_trd_yad     =  2   !: y- horizontal advection
+   INTEGER, PARAMETER ::   jptrc_trd_zad     =  3   !: z- vertical   advection
+   INTEGER, PARAMETER ::   jptrc_trd_ldf     =  4   !: lateral       diffusion
+   INTEGER, PARAMETER ::   jptrc_trd_zdf     =  5   !: vertical diffusion (Kz)
+   INTEGER, PARAMETER ::   jptrc_trd_xei     =  6   !: x- horiz. EIV advection
+   INTEGER, PARAMETER ::   jptrc_trd_yei     =  7   !: y- horiz. EIV advection
+   INTEGER, PARAMETER ::   jptrc_trd_zei     =  8   !: z- vert.  EIV advection
+   INTEGER, PARAMETER ::   jptrc_trd_bbc     =  9   !: Bottom Boundary Condition (geoth. flux)
+   INTEGER, PARAMETER ::   jptrc_trd_bbl     = 10   !: Bottom Boundary Layer (diffusive/convective)
+   INTEGER, PARAMETER ::   jptrc_trd_dmp     = 11   !: damping
+   INTEGER, PARAMETER ::   jptrc_trd_sbc     = 12   !: surface boundary condition
+#if defined key_lobster
+   INTEGER, PARAMETER ::   jptrc_trd_sms_sed = 13   !: sources m. sinks : sedimentation
+   INTEGER, PARAMETER ::   jptrc_trd_sms_bio = 14   !: sources m. sinks : bio
+   INTEGER, PARAMETER ::   jptrc_trd_sms_exp = 15   !: sources m. sinks : bio
+   INTEGER, PARAMETER ::   jptrc_trd_radn    = 16   !: corr. trn<0 in trcrad
+   INTEGER, PARAMETER ::   jptrc_trd_radb    = 17   !: corr. trb<0 in trcrad (like atf)
+   INTEGER, PARAMETER ::   jptrc_trd_atf     = 18   !: Asselin correction
+#else
+   INTEGER, PARAMETER ::   jptrc_trd_sms     = 13   !: sources m. sinks
+   INTEGER, PARAMETER ::   jptrc_trd_radn    = 14   !: corr. trn<0 in trcrad
+   INTEGER, PARAMETER ::   jptrc_trd_radb    = 15   !: corr. trb<0 in trcrad (like atf)
+   INTEGER, PARAMETER ::   jptrc_trd_atf     = 16   !: Asselin correction
+#endif
+
+#if defined key_trdmld_trc
+   !!----------------------------------------------------------------------
+   !!   'key_trdmld_trc'                     mixed layer trends diagnostics
+   !!----------------------------------------------------------------------
+
+   LOGICAL, PARAMETER ::   lk_trdmld_trc = .TRUE.    !: ML trend flag
+   LOGICAL, PARAMETER ::   l_trdtrc = .TRUE.         !: passive tracers trend flag
+
+   INTEGER, PARAMETER ::            & !: mixed layer trends indices
+        jpmld_trc_xad     =  1,     & !:     zonal      \
+        jpmld_trc_yad     =  2,     & !:     meridonal   > advection
+        jpmld_trc_zad     =  3,     & !:     vertical   /
+        jpmld_trc_ldf     =  4,     & !:     lateral diffusion (geopot. or iso-neutral)
+        jpmld_trc_zdf     =  5,     & !:     vertical diffusion (TKE)
+        jpmld_trc_xei     =  6,     & !:     zonal      \
+        jpmld_trc_yei     =  7,     & !:     meridional  > eddy induced advection
+        jpmld_trc_zei     =  8,     & !:     vertical   / 
+        jpmld_trc_bbc     =  9,     & !:     geothermal flux
+        jpmld_trc_bbl     = 10,     & !:     bottom boundary layer (advective/diffusive)
+        jpmld_trc_dmp     = 11,     & !:     internal restoring trend
+        jpmld_trc_sbc     = 12,     & !:     forcing 
+#if  defined key_lobster
+        jpmld_trc_sms_sed = 13,     & !:     sources minus sinks trend
+        jpmld_trc_sms_bio = 14,     & !:     sources minus sinks trend
+        jpmld_trc_sms_exp = 15,     & !:     sources minus sinks trend
+  !     jpmld_trc_xxx     = xx,     & !:     add here any additional trend    (** AND UPDATE JPLTRD_TRC BELOW **)
+        jpmld_trc_radn    = 16,     & !:     corr. trn<0 in trcrad
+        jpmld_trc_radb    = 17,     & !:     corr. trn<0 in trcrad (like atf) (** MUST BE BEFORE THE LAST ONE **)
+        jpmld_trc_atf     = 18        !:     asselin trend                    (** MUST BE    THE     LAST ONE **)
+#else                                                      
+        jpmld_trc_sms     = 13,     & !:     sources minus sinks trend
+  !     jpmld_trc_xxx     = xx,     & !:     add here any additional trend    (** AND UPDATE JPLTRD_TRC BELOW **)
+        jpmld_trc_radn    = 14,     & !:     corr. trn<0 in trcrad
+        jpmld_trc_radb    = 15,     & !:     corr. trb<0 in trcrad (like atf) (** MUST BE BEFORE THE LAST ONE **)
+        jpmld_trc_atf     = 16        !:     asselin trend                    (** MUST BE    THE      LAST ONE**)
+#endif
+
+   !! Trends diagnostics parameters
+   !!---------------------------------------------------------------------
+   INTEGER, PARAMETER ::            &
+#if defined key_lobster
+      jpltrd_trc = 18,                  & !: number of mixed-layer trends arrays
+#else
+      jpltrd_trc = 16,                  & !: number of mixed-layer trends arrays
+#endif
+      jpktrd_trc = jpk                    !: 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, DIMENSION(jpi,jpj) ::   &
+      nmld_trc       , &                            !: mixed layer depth indexes 
+      nbol_trc                                   !: mixed-layer depth indexes when read from file
+
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::  wkx_trc  !:
+
+   REAL(wp), DIMENSION(jpi,jpj) ::  rmld_trc     !: ML depth (m) corresponding to nmld_trc
+   REAL(wp), DIMENSION(jpi,jpj) ::  rmld_sum_trc !: needed to compute the leap-frog time mean of ML depth
+   REAL(wp), DIMENSION(jpi,jpj) ::  rmldbn_trc   !: idem
+
+   REAL(wp), DIMENSION(jpi,jpj,jptra) ::  &
+      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), DIMENSION(jpi,jpj,jptra) ::  &      
+      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), DIMENSION(jpi,jpj,jpltrd_trc,jptra) ::  &
+      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), DIMENSION(jpi,jpj,jptra) ::  &
+      tmltrdm_trc                                !: total cumulative trends over the analysis window
+
+#else
+   LOGICAL, PARAMETER ::   lk_trdmld_trc = .FALSE.   !: ML trend flag
+   LOGICAL, PARAMETER ::   l_trdtrc = .FALSE.             !: passive tracers trend flag
+#endif
+
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Empty module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+#endif
+
+
+END MODULE trdmld_trc_oce
Index: /trunk/NEMO/TOP_SRC/TRP/trdmld_trc_rst.F90
===================================================================
--- /trunk/NEMO/TOP_SRC/TRP/trdmld_trc_rst.F90	(revision 1174)
+++ /trunk/NEMO/TOP_SRC/TRP/trdmld_trc_rst.F90	(revision 1174)
@@ -0,0 +1,214 @@
+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 oce_trc
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
+   USE trcrst          ! only for lrst_trc
+   USE trc             ! for ndttrc ctrcnm
+   USE trdmld_trc_oce  ! for lk_trdmld_trc
+
+   IMPLICIT NONE
+   PRIVATE
+  
+   PUBLIC   trd_mld_trc_rst_opn
+   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
+   !!---------------------------------------------------------------------------------
+   !! OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! $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 (modipsl/doc/NEMO_CeCILL.txt)
+   !!---------------------------------------------------------------------------------
+  
+CONTAINS
+  
+    SUBROUTINE trd_mld_trc_rst_opn( kt )
+      !!--------------------------------------------------------------------------------
+      !!                  ***  SUBROUTINE trd_mld_trc_rst_opn  ***
+      !!                
+      !! ** Purpose :   open 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
+      !!--------------------------------------------------------------------------------
+      
+      IF ( lk_trdmld_trc ) THEN
+         IF( kt == nitrst - ndttrc .OR. nitend - nit000 + 1 < 2 * ndttrc ) 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))//"_restart_mld_trc"
+            IF(lwp) WRITE(numout,*) '             open ocean restart_mld_trc NetCDF file: '//clname
+            CALL iom_open( clname, nummldw_trc, ldwrt = .TRUE., kiolib = jprstlib )
+         ENDIF
+      ENDIF
+
+    END SUBROUTINE trd_mld_trc_rst_opn
+
+    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=35) :: charout
+      INTEGER ::   jk, jn               ! loop indice
+      !!--------------------------------------------------------------------------------
+
+      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
+            !                                                          ! ===========
+         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     ! loop indice
+      !!-----------------------------------------------------------------------------
+      
+      IF(lwp)  THEN
+         WRITE(numout,*)
+         WRITE(numout,*) ' trd_mld_trc_rst_read : read the NetCDF MLD restart file'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
+      ENDIF
+      
+      CALL iom_open( 'restart_mld_trc', inum, kiolib = jprstlib ) 
+      
+      IF( ln_trdmld_trc_instant ) THEN 
+         
+         DO jn = 1, jptra
+            CALL iom_get( inum, jpdom_local, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_local, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_local, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_local, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
+         END DO
+         
+      ELSE
+         CALL iom_get( inum, jpdom_local, 'rmldbn_trc', rmldbn_trc ) ! needed for rmld_sum
+         
+         !                                                          ! ===========
+         DO jn = 1, jptra                                           ! tracer loop
+            !                                                       ! ===========
+            CALL iom_get( inum, jpdom_local, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_local, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_local, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc(:,:,jn) )
+
+            CALL iom_get( inum, jpdom_local, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! needed for tml_sum
+            CALL iom_get( inum, jpdom_local, '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_local, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
+            END DO
+            
+            CALL iom_get( inum, jpdom_local, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
+                 &        tmltrd_atf_sumb_trc(:,:,jn) )
+
+            CALL iom_get( inum, jpdom_local, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
+                 &        tmltrd_rad_sumb_trc(:,:,jn) )
+            !                                                       ! ===========
+         END DO                                                     ! tracer loop
+         !                                                          ! ===========
+         
+         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
