Index: trunk/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- trunk/NEMO/TOP_SRC/trcrst.F90	(revision 1254)
+++ trunk/NEMO/TOP_SRC/trcrst.F90	(revision 1287)
@@ -24,4 +24,5 @@
    USE lib_mpp
    USE iom
+   USE daymod
    
    IMPLICIT NONE
@@ -89,7 +90,13 @@
       !! ** purpose  :   read passive tracer fields in restart files
       !!----------------------------------------------------------------------
-      INTEGER  ::   jn  
-      INTEGER  ::   iarak0
-      REAL(wp) ::   zkt, zarak0
+      INTEGER  ::  jn  
+      INTEGER  ::  iarak0
+      REAL(wp) ::  zarak0
+#if defined key_pisces 
+      INTEGER  ::  ji, jj, jk
+      REAL(wp) ::  zcaralk, zbicarb, zco3
+      REAL(wp) ::  ztmas, ztmas1
+#endif
+
       !!----------------------------------------------------------------------
 
@@ -98,43 +105,14 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
 
+      CALL iom_open( cn_trcrst_in, numrtr, kiolib = jprstlib )
+
+      ! Time domain : restart
+      ! ---------------------
+      CALL trc_rst_cal( nittrc000, 'READ' )   ! calendar
+
       IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN   ;   iarak0 = 1
       ELSE                                           ;   iarak0 = 0
       ENDIF
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '   the present run starts at the time step nit000 = ', nit000
-      IF(lwp .AND. iarak0 == 1 )   WRITE(numout,*) '   and needs previous fields for Arakawa sheme ??? '
-
-
-      ! Time domain : restart
-      ! -------------------------
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' *** passive tracer restart option'
-      SELECT CASE ( nrsttr )
-      CASE ( 0 )
-         IF(lwp) WRITE(numout,*) '    nrsttr = 0 no control of nit000'
-      CASE ( 1 )
-         IF(lwp) WRITE(numout,*) '    nrsttr = 1 we control the date of nit000'
-      CASE ( 2 )
-         IF(lwp) WRITE(numout,*) '    nrsttr = 2 the date adatrj is read in restart file'
-      CASE DEFAULT
-         IF(lwp) WRITE(numout,*) '  ===>>>> nrsttr not equal 0, 1 or 2 : no control of the date'
-         IF(lwp) WRITE(numout,*) '  =======                  ========='
-      END SELECT
-
-      CALL iom_open( cn_trcrst_in, numrtr, kiolib = jprstlib )
-
-      CALL iom_get( numrtr, 'kt'   , zkt    )
       CALL iom_get( numrtr, 'arak0', zarak0 )
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' Info on the restart file read : '
-      IF(lwp) WRITE(numout,*) '    time-step           : ', NINT( zkt    )
-      IF(lwp) WRITE(numout,*) '    arakawa option      : ', NINT( zarak0 )
-
-
-      IF( nittrc000 - NINT( zkt ) /= 1 .AND. nrsttr /= 0 )  &      ! control of date
-         &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart', &
-         &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' )
 
       IF( iarak0 /= NINT( zarak0 ) )   &                           ! Control of the scheme
@@ -142,4 +120,7 @@
          & ' it must be the same type for both restart and previous run', &
          & ' centered or euler '  )
+      IF(lwp) WRITE(numout,*)
+
+      IF(lwp) WRITE(numout,*) '    arakawa option      : ', NINT( zarak0 )
 
 
@@ -149,5 +130,5 @@
       END DO
 
-     DO jn = 1, jptra
+      DO jn = 1, jptra
          CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) ) 
       END DO
@@ -157,11 +138,34 @@
       CALL iom_get( numrtr, jpdom_autoglo, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) ) 
 #endif
-
 #if defined key_pisces
+      !
+      IF( ln_pisdmp ) CALL pis_dmp_ini  ! relaxation of some tracers
+      !
+      IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
+      ELSE
+         ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
+         ! --------------------------------------------------------
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztmas   = tmask(ji,jj,jk)
+                  ztmas1  = 1. - tmask(ji,jj,jk)
+                  zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
+                  zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
+                  zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
+!!                  write(numout,*) 'plante :',ji,jj,jk,ztmas,ztmas1,ak23(ji,jj,jk),zbicarb ,zco3
+                  hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
+               END DO
+            END DO
+         END DO
+      ENDIF
       CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) 
-      CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax', xksimax(:,:) )
-      CALL trc_rst_ini  ! Initialisation of some variables
-#endif
-
+      IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
+      ELSE
+         xksimax(:,:) = xksi(:,:)
+      ENDIF
+#endif
 #if defined key_cfc
       DO jn = jp_cfc0, jp_cfc1
@@ -169,13 +173,10 @@
       END DO
 #endif
-
 #if defined key_c14b
       CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn) , qint_c14(:,:) ) 
 #endif
-
 #if defined key_my_trc
 #endif
-
-
+      
       CALL iom_close( numrtr )
       !
@@ -191,23 +192,18 @@
       INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
       !!
-      INTEGER  :: ji, jj, jk, jn
-      !!----------------------------------------------------------------------
-
-      IF(  kt == nitrst ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
-         IF(lwp) WRITE(numout,*) '~~~~~~~'
-      ENDIF
-
-      CALL iom_rstput( kt, nitrst, numrtw, 'kt'   ,  REAL( kt, wp )  )
-      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
-         CALL iom_rstput( kt, nitrst, numrtw, 'arak0', 1. )
-      ELSE
-         CALL iom_rstput( kt, nitrst, numrtw, 'arak0', 0. )
-      ENDIF
-
-
-         ! prognostic variables
-         ! --------------------
+      INTEGER  :: jn
+      REAL(wp) :: zarak0
+      !!----------------------------------------------------------------------
+
+
+      CALL trc_rst_cal( kt, 'WRITE' )   ! calendar
+
+      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN   ;   zarak0 = 1.
+      ELSE                                           ;   zarak0 = 0.
+      ENDIF
+      CALL iom_rstput( kt, nitrst, numrtw, 'arak0', zarak0 )
+
+      ! prognostic variables
+      ! --------------------
       DO jn = 1, jptra
          CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
@@ -219,28 +215,24 @@
 
 #if defined key_lobster
-      CALL iom_rstput( kt, nitrst, numrtw, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
-      CALL iom_rstput( kt, nitrst, numrtw, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
-#endif
-
-#if defined key_pisces
-      CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
-      CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
-#endif
-
+         CALL iom_rstput( kt, nitrst, numrtw, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
+#endif
+#if defined key_pisces 
+         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
+#endif
 #if defined key_cfc
-      DO jn = jp_cfc0, jp_cfc1
-         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
-      END DO
-#endif
-
+         DO jn = jp_cfc0, jp_cfc1
+            CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
+         END DO
+#endif
 #if defined key_c14b
-      CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_c14(:,:) )
-#endif
-
+         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_c14(:,:) )
+#endif
 #if defined key_my_trc
 #endif
-
-
-       IF( kt == nitrst ) THEN
+      
+      IF( kt == nitrst ) THEN
           CALL trc_rst_stat            ! statistics
           CALL iom_close( numrtw )     ! close the restart file (only at last time step)
@@ -248,17 +240,106 @@
           lrst_trc = .FALSE.
 #endif
-       ENDIF
+      ENDIF
       !
    END SUBROUTINE trc_rst_wri
 
+   SUBROUTINE trc_rst_cal( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE trc_rst_cal  ***
+      !!
+      !!  ** Purpose : Read or write calendar in restart file:
+      !!
+      !!  WRITE(READ) mode:
+      !!       kt        : number of time step since the begining of the experiment at the
+      !!                   end of the current(previous) run
+      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
+      !!                   end of the current(previous) run (REAL -> keep fractions of day)
+      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
+      !!
+      !!   According to namelist parameter nrstdt,
+      !!       nrsttr = 0  no control on the date (nittrc000 is  arbitrary).
+      !!       nrsttr = 1  we verify that nit000 is equal to the last
+      !!                   time step of previous run + 1.
+      !!       In both those options, the  exact duration of the experiment
+      !!       since the beginning (cumulated duration of all previous restart runs)
+      !!       is not stored in the restart and is assumed to be (nit000-1)*rdt.
+      !!       This is valid is the time step has remained constant.
+      !!
+      !!       nrsttr = 2  the duration of the experiment in days (adatrj)
+      !!                    has been stored in the restart file.
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      !
+      REAL(wp) ::  zkt
+#if defined key_off_tra
+      REAL(wp) ::  zndastp
+#endif
+
+      ! Time domain : restart
+      ! ---------------------
+
+      IF( TRIM(cdrw) == 'READ' ) THEN
+         CALL iom_get ( numrtr, 'kt', zkt )   ! last time-step of previous run
+         IF(lwp) THEN
+            WRITE(numout,*) ' *** Info read in restart : '
+            WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
+            WRITE(numout,*) ' *** restart option'
+            SELECT CASE ( nrsttr )
+            CASE ( 0 )   ;   WRITE(numout,*) ' nrsttr = 0 : no control of nittrc000'
+            CASE ( 1 )   ;   WRITE(numout,*) ' nrsttr = 1 : no control the date at nit000 (use ndate0 read in the namelist)'
+            CASE ( 2 )   ;   WRITE(numout,*) ' nrsttr = 2 : calendar parameters read in restart'
+            END SELECT
+            WRITE(numout,*)
+         ENDIF
+         ! Control of date 
+         IF( nittrc000  - NINT( zkt ) /= 1 .AND.  nrsttr /= 0 )                                  &
+            &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart',                 &
+            &                  ' verify the restart file or rerun with nrsttr = 0 (namelist)' )
+#if defined key_off_tra
+         ! define ndastp and adatrj
+         IF ( nrsttr == 2 ) THEN
+            CALL iom_get( numrtr, 'ndastp', zndastp ) 
+            ndastp = NINT( zndastp )
+            CALL iom_get( numrtr, 'adatrj', adatrj  )
+         ELSE
+            ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam
+            adatrj = ( REAL( nittrc000-1, wp ) * rdttra(1) ) / rday
+            ! note this is wrong if time step has changed during run
+         ENDIF
+         !
+         IF(lwp) THEN
+           WRITE(numout,*) ' *** Info used values : '
+           WRITE(numout,*) '   date ndastp                                      : ', ndastp
+           WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
+           WRITE(numout,*)
+         ENDIF
+#endif
+
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+         !
+         IF(  kt == nitrst ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
+            IF(lwp) WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         ! calendar control
+         CALL iom_rstput( kt, nitrst, numrtw, 'kt'     , REAL( kt    , wp) )   ! time-step
+         CALL iom_rstput( kt, nitrst, numrtw, 'ndastp' , REAL( ndastp, wp) )   ! date
+         CALL iom_rstput( kt, nitrst, numrtw, 'adatrj' , adatrj            )   ! number of elapsed days since
+         !                                                                     ! the begining of the run [s]
+      ENDIF
+
+   END SUBROUTINE trc_rst_cal
+
 # if defined key_pisces 
 
-   SUBROUTINE trc_rst_ini 
-      !!----------------------------------------------------------------------
-      !!                    ***  trc_rst_ini  ***
-      !!
-      !! ** purpose  : Initialisation of some variables ( hi
-      !!----------------------------------------------------------------------
-      INTEGER  :: ji, jj, jk, jn  
+   SUBROUTINE pis_dmp_ini 
+      !!----------------------------------------------------------------------
+      !!                    ***  pis_dmp_ini  ***
+      !!
+      !! ** purpose  : Relaxation of some tracers
+      !!----------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk  
       REAL(wp) ::  &
          alkmean = 2426. ,  & ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
@@ -268,5 +349,5 @@
       
       REAL(wp) ::   zvol, ztrasum
-      REAL(wp) ::   caralk, bicarb, co3
+
 
       IF(lwp)  WRITE(numout,*)
@@ -360,20 +441,6 @@
 !     trn(:,:,:,jpnum) = trn(:,:,:,jppoc) / ( 6. * xkr_massp )
 !#endif
-      !!  Set hi (???) from  total alcalinity, borat (???), akb3 (???) and ak23 (???)
-      !!  ---------------------------------------------------------------------
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1,jpi
-               caralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / ( 1. + 1.e-8 / ( rtrn + akb3(ji,jj,jk) ) )
-               co3    = ( caralk - trn(ji,jj,jk,jpdic) ) *       tmask(ji,jj,jk)   &
-                  &   +             0.5e-3               * ( 1.- tmask(ji,jj,jk) )
-               bicarb = 2.* trn(ji,jj,jk,jpdic) - caralk
-               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * bicarb / co3 )   *       tmask(ji,jj,jk)   &
-                  &         +             1.0e-9                  * ( 1.- tmask(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
-
-   END SUBROUTINE trc_rst_ini
+
+   END SUBROUTINE pis_dmp_ini
 
 #endif
@@ -404,7 +471,7 @@
          zdiag_varmin = 0.e0
          zdiag_varmax = 0.e0
-         DO ji = 1, jpi
+         DO jk = 1, jpk
             DO jj = 1, jpj
-               DO jk = 1,jpk
+               DO ji = 1, jpi
                   zvol = cvol(ji,jj,jk)
 #  if defined key_off_degrad
