source: branches/CMIP5_IPSL/NEMO/TOP_SRC/trcrst.F90 @ 1808

MODULE trcrst
   !!======================================================================
   !!                         ***  MODULE trcrst  ***
   !! TOP :   Manage the passive tracer restart
   !!======================================================================
   !! History :    -   !  1991-03  ()  original code
   !!             1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90
   !!              -   !  2005-10 (C. Ethe) print control
   !!             2.0  !  2005-10 (C. Ethe, G. Madec) revised architecture
   !!----------------------------------------------------------------------
#if defined key_top
   !!----------------------------------------------------------------------
   !!   'key_top'                                                TOP models
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !!   trc_rst :   Restart for passive tracer
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !!   'key_top'                                                TOP models
   !!----------------------------------------------------------------------
   !!   trc_rst_opn    : open  restart file
   !!   trc_rst_read   : read  restart file
   !!   trc_rst_wri    : write restart file
   !!----------------------------------------------------------------------
   USE oce_trc
   USE trc
   USE trctrp_lec
   USE lib_mpp
   USE iom
   USE trcrst_cfc      ! CFC      
   USE trcrst_lobster  ! LOBSTER  restart
   USE trcrst_pisces   ! PISCES   restart
   USE trcrst_c14b     ! C14 bomb restart
   USE trcrst_my_trc   ! MY_TRC   restart

   IMPLICIT NONE
   PRIVATE

   PUBLIC   trc_rst_opn       ! called by ???
   PUBLIC   trc_rst_read      ! called by ???
   PUBLIC   trc_rst_wri       ! called by ???

   INTEGER, PUBLIC ::   numrtr, numrtw   !: logical unit for trc restart (read and write)

   !! * Substitutions
#  include "top_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005) 
   !! $Id$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
   
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 ! defined key_off_tra
      IF( kt == nit000 ) lrst_trc = .FALSE. 
# else
      IF( kt == nit000 ) THEN
        lrst_trc = .FALSE. 
        nitrst = nitend  
      ENDIF

      IF( MOD( kt - 1, nstock ) == 0 ) THEN
         ! we use kt - 1 and not kt - nit000 to keep the same periodicity from the beginning of the experiment
         nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
         IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
      ENDIF
# endif
     ! to get better performances with NetCDF format:
     ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*ndttrc + 1)
     ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*ndttrc + 1
     IF( kt == nitrst - 2*ndttrc + 1 .OR. nstock == ndttrc .OR. ( kt == nitend - ndttrc + 1 .AND. .NOT. lrst_trc ) ) THEN
         ! 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))//"_"//TRIM(cn_trcrst_out)
         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
      !!----------------------------------------------------------------------
      !!                    ***  trc_rst_opn  ***
      !!
      !! ** purpose  :   read passive tracer fields in restart files
      !!----------------------------------------------------------------------
      INTEGER  ::  jn     
      INTEGER  ::  iarak0 
      REAL(wp) ::  zarak0
      INTEGER  ::  jlibalt = jprstlib
      LOGICAL  ::  llok

      !!----------------------------------------------------------------------

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'trc_rst_read : read the TOP restart file'
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'

      IF ( jprstlib == jprstdimg ) THEN
        ! eventually read netcdf file (monobloc)  for restarting on different number of processors
        ! if {cn_trcrst_in}.nc exists, then set jlibalt to jpnf90 
        INQUIRE( FILE = TRIM(cn_trcrst_in)//'.nc', EXIST = llok )
        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF 
      ENDIF

      CALL iom_open( cn_trcrst_in, numrtr, kiolib = jlibalt ) 

      ! 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
      CALL iom_get( numrtr, 'arak0', zarak0 )

      IF( iarak0 /= NINT( zarak0 ) )   &                           ! Control of the scheme
         & CALL ctl_stop( ' ===>>>> : problem with advection scheme', &
         & ' 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 )

      ! READ prognostic variables and computes diagnostic variable
      DO jn = 1, jptra
         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
      END DO

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

      IF( lk_lobster )   CALL trc_rst_read_lobster( numrtr )      ! LOBSTER bio-model
      IF( lk_pisces  )   CALL trc_rst_read_pisces ( numrtr )      ! PISCES  bio-model
      IF( lk_cfc     )   CALL trc_rst_read_cfc    ( numrtr )      ! CFC     tracers
      IF( lk_c14b    )   CALL trc_rst_read_c14b   ( numrtr )      ! C14 bomb  tracer
      IF( lk_my_trc  )   CALL trc_rst_read_my_trc ( numrtr )      ! MY_TRC  tracers

      CALL iom_close( numrtr )
      !
   END SUBROUTINE trc_rst_read

   SUBROUTINE trc_rst_wri( kt )
      !!----------------------------------------------------------------------
      !!                    ***  trc_rst_wri  ***
      !!
      !! ** purpose  :   write passive tracer fields in restart files
      !!----------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
      !!
      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) )
      END DO

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

      IF( lk_lobster )   CALL trc_rst_wri_lobster( kt, nitrst, numrtw )      ! LOBSTER bio-model
      IF( lk_pisces  )   CALL trc_rst_wri_pisces ( kt, nitrst, numrtw )      ! PISCES  bio-model
      IF( lk_cfc     )   CALL trc_rst_wri_cfc    ( kt, nitrst, numrtw )      ! CFC     tracers
      IF( lk_c14b    )   CALL trc_rst_wri_c14b   ( kt, nitrst, numrtw )      ! C14 bomb  tracer
      IF( lk_my_trc  )   CALL trc_rst_wri_my_trc ( kt, nitrst, numrtw )      ! MY_TRC  tracers

      IF( kt == nitrst ) THEN
          CALL trc_rst_stat            ! statistics
          CALL iom_close( numrtw )     ! close the restart file (only at last time step)
#if ! defined key_trdmld_trc
          lrst_trc = .FALSE.
#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


   SUBROUTINE trc_rst_stat
      !!----------------------------------------------------------------------
      !!                    ***  trc_rst_stat  ***
      !!
      !! ** purpose  :   Compute tracers statistics
      !!----------------------------------------------------------------------

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


      IF( lwp ) THEN
         WRITE(numout,*) 
         WRITE(numout,*) '           ----TRACER STAT----             '
         WRITE(numout,*) 
      ENDIF
      
      zdiag_tot = 0.e0
      DO jn = 1, jptra
         zdiag_var    = 0.e0
         zdiag_varmin = 0.e0
         zdiag_varmax = 0.e0
         DO jk = 1, jpk
            DO jj = 1, jpj
               DO ji = 1, jpi
                  zvol = cvol(ji,jj,jk)
#  if defined key_off_degrad
                  zvol = zvol * facvol(ji,jj,jk)
#  endif
                  zdiag_var = zdiag_var + trn(ji,jj,jk,jn) * zvol
               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
      
      zder = ( ( zdiag_tot - trai ) / ( trai + 1.e-12 )  ) * 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, ' %'
      
   END SUBROUTINE trc_rst_stat

#else
   !!----------------------------------------------------------------------
   !!  Dummy module :                                     No passive tracer
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE trc_rst_read                      ! Empty routines
   END SUBROUTINE trc_rst_read
   SUBROUTINE trc_rst_wri( kt )
      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