MODULE trcrst
   !!======================================================================
   !!
   !!                       *** MODULE trcrst ***
   !!
   !!   Read the restart files for passive tracers
   !!
   !!======================================================================
   !!  TOP 1.0,  LOCEAN-IPSL (2005) 
   !! $Id$
   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
   !!----------------------------------------------------------------------
#if defined key_passivetrc   
   !!----------------------------------------------------------------------
   !! * Modules used
   !! ==============
   USE oce_trc
   USE trc
   USE sms
   USE trctrp_lec   
   USE lib_mpp
   USE iom
   
   IMPLICIT NONE
   PRIVATE
   
   !! * Accessibility
   PUBLIC trc_rst_opn
   PUBLIC trc_rst_read
   PUBLIC trc_rst_wri
   
   !! * Module variables
   LOGICAL, PUBLIC ::   lrst_trc         !: logical to control the trc restart write 
   INTEGER, PUBLIC ::   numrtr, numrtw   !: logical unit for trc restart (read and write)


   !! * Substitutions
#  include "passivetrc_substitute.h90"
   
CONTAINS
   
   SUBROUTINE trc_rst_opn( kt )
      !!----------------------------------------------------------------------
      !!                    ***  trc_rst_opn  ***
      !!
      !! ** purpose  :   output of sea-trc variable in a netcdf file
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt       ! number of iteration
      !
      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
      CHARACTER(LEN=50)   ::   clname   ! trc output restart file name
      !!----------------------------------------------------------------------
      !

      IF( kt == nit000 )  THEN
         lrst_trc = .FALSE.
#if defined key_off_tra
         nitrst = nitend  ! in online version, already done in rst_opn routine defined in restart.F90 module
#endif
      ENDIF
      
      IF( kt == nitrst - ndttrc .OR. nitend - nit000 + 1 < 2 * ndttrc ) THEN
         ! beware if model runs less than 2*ndttrc time step
         ! beware of the format used to write kt (default is i8.8, that should be large enough)
         IF( nitrst > 1.0e9 ) THEN   
            WRITE(clkt,*) nitrst
         ELSE
            WRITE(clkt,'(i8.8)') nitrst
         ENDIF
         ! create the file
         IF(lwp) WRITE(numout,*)
         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_trc"
         IF(lwp) WRITE(numout,*) '             open trc restart.output NetCDF file: '//clname
         CALL iom_open( clname, numrtw, ldwrt = .TRUE., kiolib = jprstlib )
         lrst_trc = .TRUE.
      ENDIF
      !
   END SUBROUTINE trc_rst_opn


   SUBROUTINE trc_rst_read 
      !!===========================================================================================
      !!
      !!                       ROUTINE trc_rst_read
      !!                       *******************
      !!
      !!  PURPOSE :
      !!  ---------
      !!     READ files for restart for passive tracer
      !!
      !!   METHOD :
      !!   -------
      !!      READ the previous fields on the FILE nutrst
      !!      the first record indicates previous characterics
      !!      after control with the present run, we READ :
      !!      - prognostic variables on the second and more record
      !!
      !!   History:
      !!   --------
      !!  original  : 96-11
      !!  00-05 (A. Estublier) TVD Limiter Scheme key_trc_tvd 
      !!  00-12 (O. Aumont, E. Kestenare) read restart file for sediments 
      !!  01-05 (O. Aumont, E. Kestenare) read restart file for calcite and silicate sediments 
      !!  05-03 (O. Aumont and A. El Moussaoui) F90           
      !!------------------------------------------------------------------------
      INTEGER ::  ji, jj, jk, jn  
      INTEGER ::  iarak0    
      REAL(wp) :: zkt, zarak0
      REAL(wp) :: caralk, bicarb, co3

#if defined key_trc_pisces 
#   if ! defined key_cfg_1d && ( defined key_orca_r4 || defined key_orca_r2 || defined key_orca_r05 || defined key_orca_r025 )
      REAL(wp) ::  ztrasum
#   endif
#endif

      !!---------------------------------------------------------------------
      !!  OPA.9 03-2005  
      !!---------------------------------------------------------------------
      !! 0. initialisations
      !!------------------


      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
         iarak0 = 1
      ELSE
         iarak0 = 0
      ENDIF


      IF(lwp) WRITE(numout,*) ' '
      IF(lwp) WRITE(numout,*) ' *** trc_rst beginning of restart for'
      IF(lwp) WRITE(numout,*) ' passive tracer'
      IF(lwp) WRITE(numout,*) ' the present run :'
      IF(lwp) WRITE(numout,*) '   with the time nit000 : ',nit000
      IF(lwp) THEN
         IF( iarak0 == 1 ) THEN
            WRITE(numout,*) '   and before fields for Arakawa sheme '
         ENDIF
         WRITE(numout,*) ' '
      ENDIF

      ! Time domain : restart
      ! -------------------------

      IF(lwp) WRITE(numout,*)
      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 ( 'restart.trc', 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(lwp) WRITE(numout,*)


      !! 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 nrstdt = 0 (namelist)' )

      !! Control of the scheme
      !! ------------------------

      IF( iarak0 /= NINT( zarak0 ) ) &
           & CALL ctl_stop( ' ===>>>> : problem with advection scheme', &
           & ' it must be the same type for both restart and previous run', &
           & ' centered or euler '  )


      !! ... READ prognostic variables and computes diagnostic variable
      !! ---------------------------------------------------------------

      DO jn = 1, jptra
         CALL iom_get( numrtr, jpdom_local, 'TRN'//ctrcnm(jn), trn(:,:,:,jn)   ) 
      END DO

      DO jn = 1, jptra
         CALL iom_get( numrtr, jpdom_local, 'TRB'//ctrcnm(jn), trb(:,:,:,jn)   ) 
      END DO

#if defined key_trc_lobster1
      CALL iom_get( numrtr, jpdom_local, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) ) 
      CALL iom_get( numrtr, jpdom_local, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) ) 

#elif defined key_trc_pisces
      CALL iom_get( numrtr, jpdom_local, 'Silicalim', xksi(:,:) ) 
      xksimax = xksi

      CALL iom_get( numrtr, jpdom_local, 'SED'//ctrcnm(jppoc), sedpoc(:,:) ) 
      CALL iom_get( numrtr, jpdom_local, 'SED'//ctrcnm(jpcal), sedcal(:,:) ) 
      CALL iom_get( numrtr, jpdom_local, 'SED'//ctrcnm(jpsil), sedsil(:,:) ) 

#elif defined key_cfc
      DO jn = 1, jptra
         CALL iom_get( numrtr, jpdom_local, 'qint'//ctrcnm(jn),qint(:,:,jn)) 
      END DO
      DO jn = 1, jptra
         CALL iom_get( numrtr, jpdom_local, 'qtr'//ctrcnm(jn) ,qtr( :,:,jn)) 
      END DO
#endif


#if defined key_trc_pisces 

#if ! defined key_cfg_1d && ( defined key_orca_r4 || defined key_orca_r2 || defined key_orca_r05 || defined key_orca_r025 ) 

      ztrasum = 0.
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               ztrasum = ztrasum + trn(ji,jj,jk,jptal) * tmask(ji,jj,jk) * tmask_i(ji,jj)  &
#if defined key_off_degrad
                  &              * facvol(ji,jj,jk)   &
#endif

                  &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
            END DO
         END DO
      END DO

      IF( lk_mpp ) THEN 
         CALL mpp_sum( ztrasum )     ! sum over the global domain  
      END IF

      WRITE(0,*) 'TALK moyen ', ztrasum/areatot*1E6
      ztrasum = ztrasum/areatot*1E6
      trn(:,:,:,jptal) = trn(:,:,:,jptal)*2391./ztrasum

      ztrasum = 0.
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               ztrasum = ztrasum + trn(ji,jj,jk,jppo4) * tmask(ji,jj,jk) * tmask_i(ji,jj)  &
#if defined key_off_degrad
                  &              * facvol(ji,jj,jk)   &
#endif

                  &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
            END DO
         END DO
      END DO

      IF( lk_mpp ) THEN 
         CALL mpp_sum( ztrasum )     ! sum over the global domain  
      END IF


      WRITE(0,*) 'PO4 moyen ', ztrasum/areatot*1E6/122.
      ztrasum = ztrasum/areatot*1E6/122.
      trn(:,:,:,jppo4) = trn(:,:,:,jppo4)*2.165/ztrasum

      ztrasum = 0.
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               ztrasum = ztrasum + trn(ji,jj,jk,jpno3) * tmask(ji,jj,jk) * tmask_i(ji,jj)  &
#if defined key_off_degrad
                  &              * facvol(ji,jj,jk)   &
#endif

                  &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
            END DO
         END DO
      END DO

      IF( lk_mpp ) THEN 
         CALL mpp_sum( ztrasum )     ! sum over the global domain  
      END IF


      WRITE(0,*) 'NO3 moyen ', ztrasum/areatot*1E6/7.6
      ztrasum = ztrasum/areatot*1E6/7.6
      trn(:,:,:,jpno3) = trn(:,:,:,jpno3)*30.9/ztrasum

      ztrasum = 0.
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               ztrasum = ztrasum + trn(ji,jj,jk,jpsil) * tmask(ji,jj,jk) * tmask_i(ji,jj)   &
#if defined key_off_degrad
                  &              * facvol(ji,jj,jk)   &
#endif

                  &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
            END DO
         END DO
      END DO

      IF( lk_mpp ) THEN 
         CALL mpp_sum( ztrasum )     ! sum over the global domain  
      END IF

      WRITE(0,*) 'SiO3 moyen ', ztrasum/areatot*1E6
      ztrasum = ztrasum/areatot*1E6
      trn(:,:,:,jpsil) = MIN( 400E-6,trn(:,:,:,jpsil)*91.51/ztrasum) 

#endif

!#if defined key_trc_kriest
!      !! Initialize number of particles from a standart restart file
!      !! The name of big organic particles jpgoc has been only change
!      !! and replace by jpnum but the values here are concentration
!      trn(:,:,:,jppoc) = trn(:,:,:,jppoc) + trn(:,:,:,jpnum)
!      trn(:,:,:,jpnum) = trn(:,:,:,jppoc) / ( 6. * xkr_massp )
!#endif
      !!  Initialization of chemical variables of the carbon cycle
      !!  --------------------------------------------------------
      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)   &
                  &        +(1.-tmask(ji,jj,jk))*.5e-3
               bicarb = (2.*trn(ji,jj,jk,jpdic)-caralk)
               hi(ji,jj,jk) = (ak23(ji,jj,jk)*bicarb/co3)     &
                  &  *tmask(ji,jj,jk)+(1.-tmask(ji,jj,jk))*1.e-9
            ENDDO
         ENDDO
      ENDDO
#endif
      trb(:,:,:,:) = trn(:,:,:,:)

      CALL iom_close( numrtr )


   END SUBROUTINE trc_rst_read

   SUBROUTINE trc_rst_wri(kt)
      !! ==================================================================================
      !!
      !!                       ROUTINE trc_rst_wri
      !!                       ******************
      !!
      !!  PURPOSE :
      !!  ---------
      !!     WRITE restart fields in nutwrs
      !!   METHOD :
      !!   -------
      !!
      !!   nutwrs FILE:
      !!   each nstock time step , SAVE fields which are necessary for
      !!   passive tracer restart
      !!
      !!
      !!   INPUT :
      !!   -----
      !!      argument
      !!              kt              : time step
      !!      COMMON
      !!            /cottrc/          : passive tracers fields (before,now
      !!                                  ,after)
      !!
      !!   OUTPUT :
      !!   ------
      !!      FILE
      !!           nutwrs          : standard restart fields OUTPUT
      !!
      !!   WORKSPACE :
      !!   ---------
      !!      ji,jj,jk,jn
      !!
      !!   History:
      !!   --------
      !!      original : 96-12
      !!      addition : 99-12 (M.-A. Foujols) NetCDF FORMAT with ioipsl
      !!      additions : 00-05 (A. Estublier)
      !!                  TVD Limiter Scheme : key_trc_tvd
      !!      additions : 01-01 (M.A Foujols, E. Kestenare) bug fix: restclo
      !!      additions : 01-01 (O. Aumont, E. Kestenare)
      !!                  write restart file for sediments
      !!      additions : 01-05 (O. Aumont, E. Kestenare)
      !!                  write restart file for calcite and silicate sediments
      !!   05-03 (O. Aumont and A. El Moussaoui) F90
      !!========================================================================================!

      !! * Arguments
      !! -----------
      INTEGER, INTENT( in ) :: kt

      !! * local declarations
      !! ====================

      INTEGER  :: ji,jj,jk,jn
      REAL(wp) :: zdiag_var, zdiag_varmin, zdiag_varmax, zdiag_tot
      REAL(wp) :: zder


      !! 1. OUTPUT of restart fields (nutwrs)
      !! ---------------------------

      IF( MOD(kt,nstock) == 0 .OR. kt == nitend ) THEN

         !! 0. initialisations
         !! ------------------

         IF(lwp) WRITE(numout,*) ' '
         IF(lwp) WRITE(numout,*) 'trc_wri : write the passive tracer restart file in NetCDF format ',   &
            'at it= ',kt,' date= ',ndastp
         IF(lwp) WRITE(numout,*) '~~~~~~~~~'


         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
         ! --------------------

         DO jn=1,jptra
            CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
         ENDDO

         DO jn=1,jptra
            CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
         END DO

#if defined key_trc_lobster1
         CALL iom_rstput( kt, nitrst, numrtw, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
         CALL iom_rstput( kt, nitrst, numrtw, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
#elif defined key_trc_pisces
         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )

         CALL iom_rstput( kt, nitrst, numrtw, 'SED'//ctrcnm(jppoc), sedpoc(:,:) )
         CALL iom_rstput( kt, nitrst, numrtw, 'SED'//ctrcnm(jpcal), sedcal(:,:) )
         CALL iom_rstput( kt, nitrst, numrtw, 'SED'//ctrcnm(jpsil), sedsil(:,:) )
         
#elif defined key_cfc
         DO jn=1,jptra
            CALL iom_rstput( kt, nitrst, numrtw, 'qint'//ctrcnm(jn), qint(:,:,jn) )
         END DO
         DO jn=1,jptra
            CALL iom_rstput( kt, nitrst, numrtw, 'qtr'//ctrcnm(jn) , qtr( :,:,jn) )
         END DO
#endif

         IF (lwp) WRITE(numout,*) '----TRACER STAT----'

         zdiag_tot=0.
         DO jn=1,jptra
            zdiag_var=0.
            zdiag_varmin=0.
            zdiag_varmax=0.

            DO ji=1, jpi
               DO jj=1, jpj
                  DO jk=1,jpk
                     zdiag_var=zdiag_var+trn(ji,jj,jk,jn)*tmask(ji,jj,jk)*tmask_i(ji,jj)  &
#if defined key_off_degrad
                        &   * facvol(ji,jj,jk)   &
#endif
                        &   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)

                  END DO
               END DO
            END DO

            zdiag_varmin=MINVAL(trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)))
            zdiag_varmax=MAXVAL(trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)))

            IF( lk_mpp ) THEN
               CALL mpp_min(zdiag_varmin)      ! min over the global domain
               CALL mpp_max(zdiag_varmax)      ! max over the global domain
               CALL mpp_sum(zdiag_var)         ! sum over the global domain
            END IF

            zdiag_tot=zdiag_tot+zdiag_var
            zdiag_var=zdiag_var/areatot

            IF (lwp) WRITE(numout,*) 'MEAN NO ',jn,ctrcnm(jn),' =',zdiag_var,'MIN= '  &
               ,zdiag_varmin,'MAX= ',zdiag_varmax

         END DO

         zdiag_tot=zdiag_tot
         zder=((zdiag_tot-trai)/trai)*100._wp
         IF (lwp) WRITE(numout,*) 'Integral of all tracers over the full domain  =',zdiag_tot
         IF (lwp) WRITE(numout,*) 'Drift of the sum of all tracers =',zder, '%'

         CALL iom_close(numrtw)

      ENDIF

   END SUBROUTINE trc_rst_wri


#else
   !!======================================================================
   !!  Empty module : No passive tracer
   !!======================================================================
CONTAINS

   SUBROUTINE trc_rst_read
      !! no passive tracers
   END SUBROUTINE trc_rst_read

   SUBROUTINE trc_rst_wri(kt)
      !! no passive tracers
      INTEGER, INTENT ( in ) :: kt
      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
   END SUBROUTINE trc_rst_wri
   
#endif
   
END MODULE trcrst