Index: trunk/NEMO/OPA_SRC/TRD/trdicp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 521)
@@ -702,11 +702,15 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )       ! Empty routine
+   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )       ! Empty routine
       REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
-      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype
+      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas     ! number of passage
+      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
+          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype, clpas
    END SUBROUTINE trd_2d
-   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype )       ! Empty routine
+   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype, clpas )       ! Empty routine
       REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
-      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype
+      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas     ! number of passage
+      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
+          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype, clpas
    END SUBROUTINE trd_3d
    SUBROUTINE trd_icp_init               ! Empty routine
Index: trunk/NEMO/OPA_SRC/TRD/trdmld.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 521)
@@ -35,5 +35,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE diadimg         ! dimg direct access file format output
-   USE trdmld_rst , ONLY : trd_mld_rst_read  ! restart for diagnosing the ML trends
+   USE trdmld_rst      ! restart for diagnosing the ML trends
    USE prtctl          ! Print control
 
@@ -230,5 +230,5 @@
       !!
       INTEGER :: ji, jj, jk, jl, ik, it
-      LOGICAL :: lldebug = .TRUE.
+      LOGICAL :: lldebug = .FALSE.
       REAL(wp) :: zavt, zfn, zfn2
       REAL(wp) ,DIMENSION(jpi,jpj) ::     &
@@ -253,20 +253,20 @@
       ! ... These terms can be estimated by flux computation at the lower boundary of the ML 
       !     (we compute (-1/h) * K_z * d_z( T ) and (-1/h) * K_z * d_z( S ))
-      DO jj = 1,jpj
-         DO ji = 1,jpi
-            ik = nmld(ji,jj)
-            zavt = avt(ji,jj,ik)
-            tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
-               &                      * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) )         &
-               &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
-            zavt = fsavs(ji,jj,ik)
-            smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
-               &                      * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) )         &
-               &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
-         END DO
-      END DO
-
-      ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
       IF( ln_traldf_iso ) THEN
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ik = nmld(ji,jj)
+               zavt = avt(ji,jj,ik)
+               tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
+                  &                      * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) )         &
+                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
+               zavt = fsavs(ji,jj,ik)
+               smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
+                  &                      * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) )         &
+                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
+            END DO
+         END DO
+
+         ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
          tmltrd(:,:,jpmld_ldf) = tmltrd(:,:,jpmld_ldf) - tmltrd(:,:,jpmld_zdf)
          smltrd(:,:,jpmld_ldf) = smltrd(:,:,jpmld_ldf) - smltrd(:,:,jpmld_zdf)
@@ -698,4 +698,10 @@
       END IF
 
+      ! ======================================================================
+      ! V. Write restart file
+      ! ======================================================================
+
+      CALL trd_mld_rst_write( kt ) 
+
    END SUBROUTINE trd_mld
 
Index: trunk/NEMO/OPA_SRC/TRD/trdmld_rst.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmld_rst.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/TRD/trdmld_rst.F90	(revision 521)
@@ -12,5 +12,6 @@
    USE in_out_manager  ! I/O manager
    USE daymod          ! calendar
-   USE ioipsl          !
+   USE iom             ! I/O module
+   USE restart         ! ocean restart 
 
    IMPLICIT NONE
@@ -32,78 +33,61 @@
       !!                
       !! ** Purpose :   Write mixed-layer diagnostics restart fields.
-      !!-----------------------------------------------------------------------------------
+      !!--------------------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!    
-      LOGICAL ::   llbon
-      CHARACTER (len=50) ::   clname, cln
-      REAL(wp) ::   zdate0
-      REAL(wp), DIMENSION(1) :: zdept
-      INTEGER ::   ic, jc, itime, inumwrs_mld
-      !!-----------------------------------------------------------------------------
+      !
+      CHARACTER (len=35) :: charout
+      INTEGER ::   jk                 ! loop indice
+      !!--------------------------------------------------------------------------------
     
       IF( ( mod( kt, nstock ) == 0 ) .OR. ( kt == nitend ) ) THEN
              
-         inumwrs_mld = 47
-         !-- Delete the restart file if it exists 
-         INQUIRE( FILE=crestart, EXIST=llbon )
-         IF(llbon) THEN
-            OPEN( UNIT=inumwrs_mld, FILE=crestart, STATUS='old' )
-            CLOSE( inumwrs_mld, STATUS='delete' )
-         ENDIF
-
-         !-- Name of the new restart file
-         ic = 1
-         DO jc = 1, 16
-            IF( cexper(jc:jc) /= ' ' )   ic = jc
-         END DO
-         WRITE(cln,'("_",i4.4,i2.2,i2.2,"_restart_mld")') nyear, nmonth, nday
-         clname = cexper(1:ic)//cln
-         ic = 1
-         DO jc = 1, 48
-            IF( clname(jc:jc) /= ' ' ) ic = jc
-         END DO
-         crestart = clname(1:ic)//".nc"
-         itime = 0
-         CALL ymds2ju( nyear, nmonth, nday, 0.e0, zdate0 )
-
          IF(lwp) THEN
             WRITE(numout,*)
             WRITE(numout,*) 'trdmld_rst: output for ML diags. restart, with trd_mld_rst_write routine'
             WRITE(numout,*) '~~~~~~~~~~'
-            WRITE(numout,*) '            in file : ', TRIM(crestart), ' at it= ', kt, ' date= ', ndastp
             WRITE(numout,*)
          ENDIF
 
-         !-- Create the NetCDF restart file and write tje appropriate fields
-         !   N.B. In this section, 3rd dimension in arrays is NOT depth 
-         zdept(1) = 1.
-         CALL restini( 'NONE', jpi, jpj, glamt, gphit, 1, zdept, clname,   &
-            &          itime, zdate0, rdt*nstock ,inumwrs_mld, domain_id=nidom )
-       
          IF( ln_trdmld_instant ) THEN 
-            CALL restput( inumwrs_mld, 'tmlbb'           , jpi, jpj,      1, 0, tmlbb           )
-            CALL restput( inumwrs_mld, 'tmlbn'           , jpi, jpj,      1, 0, tmlbn           )
-            CALL restput( inumwrs_mld, 'tmlatfb'         , jpi, jpj,      1, 0, tmlatfb         )
+            !-- Temperature
+            CALL iom_rstput( kt, nitrst, nummldw, 'tmlbb'           , tmlbb           )
+            CALL iom_rstput( kt, nitrst, nummldw, 'tmlbn'           , tmlbn           )
+            CALL iom_rstput( kt, nitrst, nummldw, 'tmlatfb'         , tmlatfb         )
 
-            CALL restput( inumwrs_mld, 'smlbb'           , jpi, jpj,      1, 0, smlbb           )
-            CALL restput( inumwrs_mld, 'smlbn'           , jpi, jpj,      1, 0, smlbn           )
-            CALL restput( inumwrs_mld, 'smlatfb'         , jpi, jpj,      1, 0, smlatfb         )
+            !-- Salinity
+            CALL iom_rstput( kt, nitrst, nummldw, 'smlbb'           , smlbb           )
+            CALL iom_rstput( kt, nitrst, nummldw, 'smlbn'           , smlbn           )
+            CALL iom_rstput( kt, nitrst, nummldw, 'smlatfb'         , smlatfb         )
          ELSE
-            CALL restput( inumwrs_mld, 'rmldbn'          , jpi, jpj,      1, 0, rmldbn          )
+            CALL iom_rstput( kt, nitrst, nummldw, 'rmldbn'          , rmldbn          )
 
             !-- Temperature
-            CALL restput( inumwrs_mld, 'tmlbn'           , jpi, jpj,      1, 0, tmlbn           )
-            CALL restput( inumwrs_mld, 'tml_sumb'        , jpi, jpj,      1, 0, tml_sumb        )
-            CALL restput( inumwrs_mld, 'tmltrd_csum_ub'  , jpi, jpj, jpltrd, 0, tmltrd_csum_ub  )
-            CALL restput( inumwrs_mld, 'tmltrd_atf_sumb' , jpi, jpj,      1, 0, tmltrd_atf_sumb )
+            CALL iom_rstput( kt, nitrst, nummldw, 'tmlbn'           , tmlbn           )
+            CALL iom_rstput( kt, nitrst, nummldw, 'tml_sumb'        , tml_sumb        )
+            DO jk = 1, jpltrd
+               IF( jk < 10 )   THEN
+                  WRITE(charout,FMT="('tmltrd_csum_ub_', I1)") jk
+               ELSE
+                  WRITE(charout,FMT="('tmltrd_csum_ub_', I2)") jk
+               ENDIF
+               CALL iom_rstput( kt, nitrst, nummldw, charout,  tmltrd_csum_ub(:,:,jk) )
+            ENDDO
+            CALL iom_rstput( kt, nitrst, nummldw, 'tmltrd_atf_sumb' , tmltrd_atf_sumb )
 
             !-- Salinity
-            CALL restput( inumwrs_mld, 'smlbn'           , jpi, jpj,      1, 0, smlbn           )
-            CALL restput( inumwrs_mld, 'sml_sumb'        , jpi, jpj,      1, 0, sml_sumb        )
-            CALL restput( inumwrs_mld, 'smltrd_csum_ub'  , jpi, jpj, jpltrd, 0, smltrd_csum_ub  )
-            CALL restput( inumwrs_mld, 'smltrd_atf_sumb' , jpi, jpj,      1, 0, smltrd_atf_sumb )
+            CALL iom_rstput( kt, nitrst, nummldw, 'smlbn'           , smlbn           )
+            CALL iom_rstput( kt, nitrst, nummldw, 'sml_sumb'        , sml_sumb        )
+            DO jk = 1, jpltrd
+               IF( jk < 10 )   THEN
+                  WRITE(charout,FMT="('smltrd_csum_ub_', I1)") jk
+               ELSE
+                  WRITE(charout,FMT="('smltrd_csum_ub_', I2)") jk
+               ENDIF
+               CALL iom_rstput( kt, nitrst, nummldw, charout , smltrd_csum_ub(:,:,jk) )
+            ENDDO
+            CALL iom_rstput( kt, nitrst, nummldw, 'smltrd_atf_sumb' , smltrd_atf_sumb )
          ENDIF
          !
-         CALL restclo( inumwrs_mld )
+         CALL iom_close( nummldw )     ! close the restart file (only at last time step)
          ! 
       ENDIF
@@ -118,12 +102,8 @@
     !! ** Purpose :   Read file for mixed-layer diagnostics restart.
     !!----------------------------------------------------------------------------
-    LOGICAL ::   llog
-    REAL(wp) ::  zlamt(jpi,jpj), zphit(jpi,jpj)
-    CHARACTER (len=8 ) ::   clvnames(30)
-    CHARACTER (len=32) ::   clname = 'restart_mld'
-    INTEGER  ::  itime, ibvar, &  
-                 inum  ! temporary logical unit
-    REAL(wp) ::   zdate0, zdt
-    REAL(wp), DIMENSION(1) :: zdept
+    INTEGER  ::  inum       ! temporary logical unit
+    !
+    CHARACTER (len=35) :: charout
+    INTEGER ::   jk         ! loop indice
     !!-----------------------------------------------------------------------------
 
@@ -134,45 +114,47 @@
     ENDIF
 
-    itime = 0
-    llog  = .FALSE.
-    zlamt(:,:) = 0.e0   ;   zphit(:,:) = 0.e0   ;   zdept(1)   = 0.e0
-    CALL restini( clname, jpi, jpj, zlamt, zphit, 1, zdept, 'NONE',   &
-       &          itime, zdate0, zdt, inum, domain_id=nidom )
-  
-    CALL ioget_vname( inum, ibvar, clvnames)
-
-    IF(lwp) THEN
-       WRITE(numout,*) 
-       WRITE(numout,*) '                    Info on the MLD restart file read : '
-       WRITE(numout,*) '                      File name           : ', clname
-       WRITE(numout,*) '                      number of variables : ', ibvar
-       WRITE(numout,*) '                      NetCDF variables    : ', clvnames
-       WRITE(numout,*)
-    ENDIF
+    inum = 10
+    CALL iom_open( 'restart_mld', inum )                       ! Open
 
     IF( ln_trdmld_instant ) THEN 
-       CALL restget( inum, 'tmlbb'           , jpi, jpj,      1, 0, llog, tmlbb          )
-       CALL restget( inum, 'tmlbn'           , jpi, jpj,      1, 0, llog, tmlbn          )
-       CALL restget( inum, 'tmlatfb'         , jpi, jpj,      1, 0, llog, tmlatfb        )
+       !-- Temperature
+       CALL iom_get( inum, jpdom_local, 'tmlbb'           , tmlbb          )
+       CALL iom_get( inum, jpdom_local, 'tmlbn'           , tmlbn          )
+       CALL iom_get( inum, jpdom_local, 'tmlatfb'         , tmlatfb        )
 
-       CALL restget( inum, 'smlbb'           , jpi, jpj,      1, 0, llog, smlbb          )
-       CALL restget( inum, 'smlbn'           , jpi, jpj,      1, 0, llog, smlbn          )
-       CALL restget( inum, 'smlatfb'         , jpi, jpj,      1, 0, llog, smlatfb        )
+       !-- Salinity
+       CALL iom_get( inum, jpdom_local, 'smlbb'           , smlbb          )
+       CALL iom_get( inum, jpdom_local, 'smlbn'           , smlbn          )
+       CALL iom_get( inum, jpdom_local, 'smlatfb'         , smlatfb        )
     ELSE
-       CALL restget( inum, 'rmldbn'          , jpi, jpj,      1, 0, llog, rmldbn         ) ! needed for rmld_sum
+       CALL iom_get( inum, jpdom_local, 'rmldbn'          , rmldbn         ) ! needed for rmld_sum
 
        !-- Temperature
-       CALL restget( inum, 'tmlbn'           , jpi, jpj,      1, 0, llog, tmlbn          ) ! needed for tml_sum
-       CALL restget( inum, 'tml_sumb'        , jpi, jpj,      1, 0, llog, tml_sumb       )
-       CALL restget( inum, 'tmltrd_csum_ub'  , jpi, jpj, jpltrd, 0, llog, tmltrd_csum_ub )
-       CALL restget( inum, 'tmltrd_atf_sumb' , jpi, jpj,      1, 0, llog, tmltrd_atf_sumb)
+       CALL iom_get( inum, jpdom_local, 'tmlbn'           , tmlbn          ) ! needed for tml_sum
+       CALL iom_get( inum, jpdom_local, 'tml_sumb'        , tml_sumb       )
+       DO jk = 1, jpltrd
+          IF( jk < 10 )   THEN
+             WRITE(charout,FMT="('tmltrd_csum_ub_', I1)") jk
+          ELSE
+             WRITE(charout,FMT="('tmltrd_csum_ub_', I2)") jk
+          ENDIF
+          CALL iom_get( inum, jpdom_local, charout, tmltrd_csum_ub(:,:,jk) )
+       ENDDO
+       CALL iom_get( inum, jpdom_local, 'tmltrd_atf_sumb' , tmltrd_atf_sumb)
 
        !-- Salinity
-       CALL restget( inum, 'smlbn'           , jpi, jpj,      1, 0, llog, smlbn          ) ! needed for sml_sum
-       CALL restget( inum, 'sml_sumb'        , jpi, jpj,      1, 0, llog, sml_sumb       )
-       CALL restget( inum, 'smltrd_csum_ub'  , jpi, jpj, jpltrd, 0, llog, smltrd_csum_ub )
-       CALL restget( inum, 'smltrd_atf_sumb' , jpi, jpj,      1, 0, llog, smltrd_atf_sumb)
+       CALL iom_get( inum, jpdom_local, 'smlbn'           , smlbn          ) ! needed for sml_sum
+       CALL iom_get( inum, jpdom_local, 'sml_sumb'        , sml_sumb       )
+       DO jk = 1, jpltrd
+          IF( jk < 10 )   THEN
+             WRITE(charout,FMT="('smltrd_csum_ub_', I1)") jk
+          ELSE
+             WRITE(charout,FMT="('smltrd_csum_ub_', I2)") jk
+          ENDIF
+          CALL iom_get( inum, jpdom_local, charout, smltrd_csum_ub(:,:,jk) )
+       ENDDO
+       CALL iom_get( inum, jpdom_local, 'smltrd_atf_sumb' , smltrd_atf_sumb)
 
-       CALL restclo( inum )
+       CALL iom_close( inum )
     ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRD/trdmod.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 507)
+++ trunk/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 521)
@@ -59,5 +59,5 @@
       INTEGER ::   ji, ikbu, ikbum1
       INTEGER ::   jj, ikbv, ikbvm1
-      CHARACTER(len=3) ::   clpas                                  ! number of passage
+      CHARACTER(len=3) ::   ccpas                                  ! number of passage
       REAL(wp) ::   zua, zva                                       ! scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   ztswu, ztswv               ! 2D workspace
@@ -69,6 +69,6 @@
 
       ! Control of optional arguments
-      clpas = 'fst'
-      IF( PRESENT(cnbpas) )  clpas = cnbpas
+      ccpas = 'fst'
+      IF( PRESENT(cnbpas) )  ccpas = cnbpas
 
       IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
@@ -97,7 +97,7 @@
             CASE ( jptra_trd_yad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_yad, ctype )   ! y- horiz adv
             CASE ( jptra_trd_zad )                                                         ! z- vertical adv 
-               CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas )   
+               CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, ccpas )   
                ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)                                                                   
-               z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)   
+               z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
                z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
                CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )                             ! 1st z- vertical adv 
@@ -214,8 +214,5 @@
          !                  - and the iso-neutral diffusion if activated 
          ! jptra_trd_zdf : called by trazdf.F90
-         !                 * in case of purely vertical diffusion (and not iso-neutral),
-         !                   we do not need to store the corresponding trend here, since it
-         !                   is recomputed later (at the basis of the ML, see trd_mld)
-         !                 * else (iso-neutral case) we store the vertical diffusion component in the 
+         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the 
          !                   lateral trend including the K_z contrib, which will be removed later (see trd_mld)
          !-----------------------------------------------------------------------------------------------
@@ -228,5 +225,9 @@
          CASE ( jptra_trd_bbl )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbl, '3D' )   ! bottom boundary layer
          CASE ( jptra_trd_zdf )
-            IF( ln_traldf_iso )       CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! vertical diffusion (K_z)
+            IF( ln_traldf_iso )   THEN
+               CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! vertical diffusion (K_z)
+            ELSE
+               CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zdf, '3D' )   ! vertical diffusion (K_z)
+            ENDIF
          CASE ( jptra_trd_dmp )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_dmp, '3D' )   ! internal 3D restoring (tradmp)
          CASE ( jptra_trd_qsr )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '3D' )   ! air-sea : penetrative sol radiat
