!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! - Ariane - (May - 2007)
!! 
!! bruno.blanke@univ-brest.fr and nicolas.grima@univ-brest.fr
!! 
!! This software is a computer program whose purpose is 
!! the computation of 3D streamlines in a given velocity field 
!! (as the output of an Ocean General Circulation Model) and 
!! subsequent water masses analyses.
!! 
!! This software is governed by the CeCILL license under French law and
!! abiding by the rules of distribution of free software.  You can  use, 
!! modify and/ or redistribute the software under the terms of the CeCILL
!! license as circulated by CEA, CNRS and INRIA at the following URL
!! "http://www.cecill.info". 
!! 
!! As a counterpart to the access to the source code and  rights to copy,
!! modify and redistribute granted by the license, users are provided only
!! with a limited warranty  and the software's author,  the holder of the
!! economic rights,  and the successive licensors  have only  limited
!! liability. 
!! 
!! In this respect, the user's attention is drawn to the risks associated
!! with loading,  using,  modifying and/or developing or reproducing the
!! software by the user in light of its specific status of free software,
!! that may mean  that it is complicated to manipulate,  and  that  also
!! therefore means  that it is reserved for developers  and  experienced
!! professionals having in-depth computer knowledge. Users are therefore
!! encouraged to load and test the software's suitability as regards their
!! requirements in conditions enabling the security of their systems and/or 
!! data to be ensured and,  more generally, to use and operate it in the 
!! same conditions as regards security. 
!! 
!! The fact that you are presently reading this means that you have had
!! knowledge of the CeCILL license and that you accept its terms.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!****h* ariane/mod_input_data
!! NAME
!!   mod_input_data (mod_input_data.f90 - Fortran90 module)
!!
!! USAGE
!!   Include 'USE mod_input_data' in the header of your Fortran 90 source 
!!   code.
!!   Then you'll have access to the subroutine:
!!      - sub_input_data
!!      - sub_input_data_main
!!      - sub_transp_alloc
!!      - sub_transp_dealloc
!!      - sub_tracer_alloc
!!      - sub_tracer_dealloc
!!
!! FUNCTION
!!   Read input data: tracers and currents (compute transport).
!!   File format is netcdf (netcdf version 3.6.0 or newer is required).
!!
!! AUTHOR
!!   * Origin  : Nicolas Grima (April-May 2005)
!! 
!! CREATION DATE
!!   * April-May 2005
!!
!! HISTORY
!!   Date (dd/mm/yyyy/) - Modification(s)
!!
!! RESULT
!!   
!!
!! EXAMPLES
!!   * USE mod_input_data
!!
!! NOTES
!!   ROBODoc header style.
!!
!! TODO
!!   
!!
!! PORTABILITY
!!         Machine-OS    - Fortran90/95 compiler
!!   * i686-pc-linux-gnu -         ifort
!!   * i686-pc-linux-gnu -          g95
!!
!! SEE ALSO
!!   
!!
!! USES
!!   * USE mod_precision
!!   * USE mod_namelist
!!   * USE mod_input_grid
!!   * USE mod_w_comput
!!   * USE mod_rhostp
!!   * USE mod_netcdf
!!   * USE reducmem
!!
!! USED BY
!!   * posini
!!   * trajec
!!
!! SOURCE
!!=========================================================================

MODULE mod_input_data 10,9

  !------------------!
  ! USE ASSOCIAITION !
  !------------------!
  USE mod_precision
  USE mod_namelist
  USE mod_input_grid
  USE mod_seq
  USE mod_w_comput
  USE mod_rhostp
  USE mod_netcdf
  USE mod_reducmem
  USE netcdf
  USE mod_B2C_grid_interpolation
  USE mod_B2C_grid_save

  !-------------!
  ! DECLARATION !
  !-------------!
  IMPLICIT NONE

  !- Global variables -!
  REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
       uu      , & ! Zonal Transport
       vv      , & ! Meridional transport
       ww          ! Vertical transport

  REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
       epr       ! + Evaporation - Precipitaiton - Runoff

  REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
       tt      , & ! Temperature
       ss      , & ! Salinity
       rr          ! Density

  !! NG: 05/08/2008 : global variables for ROMS
  !! If ROMS is not selected these arrays don't take memory...
  REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
       zeta_a, &
       zw0 , &
       e3u , &
       e3v

  REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
       tmp_array

  INTEGER(kind = iprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
       uBgridmask, &
       vBgridmask

  LOGICAL :: id_comments = .FALSE.

CONTAINS
  !!***
  !=========================================================================
  !!****f* mod_input_data/sub_input_data()
  !! NAME
  !!   sub_input_data()
  !!
  !! FUNCTION
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data() 1,6

    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)'========================================================'
    WRITE(lun_standard,*)'= READ INPUT GEOPHYSICAL DATA AND STORE THEM IN MEMORY ='
    WRITE(lun_standard,*)'========================================================'

    ! allocate uu, vv, ww
    CALL sub_transp_alloc(dims_reg(1,3),dims_reg(2,3),dims_reg(3,3),dims_reg(4,3))
    uu(:,:,:,:) = 0._rprec
    vv(:,:,:,:) = 0._rprec
    ww(:,:,:,:) = 0._rprec

    IF (key_alltracers) THEN
      ! allocate tt, ss, rr
      CALL sub_tracer_alloc(dims_reg(1,3),dims_reg(2,3),dims_reg(3,3),dims_reg(4,3))
      tt(:,:,:,:) = 0._rprec
      ss(:,:,:,:) = 0._rprec
      rr(:,:,:,:) = 0._rprec
    ENDIF

    IF (key_roms) THEN

      CALL sub_input_data_roms()

    ELSEIF (key_symphonie) THEN
      
      CALL sub_input_data_symphonie()

    ELSEIF (key_B2C_grid) THEN

       IF (B2C_grid_Z_or_Sigma=='Z') THEN
          CALL sub_input_data_B2C_gridz()
       ELSE
          STOP
          !!TODO:  CALL sub_input_data_B2CD_grids()
       ENDIF

    ELSE

      CALL sub_input_data_opa()

    ENDIF

  END SUBROUTINE sub_input_data
  !!***
  !!***
  !=========================================================================
  !!****f* mod_input_data_seq/sub_input_data_seq()
  !! NAME
  !!   sub_input_data_seq()
  !!
  !! FUNCTION
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_alloc_init_seq() 2,7

    INTEGER(kind = iprec) :: alloc_size

    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)'====================================================='
    WRITE(lun_standard,*)'= ALLOCATE INPUT GEOPHYSICAL DATA (sequential mode) ='
    WRITE(lun_standard,*)'====================================================='

    ! allocate uu, vv, ww
    CALL sub_transp_alloc(dims_reg(1,3),dims_reg(2,3),dims_reg(3,3),it_ind)
    uu(:,:,:,:) = 0._rprec
    vv(:,:,:,:) = 0._rprec
    ww(:,:,:,:) = 0._rprec

    IF (key_alltracers) THEN
       ! allocate tt, ss, rr
       CALL sub_tracer_alloc(dims_reg(1,3),dims_reg(2,3),dims_reg(3,3),it_ind)
       tt(:,:,:,:) = 0._rprec
       ss(:,:,:,:) = 0._rprec
       rr(:,:,:,:) = 0._rprec
    ENDIF

    !!NG: This is done if key_interp_temporal = .TRUE.
    !!NG: to read a first input data.
    !!NG: the last  input data time step in forward mode.
    !!NG: the first input data time step in backward mode.
    IF (key_interp_temporal) THEN

       IF (key_alltracers) THEN
          IF ((key_roms).OR.(key_symphonie))THEN
             alloc_size = 7
          ELSE
             IF (TRIM(w_surf_option) == 'E-P-R') THEN
                alloc_size = 7
             ELSE
                alloc_size = 6
             ENDIF
          ENDIF
       ELSE
          IF (TRIM(w_surf_option) == 'E-P-R') THEN
             alloc_size = 4
          ELSE
             alloc_size = 3
          ENDIF
       ENDIF

       CALL sub_seq_alloc(alloc_size)
       !!NG      CALL sub_seq_interp_init()

       IF (TRIM(forback) == 'forward' ) THEN
          forback = 'backward'
       ELSE
          forback = 'forward'
       ENDIF

       CALL sub_seq_init()

       IF (key_roms) THEN
          CALL sub_input_data_seq_roms( &
               ncids(:)               , &
               varids(:)              , &
               new_file(:)            , &
               ind_file(:)            , &
               ind_time(:)            , &
               ind_time_size(:)       , &
               sdimsorders(:,:)        )

       ELSEIF ( key_symphonie) THEN
          CALL sub_input_data_seq_symphonie( &
               ncids(:)                    , &
               varids(:)                   , &
               new_file(:)                 , &
               ind_file(:)                 , &
               ind_time(:)                 , &
               ind_time_size(:)            , &
               sdimsorders(:,:)            )

       ELSEIF (key_B2C_grid) THEN

          WRITE(lun_error,*)''
          WRITE(lun_error,*)'STOP === key_B2C_grid is not available === STOP'
          STOP

       ELSE
          CALL sub_input_data_seq_opa( &
               ncids(:)              , &
               varids(:)             , &
               new_file(:)           , &
               ind_file(:)           , &
               ind_time(:)           , &
               ind_time_size(:)      , &
               sdimsorders(:,:)       )
       ENDIF ! (key_roms)

       IF (TRIM(forback) == 'forward' ) THEN
          forback = 'backward'
       ELSE
          forback = 'forward'
       ENDIF

    ENDIF

  END SUBROUTINE sub_input_data_alloc_init_seq
  !!***

  !=========================================================================
  !!****f* mod_input_data/sub_input_data_opa()
  !! NAME
  !!   sub_input_data_opa()
  !!
  !! FUNCTION
  !!   * read OPA tracer variables.
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_opa() 1,9

    ! Variable number: U, V, [W, T, S, R]
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    INTEGER(kind  = iprec) :: kk, ll
    INTEGER(kind  = iprec) :: i,j,k,l

    !-------------!
    ! Code begins !
    !-------------!

    !======================= ZONAL TRANSPORT ==================!
    !-- Read Zonal Current --!
    CALL sub_input_data_main(c_dir_zo,c_prefix_zo,ind0_zo,indn_zo,maxsize_zo, &
         c_suffix_zo, nc_var_zo,nc_var_eivu,nc_att_mask_zo,uu(:,:,:,:))

    !-- Compute transport --!
    IF (key_partialsteps) THEN

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=1))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO
    ENDIF

    WRITE(lun_standard,*)' - Transport: max ', MAXVAL(uu), ' min ', MINVAL(uu)

    !============= MERIDIONAL TRANSPORT =======================!

    !-- Read Meridional Current --!
    CALL sub_input_data_main(c_dir_me,c_prefix_me,ind0_me,indn_me,maxsize_me, &
         c_suffix_me, nc_var_me,nc_var_eivv,nc_att_mask_me,vv(:,:,:,:))

    !-- Compute transport --
    IF (key_partialsteps) THEN

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=2))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO

    ENDIF

    WRITE(lun_standard,*)' - Transport: max ', MAXVAL(vv), ' min ', MINVAL(vv)

    !=============== VERTICAL TRANSPORT =======================!
    IF (key_computew) THEN

      !-- compute w transport from U and V transport --!
      CALL sub_w_comput(uu,vv,ww)

    ELSE

      !-- Read Vertical Current --!
      CALL sub_input_data_main(c_dir_ve,c_prefix_ve,ind0_ve,indn_ve,maxsize_ve, &
           c_suffix_ve, nc_var_ve,nc_var_eivw,nc_att_mask_ve,ww(:,:,:,:))

      !-- Compute transport --
      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ww(:,:,kk,ll) = ww(:,:,kk,ll) * e1t(:,:,1,1) * e2t(:,:,1,1)
        ENDDO
      ENDDO

    ENDIF

    !!NG: 24/07/2009 
    IF (TRIM(w_surf_option) == 'E-P-R') THEN
       !!NG: 24/07/2009 We force here W="E-P-R" at the surface.
       !!NG: 24/07/2009 Physically is the better solution,
       !!NG: 24/07/2009 but in this case the non-divergence criteria
       !!NG: 24/07/2009 is not at all respected at the sea surface. 
       !!NG: 24/07/2009 This option has to be activated
       !!NG: 24/07/2009 if you know what you are doing.

       !-- Read E-R-R --!
       CALL sub_input_data_main(c_dir_ep,c_prefix_ep,ind0_ep,indn_ep,maxsize_ep, &
           c_suffix_ep, nc_var_ep,c_none,nc_att_mask_ep,epr(:,:,:,:))

       DO ll = 1, lmt
          ww(:,:,1,ll) = epr(:,:,1,ll) * epr_coef * e1t(:,:,1,1) * e2t(:,:,1,1)
       ENDDO
       
    ELSEIF (TRIM(w_surf_option) == 'zero') THEN
       !!NG: 24/07/2009 We force here W=0 at the surface.
       !!NG: 24/07/2009 Is the same case that the precedent, but you don't
       !!NG: 24/07/2009 have the "E-P-R" fields available.
       !!NG: 24/07/2009 This option has to be activated
       !!NG: 24/07/2009 if you know what you are doing.
       ww(:,:,1,:) = 0._rprec
    ELSE
       !!NG: 24/07/2009 Nothing is changed compared to the precedent version
       !!NG: 24/07/2009 of Ariane. The non-divergence criteria is respected,
       !!NG: 24/07/2009 but a lot of particles can be intercepted at the surface.
       !!NG: 24/07/2009 To reduce the number of particles intercepted at the sea
       !!NG: 24/07/2009 surface, please activate one the two options above mentioned.
    ENDIF

    WRITE(lun_standard,*)' - Transport: max ', MAXVAL(ww), ' min ', MINVAL(ww)

    IF (key_alltracers) THEN

      !==================== TEMPERATURE =========================!
      !-- Read Temperature --!
      CALL sub_input_data_main(c_dir_te,c_prefix_te,ind0_te,indn_te,maxsize_te, &
           c_suffix_te, nc_var_te,c_none,nc_att_mask_te,tt(:,:,:,:))

      !====================== SALINITY ===========================!
      !-- Read Salinity --!
      CALL sub_input_data_main(c_dir_sa,c_prefix_sa,ind0_sa,indn_sa,maxsize_sa, &
           c_suffix_sa, nc_var_sa,c_none,nc_att_mask_sa,ss(:,:,:,:))

      !======================= DENSITY ===========================!
      IF (key_computesigma) THEN

        !-- Compute Density from Temperature and Salinity --!
        DO l = 1, lmt
          DO k = 1,  dims_reg(3,3)
            DO j = 1,  dims_reg(2,3)
              DO i = 1,  dims_reg(1,3)
                rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
              END DO
            END DO
          END DO
        ENDDO

      ELSE

        !-- Read Density --!
        CALL sub_input_data_main(c_dir_de,c_prefix_de,ind0_de,indn_de,maxsize_de, &
             c_suffix_de, nc_var_de,c_none,nc_att_mask_de,rr(:,:,:,:))

      ENDIF

    ENDIF

  END SUBROUTINE sub_input_data_opa
  !!***
  !=========================================================================
  !!****f* mod_input_data_seq/sub_input_data_seq_opa()
  !! NAME
  !!   sub_input_data_seq_opa()
  !!
  !! FUNCTION
  !!   * read OPA tracer variables.
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_seq_opa( & 4,9
       ncids                       , &
       varids                      , &
       new_file                    , &
       ind_file                    , &
       ind_time                    , &
       ind_time_size               , &
       dimsorders                  )


    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ncids
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: varids
    LOGICAL               , DIMENSION(:), INTENT(inout) :: new_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time_size
    INTEGER(kind = iprec) , DIMENSION(:,:), INTENT(inout) :: dimsorders

    ! Variable number: U, V, [W, T, S, R]
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    INTEGER(kind  = iprec) :: kk, ll
    INTEGER(kind  = iprec) :: i,j,k,l

    INTEGER(kind  = iprec) :: ind_ep
    

    !-------------!
    ! Code begins !
    !-------------!

    !======================= ZONAL TRANSPORT ==================!
    IF (key_interp_temporal) THEN
      uu(:,:,:,2) = uu(:,:,:,1) ! We are in temporal interpolation case
    ENDIF
    !-- Read Zonal Current --!
    CALL sub_input_data_seq_main(                         &
         ncids(1), varids(1)                            , &
         new_file(1), ind_file(1)                       , &
         ind_time(1), ind_time_size(1)                  , &
         dimsorders(:,1)                                , &
         c_dir_zo, c_prefix_zo, maxsize_zo, c_suffix_zo , &
         nc_var_zo, nc_var_eivu, nc_att_mask_zo,uu(:,:,:,1:1))

    !-- Compute transport --!
    IF (key_partialsteps) THEN

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=1))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO
    ENDIF

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - Transport U: max ', MAXVAL(uu), ' min ', MINVAL(uu)
    ENDIF

    !============= MERIDIONAL TRANSPORT =======================!
    IF (key_interp_temporal) THEN
      vv(:,:,:,2) = vv(:,:,:,1) ! We are in temporal interpolation case
    ENDIF
    !-- Read Meridional Current --!
    CALL sub_input_data_seq_main(                         &
         ncids(2), varids(2)                            , &
         new_file(2), ind_file(2)                       , &
         ind_time(2), ind_time_size(2)                  , &
         dimsorders(:,2)                                , &
         c_dir_me,c_prefix_me, maxsize_me, c_suffix_me  , &
         nc_var_me,nc_var_eivv,nc_att_mask_me,vv(:,:,:,1:1))

    !-- Compute transport --
    IF (key_partialsteps) THEN

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=2))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO

    ENDIF

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - Transport V: max ', MAXVAL(vv), ' min ', MINVAL(vv)
    ENDIF

    !=============== VERTICAL TRANSPORT =======================!
    IF (key_interp_temporal) THEN
      ww(:,:,:,2) = ww(:,:,:,1) ! We are in temporal interpolation case
    ENDIF

    WRITE(lun_standard,*)''

    IF (key_computew) THEN
       WRITE(lun_standard,*)'  - Vertical Transport is computing'
      !-- compute w transport from U and V transport --!
      CALL sub_w_comput(uu,vv,ww)

    ELSE

      !-- Read Vertical Current --!
      CALL sub_input_data_seq_main(                       &
           ncids(3), varids(3)                            , &
           new_file(3), ind_file(3)                       , &
           ind_time(3), ind_time_size(3)                  , &
           dimsorders(:,3)                                , &
           c_dir_ve,c_prefix_ve, maxsize_ve, c_suffix_ve  , &
           nc_var_ve,nc_var_eivw,nc_att_mask_ve,ww(:,:,:,1:1))

      !-- Compute transport --
      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ww(:,:,kk,ll) = ww(:,:,kk,ll) * e1t(:,:,1,1) * e2t(:,:,1,1)
        ENDDO
      ENDDO

    ENDIF

    !!NG: 24/07/2009 
    IF (TRIM(w_surf_option) == 'E-P-R') THEN
       !!NG: 24/07/2009 We force here W="E-P-R" at the surface.
       !!NG: 24/07/2009 Physically is the better solution,
       !!NG: 24/07/2009 but in this case the non-divergence criteria
       !!NG: 24/07/2009 is not at all respected at the sea surface. 
       !!NG: 24/07/2009 This option has to be activated
       !!NG: 24/07/2009 if you know what you are doing.
       IF (key_alltracers) THEN
          ind_ep=7
       ELSE
          ind_ep=4
       ENDIF

       !-- Read E-R-R --!
       CALL sub_input_data_seq_main(                          &
            ncids(ind_ep), varids(ind_ep)                   , &
            new_file(ind_ep), ind_file(ind_ep)              , &
            ind_time(ind_ep), ind_time_size(ind_ep)         , &
            dimsorders(:,4)                                 , &
            c_dir_ep,c_prefix_ep, maxsize_ep, c_suffix_ep   , &
            nc_var_ep,c_none,nc_att_mask_ep,epr(:,:,:,1:1))

       ww(:,:,1,1) = epr(:,:,1,1) * epr_coef * e1t(:,:,1,1) * e2t(:,:,1,1)
       
    ELSEIF (TRIM(w_surf_option) == 'zero') THEN
       !!NG: 24/07/2009 We force here W=0 at the surface.
       !!NG: 24/07/2009 Is the same case that the precedent, but you don't
       !!NG: 24/07/2009 have the "E-P-R" fields available.
       !!NG: 24/07/2009 This option has to be activated
       !!NG: 24/07/2009 if you know what you are doing.
       ww(:,:,1,1) =0._rprec
    ELSE
       !!NG: 24/07/2009 Nothing is changed compared to the precedent version
       !!NG: 24/07/2009 of Ariane. The non-divergence criteria is respected,
       !!NG: 24/07/2009 but a lot of particles can be intercepted at the surface.
       !!NG: 24/07/2009 To reduce the number of particles intercepted at the sea
       !!NG: 24/07/2009 surface, please activate one the two options above mentioned.
    ENDIF

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - Vertical Transport: max ', MAXVAL(ww), ' min ', MINVAL(ww)
    ENDIF

    IF (key_alltracers) THEN

      !==================== TEMPERATURE =========================!
      IF (key_interp_temporal) THEN
        tt(:,:,:,2) = tt(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      !-- Read Temperature --!
      CALL sub_input_data_seq_main(                       &
           ncids(4), varids(4)                            , &
           new_file(4), ind_file(4)                       , &
           ind_time(4), ind_time_size(4)                  , &
           dimsorders(:,3)                                , &
           c_dir_te, c_prefix_te, maxsize_te, c_suffix_te , &
           nc_var_te,c_none,nc_att_mask_te,tt(:,:,:,1:1))

      !====================== SALINITY ===========================!
      IF (key_interp_temporal) THEN
        ss(:,:,:,2) = ss(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      !-- Read Salinity --!
      CALL sub_input_data_seq_main(                       &
           ncids(5), varids(5)                            , &
           new_file(5), ind_file(5)                       , &
           ind_time(5), ind_time_size(5)                  , &
           dimsorders(:,5)                                , &
           c_dir_sa,c_prefix_sa, maxsize_sa, c_suffix_sa  , &
           nc_var_sa,c_none,nc_att_mask_sa,ss(:,:,:,1:1))

      !======================= DENSITY ===========================!
      IF (key_interp_temporal) THEN
        rr(:,:,:,2) = rr(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      IF (key_computesigma) THEN

        !-- Compute Density from Temperature and Salinity --!
        DO l = 1, 1
          DO k = 1,  dims_reg(3,3)
            DO j = 1,  dims_reg(2,3)
              DO i = 1,  dims_reg(1,3)
                rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
              END DO
            END DO

          END DO
        ENDDO

      ELSE

        !-- Read Density --!
        CALL sub_input_data_seq_main(                     &
             ncids(6), varids(6)                            , &
             new_file(6), ind_file(6)                       , &
             ind_time(6), ind_time_size(6)                  , &
             dimsorders(:,6)                                , &
             c_dir_de,c_prefix_de,maxsize_de, c_suffix_de   , &
             nc_var_de,c_none,nc_att_mask_de,rr(:,:,:,1:1))

      ENDIF

    ENDIF

  END SUBROUTINE sub_input_data_seq_opa
  !!***
  !=========================================================================
  !!****f* mod_input_data/sub_input_data_B2C_gridz()
  !! NAME
  !!   sub_input_data_B2C_gridz()
  !!
  !! FUNCTION
  !!   * read B-grid Model tracer variables. Depth levels are on Z (Not Sigma)
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (November 2008)
  !! 
  !! CREATION DATE
  !!   * November 2008
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_B2C_gridz() 1,17

    ! Variable number: U, V, [W, T, S, R]
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    INTEGER(kind  = iprec) :: kk, ll
    INTEGER(kind  = iprec) :: i,j,k,l

    INTEGER(kind = iprec) :: nb_i, nb_j, nb_k

    !-------------!
    ! Code begins !
    !-------------!
    nb_i=size(uu,dim=1)
    nb_j=size(uu,dim=2)
    nb_k=size(uu,dim=3)

    !======================= ZONAL TRANSPORT ==================!
    !-- Read Zonal Current on B-grid --!
    ALLOCATE(uBgridmask(nb_i,nb_j,nb_k,1))
    call sub_memory(size(uBgridmask),'i','uBgridmask','sub_input_data_B2C_gridz')

    CALL sub_input_data_main(c_dir_zo,c_prefix_zo,ind0_zo,indn_zo,maxsize_zo, &
         c_suffix_zo, nc_var_zo,nc_var_eivu,nc_att_mask_zo,uu(:,:,:,:), &
         uBgridmask(:,:,:,1:1))

    !-- Interpolate Zonal Current from B-grid to C-grid --!
    !! CALL sub_B2C_grid_interpolation(uu, e2f, e2u, 'u', uBgridmask)
    CALL sub_B2C_grid_interpolation(uu, e2f, e2u, e3f,e3t,'u', key_partialsteps)

    !-- Save U interpolated from B grid to C grid --!
    IF (key_B2C_save_data) THEN
       CALL sub_B2C_grid_save_U(uu,uBgridmask)
    ENDIF

    !-- Compute transport --!
    IF (key_partialsteps) THEN

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=1))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO
    ENDIF

    WRITE(lun_standard,*)' - Transport U: max ', MAXVAL(uu), ' min ', MINVAL(uu)

    call sub_memory(-size(uBgridmask),'i','uBgridmask','sub_input_data_B2C_gridz')
    DEALLOCATE(uBgridmask)

    !============= MERIDIONAL TRANSPORT =======================!

    !-- Read Meridional Current --!
    ALLOCATE(vBgridmask(nb_i,nb_j,nb_k,1))
    call sub_memory(size(vBgridmask),'i','vBgridmask','sub_input_data_B2C_gridz')

    CALL sub_input_data_main(c_dir_me,c_prefix_me,ind0_me,indn_me,maxsize_me, &
         c_suffix_me, nc_var_me,nc_var_eivv,nc_att_mask_me,vv(:,:,:,:), &
         vBgridmask(:,:,:,1:1))

    !-- Interpolate Zonal Current from B-grid to C-grid --!
    !! CALL sub_B2C_grid_interpolation(vv, e1f, e1v,'v', vBgridmask)
    CALL sub_B2C_grid_interpolation(vv, e1f, e1v, e3f,e3t,'v', key_partialsteps)

    !-- Save V interpolated from B grid to C grid --!
    IF (key_B2C_save_data) THEN
       CALL sub_B2C_grid_save_V(vv,vBgridmask)
    ENDIF


    !-- Compute transport --
    IF (key_partialsteps) THEN

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=2))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, lmt
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO

    ENDIF

    WRITE(lun_standard,*)' - Transport V: max ', MAXVAL(vv), ' min ', MINVAL(vv)

    call sub_memory(-size(vBgridmask),'i','vBgridmask','sub_input_data_B2C_gridz')
    DEALLOCATE(vBgridmask)

    !=============== VERTICAL TRANSPORT =======================!
    IF (key_read_w) THEN

       !-- Read Vertical Current --!
       CALL sub_input_data_main(c_dir_ve,c_prefix_ve,ind0_ve,indn_ve,maxsize_ve, &
            c_suffix_ve, nc_var_ve,nc_var_eivw,nc_att_mask_ve,ww(:,:,:,:))


       !-- Save V interpolated on from B grid to C grid --!
       IF (key_B2C_save_data) THEN
          CALL sub_B2C_grid_save_W(ww)
       ENDIF

       !-- Compute transport --
       DO ll = 1, lmt
          DO kk = 1, dims_reg(3,3)
             ww(:,:,kk,ll) = ww(:,:,kk,ll) * e1t(:,:,1,1) * e2t(:,:,1,1)
          ENDDO
       ENDDO

       WRITE(lun_standard,*)''
       WRITE(lun_standard,*)' - vertical current is read -'

    ELSE

       !-- compute w transport from U and V transport --!
       WRITE(lun_standard,*)''
       WRITE(lun_standard,*)' - vertical current is computed -'
       CALL sub_w_comput(uu,vv,ww)

    ENDIF


    WRITE(lun_standard,*)' - Transport W: max ', MAXVAL(ww), ' min ', MINVAL(ww)

    IF (key_alltracers) THEN

      !==================== TEMPERATURE =========================!
      !-- Read Temperature --!
      CALL sub_input_data_main(c_dir_te,c_prefix_te,ind0_te,indn_te,maxsize_te, &
           c_suffix_te, nc_var_te,c_none,nc_att_mask_te,tt(:,:,:,:))

      !====================== SALINITY ===========================!
      !-- Read Salinity --!
      CALL sub_input_data_main(c_dir_sa,c_prefix_sa,ind0_sa,indn_sa,maxsize_sa, &
           c_suffix_sa, nc_var_sa,c_none,nc_att_mask_sa,ss(:,:,:,:))

      !======================= DENSITY ===========================!
      IF (key_computesigma) THEN

        !-- Compute Density from Temperature and Salinity --!
        DO l = 1, lmt
          DO k = 1,  dims_reg(3,3)
            DO j = 1,  dims_reg(2,3)
              DO i = 1,  dims_reg(1,3)
                rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
              END DO
            END DO
          END DO
        ENDDO

      ELSE

        !-- Read Density --!
        CALL sub_input_data_main(c_dir_de,c_prefix_de,ind0_de,indn_de,maxsize_de, &
             c_suffix_de, nc_var_de,c_none,nc_att_mask_de,rr(:,:,:,:))

      ENDIF

    ENDIF

  END SUBROUTINE sub_input_data_B2C_gridz
  !!***
  !=========================================================================
  !!****f* mod_input_data_seq/sub_input_data_B2C_gridz_seq()
  !! NAME
  !!   ssub_input_data_B2C_gridz_seq()
  !!
  !! FUNCTION
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April 2009)
  !! 
  !! CREATION DATE
  !!   * April 2009
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   *
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * posini_seq
  !!   * trajec_seq
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_B2C_gridz_seq( & 3,15
       ncids                       , &
       varids                      , &
       new_file                    , &
       ind_file                    , &
       ind_time                    , &
       ind_time_size               , &
       dimsorders                  , &
       first_call                        )


    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ncids
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: varids
    LOGICAL               , DIMENSION(:), INTENT(inout) :: new_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time_size
    INTEGER(kind = iprec) , DIMENSION(:,:), INTENT(inout) :: dimsorders
    LOGICAL               , OPTIONAL    , INTENT(in)    :: first_call

    ! Variable number: U, V, [W, T, S, R]
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    INTEGER(kind  = iprec) :: kk, ll
    INTEGER(kind  = iprec) :: i,j,k,l
    INTEGER(kind  = iprec) :: nb_i, nb_j, nb_k

    ! To save data interpolated on C grid
    LOGICAL                :: new_file_tmp
    LOGICAL                :: close_file_tmp
    INTEGER(kind  = iprec) :: ind_file_tmp
    INTEGER(kind  = iprec) :: dimt_tmp
    
    !-------------!
    ! Code begins !
    !-------------!
    nb_i=size(uu,dim=1)
    nb_j=size(uu,dim=2)
    nb_k=size(uu,dim=3)

    !======================= ZONAL TRANSPORT ==================!
    IF (key_interp_temporal) THEN
      uu(:,:,:,2) = uu(:,:,:,1) ! We are in temporal interpolation case
    ENDIF

    !-- Read Zonal Current on B-grid --!
    IF (.NOT.ALLOCATED(uBgridmask)) THEN
       ALLOCATE(uBgridmask(nb_i,nb_j,nb_k,1))
       call sub_memory(size(uBgridmask),'i','uBgridmask','sub_input_data_B2C_gridz_seq')
    ENDIF

    new_file_tmp = new_file(1)
    
    !-- Read Zonal Current --!
    CALL sub_input_data_seq_main(                         &
         ncids(1), varids(1)                            , &
         new_file(1), ind_file(1)                       , &
         ind_time(1), ind_time_size(1)                  , &
         dimsorders(:,1)                                , &
         c_dir_zo, c_prefix_zo, maxsize_zo, c_suffix_zo , &
         nc_var_zo, nc_var_eivu, nc_att_mask_zo,uu(:,:,:,1:1), &
         uBgridmask(:,:,:,1:1))

    ind_file_tmp   = ind_file(1)
    close_file_tmp = new_file(1)
    dimt_tmp       = ind_time_size(1)

    !-- Interpolate Zonal Current from B-grid to C-grid --!
    !! CALL sub_B2C_grid_interpolation(uu, e2f, e2u, 'u', uBgridmask)
    CALL sub_B2C_grid_interpolation(uu, e2f, e2u, e3f,e3t,'u', key_partialsteps)

    !-- Save U interpolated from B grid to C grid --!
    IF (key_B2C_save_data) THEN
       IF (PRESENT(first_call)) THEN
          IF (first_call) THEN
             CALL sub_B2C_grid_save_U_seq(&
                  uu                    , &
                  new_file_tmp          , &
                  ind_file_tmp          , &
                  maxsize_zo            , &
                  dimt_tmp              , &
                  close_file_tmp          &
                  )
          ENDIF
       ENDIF
    ENDIF

    !-- Compute transport --!
    IF (key_partialsteps) THEN

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=1))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          uu(:,:,kk,ll) = uu(:,:,kk,ll) * e2u(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO
    ENDIF

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - Transport U: max ', MAXVAL(uu), ' min ', MINVAL(uu)
    ENDIF

    !============= MERIDIONAL TRANSPORT =======================!
    IF (.NOT.ALLOCATED(vBgridmask)) THEN
       ALLOCATE(vBgridmask(nb_i,nb_j,nb_k,1))
       call sub_memory(size(vBgridmask),'i','vBgridmask','sub_input_data_B2C_gridz_seq')
    ENDIF

    IF (key_interp_temporal) THEN
      vv(:,:,:,2) = vv(:,:,:,1) ! We are in temporal interpolation case
    ENDIF

    new_file_tmp = new_file(2)

    !-- Read Meridional Current --!
    CALL sub_input_data_seq_main(                         &
         ncids(2), varids(2)                            , &
         new_file(2), ind_file(2)                       , &
         ind_time(2), ind_time_size(2)                  , &
         dimsorders(:,2)                                , &
         c_dir_me,c_prefix_me, maxsize_me, c_suffix_me  , &
         nc_var_me,nc_var_eivv,nc_att_mask_me,vv(:,:,:,1:1), &
         vBgridmask(:,:,:,1:1))

    ind_file_tmp   = ind_file(2)
    close_file_tmp = new_file(2)
    dimt_tmp       = ind_time_size(2)

    !-- Interpolate Zonal Current from B-grid to C-grid --!
    !! CALL sub_B2C_grid_interpolation(vv, e1f, e1v,'v', vBgridmask)
    CALL sub_B2C_grid_interpolation(vv, e1f, e1v, e3f,e3t,'v', key_partialsteps)

    !-- Save V interpolated from B grid to C grid --!
    IF (key_B2C_save_data) THEN
       IF (PRESENT(first_call)) THEN
          IF (first_call) THEN
             CALL sub_B2C_grid_save_V_seq(&
                  vv                    , &
                  new_file_tmp          , &
                  ind_file_tmp          , &
                  maxsize_zo            , &
                  dimt_tmp              , &
                  close_file_tmp          &
                  )
          ENDIF
       ENDIF
    ENDIF

    !-- Compute transport --
    IF (key_partialsteps) THEN

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! partial steps case e3t is a 3D array
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * &
               MIN(e3t(:,:,kk,1), CSHIFT(e3t(:,:,kk,1),shift=1,dim=2))
        ENDDO
      ENDDO

    ELSE

      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ! e3t is a vector (no partial steps)
          vv(:,:,kk,ll) = vv(:,:,kk,ll) * e1v(:,:,1,1) * e3t(1,1,kk,1)
        ENDDO
      ENDDO

    ENDIF

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - Transport V: max ', MAXVAL(vv), ' min ', MINVAL(vv)
    ENDIF

    !=============== VERTICAL TRANSPORT =======================!
    IF (key_interp_temporal) THEN
      ww(:,:,:,2) = ww(:,:,:,1) ! We are in temporal interpolation case
    ENDIF

    IF (key_read_w) THEN

       new_file_tmp = new_file(3)
       
      !-- Read Vertical Current --!
      CALL sub_input_data_seq_main(                       &
           ncids(3), varids(3)                            , &
           new_file(3), ind_file(3)                       , &
           ind_time(3), ind_time_size(3)                  , &
           dimsorders(:,3)                                , &
           c_dir_ve,c_prefix_ve, maxsize_ve, c_suffix_ve  , &
           nc_var_ve,nc_var_eivw,nc_att_mask_ve,ww(:,:,:,1:1))

      ind_file_tmp   = ind_file(3)
      close_file_tmp = new_file(3)
      dimt_tmp       = ind_time_size(3)

      !-- Save W interpolated from B grid to C grid --!
      IF (key_B2C_save_data) THEN
         IF (PRESENT(first_call)) THEN
            IF (first_call) THEN
               CALL sub_B2C_grid_save_W_seq(&
                    ww                    , &
                    new_file_tmp          , &
                    ind_file_tmp          , &
                    maxsize_zo            , &
                    dimt_tmp              , &
                    close_file_tmp          &
                    )
               ENDIF
         ENDIF
      ENDIF

      !-- Compute transport --
      DO ll = 1, 1
        DO kk = 1, dims_reg(3,3)
          ww(:,:,kk,ll) = ww(:,:,kk,ll) * e1t(:,:,1,1) * e2t(:,:,1,1)
        ENDDO
      ENDDO

   ELSE

      !-- compute w transport from U and V transport --!
       WRITE(lun_standard,*)''
       WRITE(lun_standard,*)' - vertical current is computed -'
       CALL sub_w_comput(uu,vv,ww)

    ENDIF

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - Transport W: max ', MAXVAL(ww), ' min ', MINVAL(ww)
    ENDIF

    IF (key_alltracers) THEN

      !==================== TEMPERATURE =========================!
      IF (key_interp_temporal) THEN
        tt(:,:,:,2) = tt(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      !-- Read Temperature --!
      CALL sub_input_data_seq_main(                       &
           ncids(4), varids(4)                            , &
           new_file(4), ind_file(4)                       , &
           ind_time(4), ind_time_size(4)                  , &
           dimsorders(:,3)                                , &
           c_dir_te, c_prefix_te, maxsize_te, c_suffix_te , &
           nc_var_te,c_none,nc_att_mask_te,tt(:,:,:,1:1))

      !====================== SALINITY ===========================!
      IF (key_interp_temporal) THEN
        ss(:,:,:,2) = ss(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      !-- Read Salinity --!
      CALL sub_input_data_seq_main(                       &
           ncids(5), varids(5)                            , &
           new_file(5), ind_file(5)                       , &
           ind_time(5), ind_time_size(5)                  , &
           dimsorders(:,5)                                , &
           c_dir_sa,c_prefix_sa, maxsize_sa, c_suffix_sa  , &
           nc_var_sa,c_none,nc_att_mask_sa,ss(:,:,:,1:1))

      !======================= DENSITY ===========================!
      IF (key_interp_temporal) THEN
        rr(:,:,:,2) = rr(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      IF (key_computesigma) THEN

        !-- Compute Density from Temperature and Salinity --!
        DO l = 1, 1
          DO k = 1,  dims_reg(3,3)
            DO j = 1,  dims_reg(2,3)
              DO i = 1,  dims_reg(1,3)
                rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
              END DO
            END DO

          END DO
        ENDDO

      ELSE

        !-- Read Density --!
        CALL sub_input_data_seq_main(                     &
             ncids(6), varids(6)                            , &
             new_file(6), ind_file(6)                       , &
             ind_time(6), ind_time_size(6)                  , &
             dimsorders(:,6)                                , &
             c_dir_de,c_prefix_de,maxsize_de, c_suffix_de   , &
             nc_var_de,c_none,nc_att_mask_de,rr(:,:,:,1:1))

      ENDIF

    ENDIF

  END SUBROUTINE sub_input_data_B2C_gridz_seq
  !!***

  !=========================================================================
  !!****f* mod_input_data/sub_input_data_roms()
  !! NAME
  !!   sub_input_data_roms()
  !!
  !! FUNCTION
  !!   Data used in Ariane are on C grid and they are based on 
  !!   OPA conventions. ROMS data are on a C grid but don't respect
  !!   all OPA conventions.
  !!   The vertical axis is reversed.
  !!   U is defined only on (imt-1, jmt  , kmt-1)
  !!   V is defined only on (imt  , jmt-1, kmt-1)
  !!   W must be computed
  !!   T and S are defined on (imt, jmt, kmt-1)
  !!
  !!   * read ROMS tracer variables.
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute scale factor e3t.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (November 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_roms() 1,29

    INTEGER(kind=iprec) :: i,j,k,l
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    !!NG: 05/08/2008 REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
    !!NG: 05/08/2008      zeta, &
    !!NG: 05/08/2008      zw0 , &
    !!NG: 05/08/2008      e3u , &
    !!NG: 05/08/2008      e3v

    INTEGER(kind=iprec) :: ncid
    INTEGER(kind=iprec) :: varid

    LOGICAL :: key_ijt

    REAL(kind = rprec) :: hc ! S-coordinate parameter, critical depth (m)
    REAL(kind = rprec), DIMENSION(:), ALLOCATABLE :: &
         sc_w , & ! S-coordinate at W-points
         cs_w     ! S-coordinate stretching curves at W-points
    CHARACTER(len = 128)  :: c_filename   ! netcdf file name

    !======================= ZONAL CURRENT ==================!
    !- Read U on ROMS grid -!
    CALL sub_input_data_main(c_dir_zo,c_prefix_zo, &
         ind0_zo,indn_zo,maxsize_zo, &
         c_suffix_zo, nc_var_zo,nc_var_eivu,nc_att_mask_zo, &
         uu(1:dims_reg(1,3)-1,:,dims_reg(3,3)-1:1:-1,:))

    !============= MERIDIONAL CURRENT =======================!
    !- Read V on ROMS grid -!
    CALL sub_input_data_main(c_dir_me,c_prefix_me,&
         ind0_me,indn_me,maxsize_me, &
         c_suffix_me, nc_var_me,nc_var_eivv,nc_att_mask_me,&
         vv(:,1:dims_reg(2,3)-1,dims_reg(3,3)-1:1:-1,:))

    IF (key_alltracers) THEN
      !==================== TEMPERATURE =========================!
      !- Read T on ROMS grid -!
      CALL sub_input_data_main(c_dir_te,c_prefix_te,&
           ind0_te,indn_te,maxsize_te, &
           c_suffix_te, nc_var_te,c_none,nc_att_mask_te,&
           tt(:,:,dims_reg(3,3)-1:1:-1,:))

      !====================== SALINITY ===========================!
      !- Read S on ROMS grid -!
      CALL sub_input_data_main(c_dir_sa,c_prefix_sa, &
           ind0_sa,indn_sa,maxsize_sa, &
           c_suffix_sa, nc_var_sa,c_none,nc_att_mask_sa, &
           ss(:,:,dims_reg(3,3)-1:1:-1,:))

      !====================== COMPUTE RR =========================!
      WRITE(lun_standard,*)''
      WRITE(lun_standard,*)'------ Density is computing -------'
      !-- Compute Density from Temperature and Salinity --!
      DO l = 1, lmt
        DO k = 1,  dims_reg(3,3)
          DO j = 1,  dims_reg(2,3)
            DO i = 1,  dims_reg(1,3)
              rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
            END DO
          END DO
        END DO
      ENDDO
      !- Comments -!
      WRITE(lun_standard,*)' - Density: max ', MAXVAL(rr), ' min ', MINVAL(rr)

    ENDIF

    !====================== ZETA 2D + Time =======================!
    !- Read Zeta on ROMS grid -!
    ALLOCATE(zeta_a(dims_reg(1,3),dims_reg(2,3), 1, dims_reg(4,3)))
    call sub_memory(size(zeta_a),'r','zeta_a','sub_input_data_roms')

    key_ijt = .TRUE.

    CALL sub_input_data_main(c_dir_ze,c_prefix_ze, &
         ind0_ze,indn_ze,maxsize_ze, &
         c_suffix_ze, nc_var_ze,c_none,nc_att_mask_ze, &
         zeta_a(:,:,:,:), key_ijt = key_ijt)

    !========================= ZW0 ==============================!
    !- Read global attributes hc (scalar) sc_w and cs_w (vectors) -!
    !- in AVG files.
    !- Build file name from temperature file -!
    CALL sub_build_filename(ind0_te, maxsize_te, &
         c_prefix_te, c_suffix_te, c_filename)

    !======== Read ROMS global attributes hc, sc_w and Cs_w ======!
    ALLOCATE(sc_w(dims_reg(3,3)))
    call sub_memory(size(sc_w),'r','sc_w','sub_input_data_roms')

    ALLOCATE(cs_w(dims_reg(3,3)))
    call sub_memory(size(cs_w),'r','cs_w','sub_input_data_roms')

    !! NG : 29/09/2009
    IF (TRIM(which_roms) == 'agrif'.OR.TRIM(which_roms) == 'AGRIF') THEN 
       CALL sub_read_netcdf_global_att(                   &
            dir_glbatt, fn_glbatt                       , &
            nc_glbatt_hc, nc_glbatt_sc_w, nc_glbatt_Cs_w, &
            hc, sc_w, Cs_w)
    ELSE

       !- Open NetCDF file -!
       CALL sub_open_netcdf_file(dir_glbatt, fn_glbatt, ncid)

       !- Read hc -!
       CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_hc, varid)
       CALL sub_read_netcdf_var( ncid, varid, hc)

       !- Read sc_w -!
       CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_sc_w, varid)
       CALL sub_read_netcdf_var1d( ncid, varid, sc_w(dims_reg(3,3):1:-1))

       !- Read Cs_w -!
       CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_Cs_w, varid)
       CALL sub_read_netcdf_var1d( ncid, varid, Cs_w(dims_reg(3,3):1:-1))

       !- Close netcdf file -!
       CALL sub_close_netcdf_file(ncid)

    ENDIF

    !============== ZW0 (3D - constant in time) =================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZW0 is computing ------'
    WRITE(lun_standard,*)' - hc: '  , hc
    WRITE(lun_standard,*)' - SC_W: ', sc_w
    WRITE(lun_standard,*)' - CS_w: ', cs_w
    ALLOCATE(zw0(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
    call sub_memory(size(zw0),'r','zw0','sub_input_data_roms')

    DO k = 1, dims_reg(3,3)
      DO j = 1, dims_reg(2,3)
        DO i = 1, dims_reg(1,3)
          zw0(i,j,k,1) = hc * sc_w(k) + &
               (h_roms(i,j,1,1) - hc) * cs_w(k)
        ENDDO
      ENDDO
    ENDDO
    WRITE(lun_standard,*)' - ZW0 : max ', MAXVAL(zw0), ' min ', MINVAL(zw0)

    call sub_memory(-size(sc_w),'r','sc_w','sub_input_data_roms')
    DEALLOCATE(sc_w)

    call sub_memory(-size(cs_w),'r','cs_w','sub_input_data_roms')
    DEALLOCATE(cs_w)


    !==================== ZZ_WW (3D + time) ========================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZZ_WW is computing ------'

    ALLOCATE(zz_ww(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(zz_ww),'r','zz_ww','sub_input_data_roms')

    DO l = 1, dims_reg(4,3)
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3) 
          DO i = 1, dims_reg(1,3)
            zz_ww(i,j,k,l) = &
                 zeta_a(i,j,1,l) + zw0(i,j,k,1) * ( 1._rprec + zeta_a(i,j,1,l) / h_roms(i,j,1,1))
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    WRITE(lun_standard,*)' - ZZ_WW : max ', MAXVAL(zz_ww), ' min ', MINVAL(zz_ww)

    call sub_memory(-size(zw0),'r','zw0','sub_input_data_roms')
    DEALLOCATE(zw0)

    !=================== E3T (3D + time) ========================!
    !- Compute E3T from ZZ_WW -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ E3T scale factor is computing ------'
    ALLOCATE(e3t(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(e3t),'r','e3t','sub_input_data_roms')

    e3t(:,:,1:dims_reg(3,3)-1,:)= &
         zz_ww(:,:,1:dims_reg(3,3)-1,:) - zz_ww(:,:,2:dims_reg(3,3),:)

    e3t(:,:,dims_reg(3,3),:) = e3t(:,:,dims_reg(3,3)-1,:)
    WRITE(lun_standard,*)' - E3T : max ', MAXVAL(e3t), ' min ', MINVAL(e3t)

    !========================= E3U ==============================!
    ALLOCATE(e3u(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(e3u),'r','e3u','sub_input_data_roms')

    e3u(1:dims_reg(1,3)-1,:,:,:) = &
         .5_rprec *( e3t(1:dims_reg(1,3)-1,:,:,:) + e3t(2:dims_reg(1,3),:,:,:) )

    e3u(dims_reg(1,3),:,:,:) = e3u(dims_reg(1,3)-1,:,:,:)

    !========================= E3V ==============================!
    ALLOCATE(e3v(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(e3v),'r','e3v','sub_input_data_roms')

    e3v(:,1:dims_reg(2,3)-1,:,:) = &
         .5_rprec *( e3t(:,1:dims_reg(2,3)-1,:,:) + e3t(:,2:dims_reg(2,3),:,:) )

    e3v(:,dims_reg(2,3),:,:) = e3v(:,dims_reg(2,3)-1,:,:)

    !====================== TRANSPORT ===========================!
    !=== U ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Zonal transport is computing ------'
    DO l = 1, dims_reg(4,3)
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            uu(i,j,k,l) =  uu(i,j,k,l) * e3u(i,j,k,l) * e2u(i,j,1,1)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    WRITE(lun_standard,*)' - U Transport: max ', MAXVAL(uu), ' min ', MINVAL(uu)

    call sub_memory(-size(e3u),'r','e3u','sub_input_data_roms')
    DEALLOCATE(e3u)

    !=== V ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Meridional transport is computing ------'
    DO l = 1, dims_reg(4,3)
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            vv(i,j,k,l) =  vv(i,j,k,l) * e3v(i,j,k,l) * e1v(i,j,1,1)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    WRITE(lun_standard,*)' - V Transport: max ', MAXVAL(vv), ' min ', MINVAL(vv)

    call sub_memory(-size(e3v),'r','e3v','sub_input_data_roms')
    DEALLOCATE(e3v)


    !=== W ===!
    !- Must be no divergent -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Vertical transport is computing ------'
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    DO l = 1, dims_reg(4,3)
      DO k = dims_reg(3,3)-1, 1, -1
        DO j = 2, dims_reg(2,3)
          DO i = 2, dims_reg(1,3)
            ww(i,j,k,l) = ww(i,j,k+1,l) &
                 + uu(i-1,   j, k, l) - uu(i, j, k, l) &
                 + vv(  i, j-1, k, l) - vv(i, j, k, l)
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    DO l = 1, dims_reg(4,3)
      DO k = 2, dims_reg(3,3)-1
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            ww(i,j,k,l) = &
                 +   ww(i,j,k,l) &
                 - ( ww(i,j,1,l) / ( zz_ww(i,j,1,l) - zz_ww(i,j,dims_reg(3,3),l) ) ) &
                 *                 ( zz_ww(i,j,k,l) - zz_ww(i,j,dims_reg(3,3),l) )
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    ww(:,:, 1,:) = 0._rprec
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    WRITE(lun_standard,*)' - W Transport: max ', MAXVAL(ww), ' min ', MINVAL(ww)

  END SUBROUTINE sub_input_data_roms
  !!***
  !=========================================================================
  !!****f* mod_input_data_seq/sub_input_data_seq_roms()
  !! NAME
  !!   sub_input_data_seq_roms()
  !!
  !! FUNCTION
  !!   Data used in Ariane are on C grid and they are based on 
  !!   OPA conventions. ROMS data are on a C grid but don't respect
  !!   all OPA conventions.
  !!   The vertical axis is reversed.
  !!   U is defined only on (imt-1, jmt  , kmt-1)
  !!   V is defined only on (imt  , jmt-1, kmt-1)
  !!   W must be computed
  !!   T and S are defined on (imt, jmt, kmt-1)
  !!
  !!   * read ROMS tracer variables.
  !!   * loop on current and tracer variables: U, V, [W, T, S, R].
  !!   * call read netcdf.
  !!   * compute scale factor e3t.
  !!   * compute transport from currents.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (November 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_seq_roms(& 4,25
       ncids                       , &
       varids                      , &
       new_file                    , &
       ind_file                    , &
       ind_time                    , &
       ind_time_size               , &
       dimsorders                  )


    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ncids
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: varids
    LOGICAL               , DIMENSION(:), INTENT(inout) :: new_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time_size
    INTEGER(kind = iprec) , DIMENSION(:,:), INTENT(inout) :: dimsorders

    INTEGER(kind=iprec) :: i,j,k,l
    INTEGER(kind=iprec) :: ncid
    INTEGER(kind=iprec) :: varid
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    LOGICAL :: key_ijt

    !!NG: 05/08/2008 REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
    !!NG: 05/08/2008         zeta, &
    !!NG: 05/08/2008         zw0 , &
    !!NG: 05/08/2008         e3u , &
    !!NG: 05/08/2008         e3v

    REAL(kind = rprec) :: hc ! S-coordinate parameter, critical depth (m)
    REAL(kind = rprec), DIMENSION(:), ALLOCATABLE :: &
         sc_w , & ! S-coordinate at W-points
         cs_w     ! S-coordinate stretching curves at W-points
    !!NG: 05/08/2008 CHARACTER(len = 128)  :: c_filename   ! netcdf file name

    !======================= ZONAL CURRENT ==================!
    IF (key_interp_temporal) THEN
      uu(:,:,:,2) = uu(:,:,:,1) ! We are in temporal interpolation case
    ENDIF
    !- Read U on ROMS grid -!
    CALL sub_input_data_seq_main(                         &
         ncids(1), varids(1)                            , &
         new_file(1), ind_file(1)                       , &
         ind_time(1), ind_time_size(1)                  , &
         dimsorders(:,1)                                , &
         c_dir_zo,c_prefix_zo, maxsize_zo, c_suffix_zo  , &
         nc_var_zo,nc_var_eivu,nc_att_mask_zo           , &
         uu(1:dims_reg(1,3)-1,:,dims_reg(3,3)-1:1:-1,1:1))

    !============= MERIDIONAL CURRENT =======================!
    IF (key_interp_temporal) THEN
      vv(:,:,:,2) = vv(:,:,:,1) ! We are in temporal interpolation case
    ENDIF
    !- Read V on ROMS grid -!
    CALL sub_input_data_seq_main(                         &
         ncids(2), varids(2)                            , &
         new_file(2), ind_file(2)                       , &
         ind_time(2), ind_time_size(2)                  , &
         dimsorders(:,2)                                , &
         c_dir_me,c_prefix_me,maxsize_me, c_suffix_me   , &
         nc_var_me,nc_var_eivv,nc_att_mask_me           , &
         vv(:,1:dims_reg(2,3)-1,dims_reg(3,3)-1:1:-1,1:1))

    IF (key_alltracers) THEN
      !==================== TEMPERATURE =========================!
      IF (key_interp_temporal) THEN
        tt(:,:,:,2) = tt(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      !- Read T on ROMS grid -!
      CALL sub_input_data_seq_main(                       &
           ncids(4), varids(4)                            , &
           new_file(4), ind_file(4)                       , &
           ind_time(4), ind_time_size(4)                  , &
           dimsorders(:,4)                                , &
           c_dir_te,c_prefix_te,maxsize_te, c_suffix_te   , &
           nc_var_te,c_none,nc_att_mask_te                , &
           tt(:,:,dims_reg(3,3)-1:1:-1,1:1))

      !====================== SALINITY ===========================!
      IF (key_interp_temporal) THEN
        ss(:,:,:,2) = ss(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      !- Read S on ROMS grid -!
      CALL sub_input_data_seq_main(                       &
           ncids(5), varids(5)                            , &
           new_file(5), ind_file(5)                       , &
           ind_time(5), ind_time_size(5)                  , &
           dimsorders(:,5)                                , &
           c_dir_sa,c_prefix_sa, maxsize_sa, c_suffix_sa  , &
           nc_var_sa,c_none,nc_att_mask_sa                , &
           ss(:,:,dims_reg(3,3)-1:1:-1,1:1))

      !====================== COMPUTE RR =========================!
      IF (key_interp_temporal) THEN
        rr(:,:,:,2) = rr(:,:,:,1) ! We are in temporal interpolation case
      ENDIF
      WRITE(lun_standard,*)''
      WRITE(lun_standard,*)'------ Density is computing -------'
      !-- Compute Density from Temperature and Salinity --!
      DO l = 1, 1
        DO k = 1,  dims_reg(3,3)
          DO j = 1,  dims_reg(2,3)
            DO i = 1,  dims_reg(1,3)
              rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
            END DO
          END DO
        END DO
      ENDDO
      !- Comments -!
      IF (id_comments) THEN
        WRITE(lun_standard,*)' - Density: max ', MAXVAL(rr), ' min ', MINVAL(rr)
      ENDIF
    ENDIF

    !====================== ZETA 2D + Time =======================!
    !- Read Zeta on ROMS grid -!
    IF (.NOT.ALLOCATED(zeta_a)) THEN
       ALLOCATE(zeta_a(dims_reg(1,3),dims_reg(2,3), 1, 1))
       call sub_memory(size(zeta_a),'r','zeta_a','sub_input_data_seq_roms')
    ENDIF

    key_ijt = .TRUE.

    CALL sub_input_data_seq_main(                         &
         ncids(7), varids(7)                            , &
         new_file(7), ind_file(7)                       , &
         ind_time(7), ind_time_size(7)                  , &
         dimsorders(:,7)                                , &
         c_dir_ze,c_prefix_ze, maxsize_ze, c_suffix_ze  , &
         nc_var_ze,c_none,nc_att_mask_ze, &
         zeta_a(:,:,:,:), key_ijt = key_ijt)

    !! NG: 05/08/2008: We assume that ZW0 is constant in time
    !! NG: 05/08/2008: It will be read only one time...

    IF (.NOT.ALLOCATED(zw0)) THEN

       !! NG: 05/08/2008!========================= ZW0 ==============================!
       !! NG: 05/08/2008!- Read global attributes hc (scalar) sc_w and cs_w (vectors) -!
       !! NG: 05/08/2008!- in AVG files.
       !! NG: 05/08/2008!- Build file name from temperature file -!
       !! NG: 05/08/2008 CALL sub_build_filename(ind0_te, maxsize_te, &
       !! NG: 05/08/2008     c_prefix_te, c_suffix_te, c_filename)

       !======== Read ROMS global attributes hc, sc_w and Cs_w ======!
       IF (.NOT.ALLOCATED(sc_w)) THEN
          ALLOCATE(sc_w(dims_reg(3,3)))
          !!NG 1dec2009 write(*,*)' sc_w: ',size(sc_w)
          call sub_memory(size(sc_w),'r','sc_w','sub_input_data_seq_roms')
       ENDIF

       IF (.NOT.ALLOCATED(cs_w)) THEN
          ALLOCATE(cs_w(dims_reg(3,3)))
          !!NG 1dec2009 write(*,*)' cs_w: ',size(cs_w)
          call sub_memory(size(cs_w),'r','cs_w','sub_input_data_seq_roms')
       ENDIF

       !! NG : 29/09/2009
       IF (TRIM(which_roms) == 'agrif'.OR.TRIM(which_roms) == 'AGRIF') THEN 

          CALL sub_read_netcdf_global_att(                   &
               dir_glbatt, fn_glbatt                       , &
               nc_glbatt_hc, nc_glbatt_sc_w, nc_glbatt_Cs_w, &
               hc, sc_w, Cs_w)

       ELSE

          !- Open NetCDF file -!
          CALL sub_open_netcdf_file(dir_glbatt, fn_glbatt, ncid)

          !- Read hc -!
          CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_hc, varid)
          CALL sub_read_netcdf_var( ncid, varid, hc)

          !- Read sc_w -!
          CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_sc_w, varid)
          CALL sub_read_netcdf_var1d( ncid, varid, sc_w(dims_reg(3,3):1:-1))

          !- Read Cs_w -!
          CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_Cs_w, varid)
          CALL sub_read_netcdf_var1d( ncid, varid, Cs_w(dims_reg(3,3):1:-1))

          !- Close netcdf file -!
          CALL sub_close_netcdf_file(ncid)

       ENDIF

       !============== ZW0 (3D - constant in time) =================!
       WRITE(lun_standard,*)''
       WRITE(lun_standard,*)' ------ ZW0 is computing ------'
       WRITE(lun_standard,*)' - hc: ',hc
       WRITE(lun_standard,*)' - SC_W: ',sc_w
       WRITE(lun_standard,*)' - CS_w: ',cs_w

       IF (.NOT.ALLOCATED(zw0)) THEN
          ALLOCATE(zw0(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
          call sub_memory(size(zw0),'r','zw0','sub_input_data_seq_roms')
       ENDIF

       DO k = 1, dims_reg(3,3)
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                zw0(i,j,k,1) = hc * sc_w(k) + &
                     (h_roms(i,j,1,1) - hc) * cs_w(k)
             ENDDO
          ENDDO
       ENDDO
       !!    DEALLOCATE(sc_w)
       !!    DEALLOCATE(cs_w)
       !! NG:10/07/1009
       call sub_memory(-size(sc_w),'r','sc_w','sub_input_approx_zz_ww_roms')
       DEALLOCATE(sc_w)
       call sub_memory(-size(cs_w),'r','cs_w','sub_input_approx_zz_ww_roms')
       DEALLOCATE(cs_w)
       !! NG:10/07/1009

    ELSE

       WRITE(lun_standard,*)''
       WRITE(lun_standard,*) &
            '  ------  We assume that ZW0 is constant in time. ------'
       
    ENDIF

    IF (id_comments) THEN
       WRITE(lun_standard,*)' - ZW0 : max ', MAXVAL(zw0), ' min ', MINVAL(zw0)
    ENDIF

    !==================== ZZ_WW (3D + time) ========================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZZ_WW is computing ------'

    IF (.NOT.ALLOCATED(zz_ww)) THEN
      ALLOCATE(zz_ww(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3),1))
      call sub_memory(size(zz_ww),'r','zz_ww','sub_input_data_seq_roms')
    ENDIF

    DO l = 1, 1
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            zz_ww(i,j,k,l) = &
                 zeta_a(i,j,1,l) + zw0(i,j,k,1) * ( 1. + zeta_a(i,j,1,l) / h_roms(i,j,1,1))
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - ZZ_WW : max ', MAXVAL(zz_ww), ' min ', MINVAL(zz_ww)
    ENDIF
    !!    DEALLOCATE(zw0)

    !=================== E3T (3D + time) ========================!
    !- Compute E3T from ZZ_WW -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ E3T scale factor is computing ------'

    IF (.NOT.ALLOCATED(e3t)) THEN
      ALLOCATE(e3t(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
      call sub_memory(size(e3t),'r','e3t','sub_input_data_seq_roms')
    ENDIF

    e3t(:,:,1:dims_reg(3,3)-1,:)= &
         zz_ww(:,:,1:dims_reg(3,3)-1,:) - zz_ww(:,:,2:dims_reg(3,3),:)

    e3t(:,:,dims_reg(3,3),:) = e3t(:,:,dims_reg(3,3)-1,:)

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - E3T : max ', MAXVAL(e3t), ' min ', MINVAL(e3t)
    ENDIF
    !========================= E3U ==============================!
    IF (.NOT.ALLOCATED(e3u)) THEN
      ALLOCATE(e3u(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
      call sub_memory(size(e3u),'r','e3u','sub_input_data_seq_roms')
    ENDIF

    e3u(1:dims_reg(1,3)-1,:,:,:) = &
         .5_rprec *( e3t(1:dims_reg(1,3)-1,:,:,:) + e3t(2:dims_reg(1,3),:,:,:) )

    e3u(dims_reg(1,3),:,:,:) = e3u(dims_reg(1,3)-1,:,:,:)

    !========================= E3V ==============================!
    IF (.NOT.ALLOCATED(e3v)) THEN
      ALLOCATE(e3v(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
      call sub_memory(size(e3v),'r','e3u','sub_input_data_seq_roms')
    ENDIF

    e3v(:,1:dims_reg(2,3)-1,:,:) = &
         .5_rprec *( e3t(:,1:dims_reg(2,3)-1,:,:) + e3t(:,2:dims_reg(2,3),:,:) )

    e3v(:,dims_reg(2,3),:,:) = e3v(:,dims_reg(2,3)-1,:,:)

    !====================== TRANSPORT ===========================!
    !=== U ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Zonal transport is computing ------'
    DO l = 1, 1
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            uu(i,j,k,l) =  uu(i,j,k,l) * e3u(i,j,k,l) * e2u(i,j,1,1)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    IF (id_comments) THEN
      WRITE(lun_standard,*)' - U Transport: max ', MAXVAL(uu), ' min ', MINVAL(uu)
    ENDIF
    !!    DEALLOCATE(e3u)

    !=== V ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Meridional transport is computing ------'
    DO l = 1, 1
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            vv(i,j,k,l) =  vv(i,j,k,l) * e3v(i,j,k,l) * e1v(i,j,1,1)
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - V Transport: max ', MAXVAL(vv), ' min ', MINVAL(vv)
    ENDIF
    !!    DEALLOCATE(e3v)

    !=== W ===!
    IF (key_interp_temporal) THEN
      ww(:,:,:,2) = ww(:,:,:,1) ! We are in temporal interpolation case
    ENDIF
    !- Must be no divergent -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Vertical transport is computing ------'
    ww(:,:,dims_reg(3,3),1:1) = 0._rprec

    DO l = 1, 1
      DO k = dims_reg(3,3)-1, 1, -1
        DO j = 2, dims_reg(2,3)
          DO i = 2, dims_reg(1,3)
            ww(i,j,k,l) = ww(i,j,k+1,l) &
                 + uu(i-1,   j, k, l) - uu(i, j, k, l) &
                 + vv(  i, j-1, k, l) - vv(i, j, k, l)
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    DO l = 1, 1
      DO k = 2, dims_reg(3,3)-1
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            ww(i,j,k,l) = &
                 +   ww(i,j,k,l) &
                 - ( ww(i,j,1,l) / ( zz_ww(i,j,1,l) - zz_ww(i,j,dims_reg(3,3),l) ) ) &
                 *                 ( zz_ww(i,j,k,l) - zz_ww(i,j,dims_reg(3,3),l) )
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    ww(:,:, 1,:) = 0._rprec
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - W Transport: max ', MAXVAL(ww), ' min ', MINVAL(ww)
    ENDIF

  END SUBROUTINE sub_input_data_seq_roms
 !!***
  !=========================================================================
  !!****f* mod_input_data/sub_input_approx_zz_ww_roms()
  !! NAME
  !!  sub_input_approx_zz_ww_roms ()
  !!
  !! FUNCTION
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (February 2008)
  !! 
  !! CREATION DATE
  !!   * February 2008
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec_seq
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_approx_zz_ww_roms() 2,17

    !! WE ASSUME HERE THAT ZETA=0. !!

    REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: zw0

    REAL(kind = rprec) :: hc ! S-coordinate parameter, critical depth (m)
    REAL(kind = rprec), DIMENSION(:), ALLOCATABLE :: &
         sc_w , & ! S-coordinate at W-points
         cs_w     ! S-coordinate stretching curves at W-points
    CHARACTER(len = 128)  :: c_filename   ! netcdf file name

    INTEGER(kind=iprec) :: i,j,k,l
    INTEGER(kind=iprec) :: ncid
    INTEGER(kind=iprec) :: varid

    !========================= ZW0 ==============================!
    !- Read global attributes hc (scalar) sc_w and cs_w (vectors) -!
    !- in AVG files.
    !- Build file name from temperature file -!
    CALL sub_build_filename(ind0_te, maxsize_te, &
         c_prefix_te, c_suffix_te, c_filename)

    !======== Read ROMS global attributes hc, sc_w and Cs_w ======!
    IF (.NOT.ALLOCATED(sc_w)) THEN
      ALLOCATE(sc_w(dims_reg(3,3)))
      call sub_memory(size(sc_w),'r','sc_w','sub_input_approx_zz_ww_roms')
    ENDIF

    IF (.NOT.ALLOCATED(cs_w)) THEN
      ALLOCATE(cs_w(dims_reg(3,3)))
      call sub_memory(size(cs_w),'r','cs_w','sub_input_approx_zz_ww_roms')
    ENDIF

    !! NG : 29/09/2009
    IF (TRIM(which_roms) == 'agrif'.OR.TRIM(which_roms) == 'AGRIF') THEN 
       CALL sub_read_netcdf_global_att(                   &
            dir_glbatt, fn_glbatt                       , &
            nc_glbatt_hc, nc_glbatt_sc_w, nc_glbatt_Cs_w, &
            hc, sc_w, Cs_w)
    ELSE

       !- Open NetCDF file -!
       CALL sub_open_netcdf_file(dir_glbatt, fn_glbatt, ncid)

       !- Read hc -!
       CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_hc, varid)
       CALL sub_read_netcdf_var( ncid, varid, hc)

       !- Read sc_w -!
       CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_sc_w, varid)
       CALL sub_read_netcdf_var1d( ncid, varid, sc_w)

       !- Read Cs_w -!
       CALL sub_read_netcdf_varid_ndims(ncid, nc_glbatt_Cs_w, varid)
       CALL sub_read_netcdf_var1d( ncid, varid, Cs_w)

       !- Close netcdf file -!
       CALL sub_close_netcdf_file(ncid)

    ENDIF

    !============== ZW0 (3D - constant in time) =================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZW0 is computing ------'
    WRITE(lun_standard,*)' - hc: ',hc
    WRITE(lun_standard,*)' - SC_W: ',sc_w
    WRITE(lun_standard,*)' - CS_w: ',cs_w

    IF (.NOT.ALLOCATED(zw0)) THEN
      ALLOCATE(zw0(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
      call sub_memory(size(zw0),'r','zw0','sub_input_approx_zz_ww_roms')
    ENDIF

    DO k = 1, dims_reg(3,3)
      DO j = 1, dims_reg(2,3)
        DO i = 1, dims_reg(1,3)
          zw0(i,j,k,1) = hc * sc_w(k) + &
               (h_roms(i,j,1,1) - hc) * cs_w(k)
        ENDDO
      ENDDO
    ENDDO

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - ZW0 : max ', MAXVAL(zw0), ' min ', MINVAL(zw0)
    ENDIF

    !==================== ZZ_WW (3D + time) ========================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZZ_WW is computing ------'

    IF (.NOT.ALLOCATED(zz_ww)) THEN
      ALLOCATE(zz_ww(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3),1))
      call sub_memory(size(zz_ww),'r','zz_ww','sub_input_approx_zz_ww_roms')
    ENDIF

    DO l = 1, 1
      DO k = 1, dims_reg(3,3)
        DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
            zz_ww(i,j,k,l) = zw0(i,j,k,1)
          ENDDO
        ENDDO
      ENDDO
    ENDDO

    IF (id_comments) THEN
      WRITE(lun_standard,*)' - ZZ_WW : max ', MAXVAL(zz_ww), ' min ', MINVAL(zz_ww)
    ENDIF

    call sub_memory(-size(zw0),'r','zw0','sub_input_approx_zz_ww_roms')
    DEALLOCATE(zw0)
    call sub_memory(-size(sc_w),'r','sc_w','sub_input_approx_zz_ww_roms')
    DEALLOCATE(sc_w)
    call sub_memory(-size(cs_w),'r','cs_w','sub_input_approx_zz_ww_roms')
    DEALLOCATE(cs_w)

    
  END SUBROUTINE sub_input_approx_zz_ww_roms
  !!***
  !=========================================================================
  !!****f* mod_input_data/sub_input_data_symphonie()
  !! NAME
  !!   sub_input_data_symphonie()
  !!
  !! FUNCTION
  !!   Data used in Ariane are on C grid and they are based on 
  !!   OPA conventions. SYMPHONIE data are on a C grid but don't respect
  !!   all OPA conventions.
  !!   The vertical axis is reversed.
  !!   U is defined only on (imt-1, jmt  , kmt-1)
  !!   V is defined only on (imt  , jmt-1, kmt-1)
  !!   W must be computed
  !!   T and S are defined on (imt, jmt, kmt-1)
  !!
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April 2007)
  !! 
  !! CREATION DATE
  !!   * April 2007
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_symphonie() 1,19

    INTEGER(kind=iprec) :: i,j,k,l
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    INTEGER(kind = iprec) :: kmax

    REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
         sse  , &
         zw0, &
         e3u  , &
         e3v

    REAL(kind=rprec), DIMENSION(:,:), ALLOCATABLE :: depth_max

    kmax =  dims_reg(3,3)


    !======================= ZONAL CURRENT ==================!
    !- Read U on SYMPHONIE grid -!
    CALL sub_input_data_main(c_dir_zo,c_prefix_zo, &
         ind0_zo,indn_zo,maxsize_zo, &
         c_suffix_zo, nc_var_zo,nc_var_eivu,nc_att_mask_zo, &
         uu(:,:,dims_reg(3,3)-1:1:-1,:))

    uu(:,:,:,:) = CSHIFT(uu(:,:,:,:), shift=1, dim=1)
    uu(dims_reg(1,3),:,:,:) = 0._rprec

    !============= MERIDIONAL CURRENT =======================!
    !- Read V on SYMPHONIE grid -!
    CALL sub_input_data_main(c_dir_me,c_prefix_me,&
         ind0_me,indn_me,maxsize_me, &
         c_suffix_me, nc_var_me,nc_var_eivv,nc_att_mask_me,&
         vv(:,:,dims_reg(3,3)-1:1:-1,:))

    vv(:,:,:,:) = CSHIFT(vv(:,:,:,:), shift=1, dim=2)
    vv(:,dims_reg(2,3),:,:) = 0._rprec

    IF (key_alltracers) THEN
       !==================== TEMPERATURE =========================!
       !- Read T on SYMPHONIE grid -!
       CALL sub_input_data_main(c_dir_te,c_prefix_te,&
            ind0_te,indn_te,maxsize_te, &
            c_suffix_te, nc_var_te,c_none,nc_att_mask_te,&
            tt(:,:,dims_reg(3,3)-1:1:-1,:))

       !====================== SALINITY ===========================!
       !- Read S on SYMPHONIE grid -!
       CALL sub_input_data_main(c_dir_sa,c_prefix_sa, &
            ind0_sa,indn_sa,maxsize_sa, &
            c_suffix_sa, nc_var_sa,c_none,nc_att_mask_sa, &
            ss(:,:,dims_reg(3,3)-1:1:-1,:))



       IF (key_computesigma) THEN

          !====================== COMPUTE RR =========================!
          WRITE(lun_standard,*)''
          WRITE(lun_standard,*)'------ Density is computing -------'
          !-- Compute Density from Temperature and Salinity --!
          DO l = 1, lmt
             DO k = 1,  dims_reg(3,3)-1
                DO j = 1,  dims_reg(2,3)
                   DO i = 1,  dims_reg(1,3)
                      rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
                   END DO
                END DO
             END DO
          ENDDO

       ELSE

          !-- Read Density --!
          CALL sub_input_data_main(                           &
               c_dir_de,c_prefix_de                         , &
               ind0_de,indn_de, maxsize_de                  , &
               c_suffix_de, nc_var_de,c_none,nc_att_mask_de , &
               rr(:,:,dims_reg(3,3)-1:1:-1,1:1) )

       ENDIF
    ENDIF

    !- Comments -!
    WRITE(lun_standard,*)' - Density: max ', MAXVAL(rr), ' min ', MINVAL(rr)

    !====================== READ ZZ_TT =======================!
    !- Read ZZ_TT on SYMPHONIE grid -!
    !- surface is at indice kmax and deep at 0.

    ALLOCATE(zw0(dims_reg(1,3),dims_reg(2,3),dims_reg(3,3), 1))
    call sub_memory(size(zw0),'r','zw0','sub_input_data_symphonie')

    zw0(:,:,:,:) = zz_ww(:,:,:,:)

    !====================== READ SSE 2D + Time =======================!
    !- Read SSE on SYMPHONIE grid -!
    ALLOCATE(sse(dims_reg(1,3),dims_reg(2,3), 1, dims_reg(4,3)))
    call sub_memory(size(sse),'r','sse','sub_input_data_symphonie')

    CALL sub_input_data_main(c_dir_sse,c_prefix_sse, &
         ind0_sse,indn_sse,maxsize_sse, &
         c_suffix_sse, nc_var_sse,c_none,nc_att_mask_sse, &
         sse(:,:,:,:))

    !==================== COMPUTE ZZ_WW (3D + time) =====================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZZ_WW is computing ------'
    ALLOCATE(zz_ww(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(zz_ww),'r','zz_ww','sub_input_data_symphonie')

    zz_ww(:,:,:,:)=0._rprec

    ALLOCATE(depth_max(dims_reg(1,3),dims_reg(2,3)))
    call sub_memory(size(depth_max),'r','depth_max','sub_input_data_symphonie')

    depth_max(:,:) = MAXVAL(ABS(zw0(:,:,1:kmax,1)), dim = 3)

    DO l = 1, dims_reg(4,3)
       DO k = 1, kmax
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                IF (depth_max(i,j) /= 0._RPREC) THEN
                   zz_ww(i,j,k,l) = sse(i,j,1,l) + zw0(i,j,k,1) * &
                        (1 + ( sse(i,j,1,l)/depth_max(i,j) ) )
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    call sub_memory(-size(depth_max),'r','depth_max','sub_input_data_symphonie')
    DEALLOCATE(depth_max)

    WRITE(lun_standard,*)' - ZZ_WW : max ', MAXVAL(zz_ww), ' min ', MINVAL(zz_ww)

    !=================== E3T (3D + time) ========================!
    !- Compute E3T from ZZ_WW -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ E3T scale factor is computing ------'
    ALLOCATE(e3t(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(e3t),'r','e3t','sub_input_data_symphonie')

    e3t(:,:,1:dims_reg(3,3)-1,:)= &
         zz_ww(:,:,1:dims_reg(3,3)-1,:) - zz_ww(:,:,2:dims_reg(3,3),:)

    e3t(:,:,dims_reg(3,3),:) = e3t(:,:,dims_reg(3,3)-1,:)

    e3t(:,:,:,:) = e3t(:,:,:,:) * tmask(:,:,:,:)

    WRITE(lun_standard,*)' - E3T : max ', MAXVAL(e3t), ' min ', MINVAL(e3t)

    !========================= E3U ==============================!
    ALLOCATE(e3u(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(e3u),'r','e3u','sub_input_data_symphonie')
    
    e3u(1:dims_reg(1,3)-1,:,:,:) = &
         .5_rprec *( e3t(1:dims_reg(1,3)-1,:,:,:) + e3t(2:dims_reg(1,3),:,:,:) )

    e3u(dims_reg(1,3),:,:,:) = e3u(dims_reg(1,3)-1,:,:,:)

    !========================= E3V ==============================!
    ALLOCATE(e3v(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), dims_reg(4,3)))
    call sub_memory(size(e3u),'r','e3u','sub_input_data_symphonie')

    e3v(:,1:dims_reg(2,3)-1,:,:) = &
         .5_rprec *( e3t(:,1:dims_reg(2,3)-1,:,:) + e3t(:,2:dims_reg(2,3),:,:) )

    e3v(:,dims_reg(2,3),:,:) = e3v(:,dims_reg(2,3)-1,:,:)

    !====================== TRANSPORT ===========================!
    !=== U ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Zonal transport is computing ------'
    DO l = 1, dims_reg(4,3)
       DO k = 1, dims_reg(3,3)
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                uu(i,j,k,l) =  uu(i,j,k,l) * e3u(i,j,k,l) * e2u(i,j,1,1)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    WRITE(lun_standard,*)' - U Transport: max ', MAXVAL(uu), ' min ', MINVAL(uu)

    call sub_memory(-size(e3u),'r','e3u','sub_input_data_symphonie')
    DEALLOCATE(e3u)
    

    !=== V ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Meridional transport is computing ------'
    DO l = 1, dims_reg(4,3)
       DO k = 1, dims_reg(3,3)
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                vv(i,j,k,l) =  vv(i,j,k,l) * e3v(i,j,k,l) * e1v(i,j,1,1)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    WRITE(lun_standard,*)' - V Transport: max ', MAXVAL(vv), ' min ', MINVAL(vv)

    call sub_memory(-size(e3v),'r','e3v','sub_input_data_symphonie')
    DEALLOCATE(e3v)


    !=== W ===!
    !- Must be no divergent -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Vertical transport is computing ------'
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    DO l = 1, dims_reg(4,3)
       DO k = dims_reg(3,3)-1, 1, -1
          DO j = 2, dims_reg(2,3)-1
             DO i = 2, dims_reg(1,3)-1
                ww(i,j,k,l) = ww(i,j,k+1,l) &
                     + uu(i-1,   j, k, l) - uu(i, j, k, l) &
                     + vv(  i, j-1, k, l) - vv(i, j, k, l)
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    DO l = 1, dims_reg(4,3)
       DO k = 2, dims_reg(3,3)-1
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                IF (zz_ww(i,j,1,l) - zz_ww(i,j,dims_reg(3,3),l) /= 0._rprec) THEN
                   ww(i,j,k,l) = &
                        +   ww(i,j,k,l) &
                        - ( ww(i,j,1,l) / ( zz_ww(i,j,1,l) - zz_ww(i,j,dims_reg(3,3),l) ) ) &
                        *                 ( zz_ww(i,j,k,l) - zz_ww(i,j,dims_reg(3,3),l) )
                ELSE
                   ww(i,j,k,l) = 0._rprec
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    ww(:,:, 1,:) = 0._rprec
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    WRITE(lun_standard,*)' - W Transport: max ', MAXVAL(ww), ' min ', MINVAL(ww)

    zz_ww(:,:,:,:) = zw0(:,:,:,:)

    call sub_memory(-size(zw0),'r','zw0','sub_input_data_symphonie')
    DEALLOCATE(zw0)

    call sub_memory(-size(sse),'r','sse','sub_input_data_symphonie')
    DEALLOCATE(sse)


  END SUBROUTINE sub_input_data_symphonie
  !!***
  !=========================================================================
  !!****f* mod_input_data/sub_input_data_seq_symphonie()
  !! NAME
  !!   sub_input_data_seq_symphonie()
  !!
  !! FUNCTION
  !!   Data used in Ariane are on C grid and they are based on 
  !!   OPA conventions. SYMPHONIE data are on a C grid but don't respect
  !!   all OPA conventions.
  !!   The vertical axis is reversed.
  !!   U is defined only on (imt-1, jmt  , kmt-1)
  !!   V is defined only on (imt  , jmt-1, kmt-1)
  !!   W must be computed
  !!   T and S are defined on (imt, jmt, kmt-1)
  !!
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April 2007)
  !! 
  !! CREATION DATE
  !!   * April 2007
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * no arguments
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_seq_symphonie( & 4,21
       ncids                             , &
       varids                            , &
       new_file                          , &
       ind_file                          , &
       ind_time                          , &
       ind_time_size                     , &
       dimsorders                        )

    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ncids
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: varids
    LOGICAL               , DIMENSION(:), INTENT(inout) :: new_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_file
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time
    INTEGER(kind = iprec) , DIMENSION(:), INTENT(inout) :: ind_time_size
    INTEGER(kind = iprec) , DIMENSION(:,:), INTENT(inout) :: dimsorders

    INTEGER(kind=iprec) :: i,j,k,l
    CHARACTER(len=4), PARAMETER :: c_none='NONE'

    INTEGER(kind = iprec) :: kmax

    REAL(kind=rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: &
         sse  , &
         zw0, &
         e3u  , &
         e3v

    REAL(kind=rprec), DIMENSION(:,:), ALLOCATABLE :: depth_max
    INTEGER(kind = rprec) , DIMENSION(:,:), ALLOCATABLE :: ind_depth_max
    REAL(kind=rprec):: miss_val


    kmax =  dims_reg(3,3)

    !======================= ZONAL CURRENT ==================!
    !- Read U on SYMPHONIE grid -!
    CALL sub_input_data_seq_main(                         &
         ncids(1), varids(1)                            , &
         new_file(1), ind_file(1)                       , &
         ind_time(1), ind_time_size(1)                  , &
         dimsorders(:,1)                                , &
         c_dir_zo,c_prefix_zo, maxsize_zo, c_suffix_zo  , &
         nc_var_zo,nc_var_eivu,nc_att_mask_zo           , &
         uu(:,:,dims_reg(3,3)-1:1:-1,1:1))

    uu(:,:,:,:) = CSHIFT(uu(:,:,:,:), shift=1, dim=1)
    uu(dims_reg(1,3),:,:,:) = 0._rprec

    !============= MERIDIONAL CURRENT =======================!
    !- Read V on SYMPHONIE grid -!
    CALL sub_input_data_seq_main(                         &
         ncids(2), varids(2)                            , &
         new_file(2), ind_file(2)                       , &
         ind_time(2), ind_time_size(2)                  , &
         dimsorders(:,2)                                , &
         c_dir_me,c_prefix_me,maxsize_me, c_suffix_me   , &
         nc_var_me,nc_var_eivv,nc_att_mask_me           , &
         vv(:,:,dims_reg(3,3)-1:1:-1,1:1))

    vv(:,:,:,:) = CSHIFT(vv(:,:,:,:), shift=1, dim=2)
    vv(:,dims_reg(2,3),:,:) = 0._rprec

    IF (key_alltracers) THEN
       !==================== TEMPERATURE =========================!
       !- Read T on SYMPHONIE grid -!
       CALL sub_input_data_seq_main(                       &
            ncids(4), varids(4)                          , &
            new_file(4), ind_file(4)                     , &
            ind_time(4), ind_time_size(4)                , &
            dimsorders(:,4)                              , &
            c_dir_te,c_prefix_te,maxsize_te, c_suffix_te , &
            nc_var_te,c_none,nc_att_mask_te              , &
            tt(:,:,dims_reg(3,3)-1:1:-1,1:1))

       !- Read mask value in netcdf file -!
       CALL sub_read_netcdf_att_val(ncids(4), varids(4), &
            TRIM(nc_att_mask_te), miss_val,0._rprec)

       !====================== SALINITY ===========================!
       !- Read S on SYMPHONIE grid -!
       CALL sub_input_data_seq_main(                       &
            ncids(5), varids(5)                          , &
            new_file(5), ind_file(5)                     , &
            ind_time(5), ind_time_size(5)                , &
            dimsorders(:,5)                              , &
            c_dir_sa,c_prefix_sa, maxsize_sa, c_suffix_sa, &
            nc_var_sa,c_none,nc_att_mask_sa              , &
            ss(:,:,dims_reg(3,3)-1:1:-1,1:1))


       IF (key_computesigma) THEN
          !====================== COMPUTE RR =========================!
          WRITE(lun_standard,*)''
          WRITE(lun_standard,*)'------ Density is computing -------'
          !-- Compute Density from Temperature and Salinity --!
          DO l = 1, 1
             DO k = 1,  dims_reg(3,3)-1
                DO j = 1,  dims_reg(2,3)
                   DO i = 1,  dims_reg(1,3)
                      rr(i,j,k,l)=rhostp(zsigma,ss(i,j,k,l),tt(i,j,k,l),-ABS(zsigma))
                   END DO
                END DO
             END DO
          ENDDO

       ELSE

          !-- Read Density --!
          CALL sub_input_data_seq_main(                         &
               ncids(6), varids(6)                            , &
               new_file(6), ind_file(6)                       , &
               ind_time(6), ind_time_size(6)                  , &
               dimsorders(:,6)                                , &
               c_dir_de,c_prefix_de,maxsize_de, c_suffix_de   , &
               nc_var_de,c_none,nc_att_mask_de                , &
               rr(:,:,dims_reg(3,3)-1:1:-1,1:1))

       ENDIF

       !- Comments -!
       IF (id_comments) THEN
          WRITE(lun_standard,*)' - Density: max ', MAXVAL(rr), ' min ', MINVAL(rr)
       ENDIF
    ENDIF

    !====================== READ ZZ_WW =======================!
    !- Read ZZ_WW on SYMPHONIE grid -!
    !- surface is at indice kmax and deep at 0.

    ALLOCATE(zw0(dims_reg(1,3),dims_reg(2,3),dims_reg(3,3), 1))
    call sub_memory(size(zw0),'r','zw0','sub_input_data_seq_symphonie')

    zw0(:,:,:,:) = zz_ww(:,:,:,:)

    !====================== READ SSE 2D + Time =======================!
    !- Read SSE on SYMPHONIE grid -!
    ALLOCATE(sse(dims_reg(1,3),dims_reg(2,3), 1, 1))
    call sub_memory(size(sse),'r','sse','sub_input_data_seq_symphonie')

    CALL sub_input_data_seq_main(                               &
         ncids(7), varids(7)                                  , &
         new_file(7), ind_file(7)                             , &
         ind_time(7), ind_time_size(7)                        , &
         dimsorders(:,7)                                      , &
         c_dir_sse , c_prefix_sse, maxsize_sse, c_suffix_sse  , &
         nc_var_sse, c_none,nc_att_mask_sse, &
         sse(:,:,:,:))

    !==================== COMPUTE ZZ_WW (3D + time) =====================!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ ZZ_WW is computing ------'

    zz_ww(:,:,:,1) = 0._rprec

    ALLOCATE(depth_max(dims_reg(1,3),dims_reg(2,3)))
    call sub_memory(size(depth_max),'r','depth_max','sub_input_data_seq_symphonie')

    depth_max(:,:) = MAXVAL(ABS(zw0(:,:,1:kmax,1)), dim = 3)

    !!NG !! Save bathymetry for Symphonie data.
    !!NG  IF (ind_time(1) == 2) THEN
    !!NG    DO j = 1, dims_reg(2,3)
    !!NG      WRITE(22,'(340(1x,f0.2))') (depth_max(i,j), i= 1, dims_reg(1,3))
    !!NG    ENDDO
    !!NG  ENDIF

    ALLOCATE(ind_depth_max(dims_reg(1,3),dims_reg(2,3)))
    call sub_memory(size(ind_depth_max),'i','ind_depth_max','sub_input_data_seq_symphonie')

    ind_depth_max(:,:) = MAXLOC(ABS(zw0(:,:,1:kmax,1)), dim = 3)

    DO k = 1, kmax
       DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
             IF (depth_max(i,j) /= 0._rprec) THEN
                zz_ww(i,j,k,1) = sse(i,j,1,1) + zw0(i,j,k,1) * &
                     (1 + ( sse(i,j,1,1)/depth_max(i,j) ) )
             ENDIF
          ENDDO
       ENDDO
    ENDDO

    WRITE(lun_standard,*)' - ZZ_WW : max ', MAXVAL(zz_ww), ' min ', MINVAL(zz_ww)

    !=================== E3T (3D + time) ========================!
    !- Compute E3T from ZZ_WW -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ E3T scale factor is computing ------'
    IF (.NOT.ALLOCATED(e3t)) THEN
       ALLOCATE(e3t(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
       call sub_memory(size(e3t),'r','e3t','sub_input_data_seq_symphonie')
    ENDIF

    e3t(:,:,1:dims_reg(3,3)-1,:)= &
         zz_ww(:,:,1:dims_reg(3,3)-1,:) - zz_ww(:,:,2:dims_reg(3,3),:)


    e3t(:,:,:,:) = e3t(:,:,:,:) * tmask(:,:,:,:)

    !!  WRITE(0,*)MINLOC(e3t(:,:,:,:), mask = e3t(:,:,:,:) < 0._rprec)

    WRITE(lun_standard,*)' - E3T : max ', MAXVAL(e3t), ' min ', MINVAL(e3t)

    !========================= E3U ==============================!
    IF (.NOT.ALLOCATED(e3u)) THEN
       ALLOCATE(e3u(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
       call sub_memory(size(e3u),'r','e3u','sub_input_data_seq_symphonie')
    ENDIF
    e3u(1:dims_reg(1,3)-1,:,:,:) = &
         .5_rprec *( e3t(1:dims_reg(1,3)-1,:,:,:) + e3t(2:dims_reg(1,3),:,:,:) )

    e3u(dims_reg(1,3),:,:,:) = e3u(dims_reg(1,3)-1,:,:,:)

    !========================= E3V ==============================!
    IF (.NOT.ALLOCATED(e3v)) THEN
       ALLOCATE(e3v(dims_reg(1,3), dims_reg(2,3), dims_reg(3,3), 1))
       call sub_memory(size(e3v),'r','e3v','sub_input_data_seq_symphonie')
    ENDIF
    e3v(:,1:dims_reg(2,3)-1,:,:) = &
         .5_rprec *( e3t(:,1:dims_reg(2,3)-1,:,:) + e3t(:,2:dims_reg(2,3),:,:) )

    e3v(:,dims_reg(2,3),:,:) = e3v(:,dims_reg(2,3)-1,:,:)

    !====================== TRANSPORT ===========================!
    !=== U ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Zonal transport is computing ------'
    DO l = 1, 1
       DO k = 1, dims_reg(3,3)
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                uu(i,j,k,l) =  uu(i,j,k,l) * e3u(i,j,k,l) * e2u(i,j,1,1)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    WRITE(lun_standard,*)' - U Transport: max ', MAXVAL(uu), ' min ', MINVAL(uu)

    call sub_memory(-size(e3u),'r','e3u','sub_input_data_seq_symphonie')
    DEALLOCATE(e3u)


    !=== V ===!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Meridional transport is computing ------'
    DO l = 1, 1
       DO k = 1, dims_reg(3,3)
          DO j = 1, dims_reg(2,3)
             DO i = 1, dims_reg(1,3)
                vv(i,j,k,l) =  vv(i,j,k,l) * e3v(i,j,k,l) * e1v(i,j,1,1)
             ENDDO
          ENDDO
       ENDDO
    ENDDO
    WRITE(lun_standard,*)' - V Transport: max ', MAXVAL(vv), ' min ', MINVAL(vv)

    call sub_memory(-size(e3v),'r','e3v','sub_input_data_seq_symphonie')
    DEALLOCATE(e3v)


    !=== W ===!
    !- Must be no divergent -!
    WRITE(lun_standard,*)''
    WRITE(lun_standard,*)' ------ Vertical transport is computing ------'
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    DO l = 1,1
       DO k = dims_reg(3,3)-1, 1, -1
          DO j = 2, dims_reg(2,3)-1 
             DO i = 2, dims_reg(1,3)-1
                ww(i,j,k,l) = ww(i,j,k+1,l) &
                     + uu(i-1,   j, k, l) - uu(i, j, k, l) &
                     + vv(  i, j-1, k, l) - vv(i, j, k, l)
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    DO l = 1, 1
       DO j = 1, dims_reg(2,3)
          DO i = 1, dims_reg(1,3)
             DO k = 2, ind_depth_max(i,j) - 1
                ww(i,j,k,l) = &
                     +   ww(i,j,k,l) &
                     - ( ww(i,j,1,l) / ( zz_ww(i,j,1,l) - depth_max(i,j) ) ) &
                     *                 ( zz_ww(i,j,k,l) - depth_max(i,j) )
             ENDDO
          ENDDO
       ENDDO
    ENDDO

    ww(:,:, 1,:) = 0._rprec
    ww(:,:,dims_reg(3,3),:) = 0._rprec

    WRITE(lun_standard,*)' - W Transport: max ', MAXVAL(ww), ' min ', MINVAL(ww)

    zz_ww(:,:,:,:) = zw0(:,:,:,:)

    call sub_memory(-size(depth_max),'r','depth_max','sub_input_data_seq_symphonie')
    DEALLOCATE(depth_max)
    call sub_memory(-size(ind_depth_max),'i','ind_depth_max','sub_input_data_seq_symphonie')
    DEALLOCATE(ind_depth_max)
    call sub_memory(-size(zw0),'r','zw0','sub_input_data_seq_symphonie')
    DEALLOCATE(zw0)
    call sub_memory(-size(sse),'r','sse','sub_input_data_seq_symphonie')
    DEALLOCATE(sse)


  END SUBROUTINE sub_input_data_seq_symphonie
!!***
!=========================================================================
!!****f* mod_input_data/sub_input_data_main()
!! NAME
!!   sub_input_data_main()
!!
!! FUNCTION
!!   These subroutine reads 3d data time series in different netcdf 
!!   files.
  !!   We only assume that netcdf file_name are build on this form:
  !!                  file_name = prefix_XXXX_suffix
  !!   Where XXXX, is an integer from "ind0" to "indn" on "ndigits" (here 4).
  !!   
  !!   The number of time step in each netcdf file could be different.
  !!
  !!   If time series is not completed, program stop.
  !!
  !!   WARNING: Today, we assume that the first data time step correspond
  !!            to the first record in the netcdf file. (start_time = 1)
  !!   
  !!   * Begin loop on netcdf files for a variable (current or tracer):
  !!       - build file name,
  !!       - open netcdf file,
  !!       - read dimensions,
  !!       - verify if time dimension is reached,
  !!       - read data,
  !!       - If needed, read EIV current components,
  !!       - if needed, masked values are set to zero,
  !!       - close netcdf file
  !!   * END loop.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   Except "var", all these arguments are red from a namelist file.
  !!   * cdir   : netcdf file data directory
  !!   * prefix : prefix file name
  !!   * ind0   : first indice
  !!   * indn   : last indice
  !!   * ndigits: number of digits of indices
  !!   * suffix : suffix file name
  !!   * ncvar  : netcdf variable name
  !!   * nceiv  : netcdf EIV variable name for currents (or NONE)
  !!   * ncmask : netcdf mask name                      (or NONE)
  !!
  !!   * var    : 4d array where netcdf values will be affected
  !!
  !! TODO
  !!   * add the possibilty to start time step not at the first netcdf
  !!     file record.
  !!
  !! USES
  !!   * sub_build_filename      (mod_netcdf.f90)
  !!   * sub_open_netcdf_file    (mod_netcdf.f90)
  !!   * sub_select_var_dims     (mod_netcdf.f90)
  !!   * sub_read_netcdf_var4d   (mod_netcdf.f90)
  !!   * sub_read_netcdf_att_val (mod_netcdf.f90)
  !!   * sub_close_netcdf_file   (mod_netcdf.f90)
  !!
  !! USED BY
  !!   * sub_input_data (mod_input_data.f90)
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_main(cdir,prefix,ind0,indn,ndigits,suffix, & 24,15
       ncvar,nceiv,ncmask,var,mask, key_ijt)

    !-------------!
    ! Declaration !
    !-------------!
    !- arguments -!
    CHARACTER(len = *)    , INTENT(in) :: cdir, prefix, suffix,  ncvar, nceiv, ncmask
    INTEGER(kind  = iprec), INTENT(in) :: ind0, indn, ndigits
    REAL(kind     = rprec), DIMENSION(:,:,:,:), INTENT(out) :: var
    INTEGER(kind  = iprec), DIMENSION(:,:,:,:), INTENT(out), OPTIONAL :: mask
    LOGICAL, INTENT(in), OPTIONAL :: key_ijt


    !- local variables -!
    CHARACTER(len = 4), PARAMETER :: c_none='NONE' ! 
    INTEGER(kind = iprec) :: ii           ! Used for do loop
    INTEGER(kind = iprec) :: start_time ! First time step, dim. 4 of var
    INTEGER(kind = iprec) :: end_time   ! Last time step , dim. 4 of var
    INTEGER(kind = iprec) :: ncid         ! netcdf file ID
    INTEGER(kind = iprec) :: varid        ! netcdf variable ID
    INTEGER(kind = iprec) :: indmin       ! minimum indice
    INTEGER(kind = iprec) :: indmax       ! maximum indice
    INTEGER(kind = iprec) :: dimx         ! dimension in x (i)
    INTEGER(kind = iprec) :: dimy         ! dimension in y (j)
    INTEGER(kind = iprec) :: dimz         ! dimension in z (k)
    INTEGER(kind = iprec) :: dimt         ! dimension in t (l)
    INTEGER(kind = iprec) :: tmp_var_id   ! temporaly variable ID
    REAL(kind    = rprec) :: att_mask_val ! attribute mask value
    CHARACTER(len = 128)  :: c_filename   ! file name
    REAL(kind = rprec) :: sfval
    REAL(kind = rprec) :: aoval

    !!NG: 10/07/2009 ! tmp_array: temporaly array
    !!NG: 10/07/2009 REAL(kind    = rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: tmp_array

    ! l_esc: a logical to exit do loop if time series is completed before 
    !        the end of the netcdf file.
    LOGICAL                :: l_esc = .FALSE. ! a logical to exit do loop

    ! dimsorder: array where netcdf variable dimensions order are stored.
    !            This, because sometime dimensions order could be different 
    !            than (x, y, z, t). It could be (x, y, t, nothing) or 
    !            (z, t, nothing, nothing).
    !            More information: see sub_select_var_dims and 
    !            sub_read_netcdf_var4d in mod_netcdf.f90 file.
    INTEGER(kind = iprec), DIMENSION(4) :: dimsorder 

    !-------------!
    ! Code begins !
    !-------------!
    !- Initialization -!
    start_time = 1 ! first time step corresponds to the first netcdf record !
    end_time   = 0 ! end_time is intialized.
    l_esc      = .FALSE. ! l_esc is intialized.
    IF (PRESENT(mask)) THEN
       mask(:,:,:,:)=1_iprec
    ENDIF

    !- If there is only one netcdf file, ind0 must be < 0 (by default -1). -!
    !- and DoLoop is executed only once.                                   -!
    !- Else Doloop is executed indn - ind0 + 1 times.                      -!
    !- We test it here.                                                    -!
    IF (ind0 < 0) THEN
      indmin = -1
      indmax = -1
    ELSE
      indmin = ind0
      indmax = indn
    ENDIF

    !---------------------------!
    !- main DoLoop starts here -!
    !---------------------------!
    DO ii = indmin, indmax

      !- Build file name -!
      CALL sub_build_filename(ii, ndigits, prefix, suffix, c_filename)

      !- Open Netcdf file -!
      CALL sub_open_netcdf_file(cdir, c_filename, ncid)

      !- Read variable dimensions -!
      IF(PRESENT(key_ijt)) THEN
         CALL sub_select_var_dims(ncid, ncvar, varid, dimsorder, &
              dimx, dimy, dimz, dimt, key_ijt)
      ELSE
         CALL sub_select_var_dims(ncid, ncvar, varid, dimsorder, &
              dimx, dimy, dimz, dimt)
      ENDIF

      !- From the precedent call we increment end_time -!
      end_time = end_time + dimt

      !- We test here, if end_time is greater than the last time dimension value -!
      IF ( end_time >= lmt ) THEN !! NG: 10/07/2009
        end_time = lmt
        l_esc = .TRUE.
      ENDIF

      !-----------------------------------------------------------------!
      !- Read netcdf 3d variable and store it in our 4d array variable -!
      !- Fourth dimension is time -!------------------------------------!
      !----------------------------!
      !!NG: bug 15/04/2009 CALL sub_read_netcdf_var4d(ncid, varid, var(:,:,1:dimz,start_time:end_time), &
      !!NG: bug 15/04/2009     dimsorder, dims_reg)
      CALL sub_read_netcdf_var4d(ncid, varid,   &
           var(:,:,:,start_time:end_time), &
           dimsorder, dims_reg)

      !- Test if mask value must be read -!
      IF ( ncmask /= c_none) THEN

        !- Read mask value in netcdf file -!
        CALL sub_read_netcdf_att_val( &
             ncid                   , &
             varid                  , &
             ncmask                 , &
             att_mask_val           , &
             0._rprec               )

       !-----------------------------------!
       !- Scale_factor is red or set to 1 -!
       !-----------------------------------!
       CALL sub_read_netcdf_att_val( &
            ncid                   , &
            varid                  , &
            'scale_factor'         , &
            sfval                  , &
            1._rprec               )

       att_mask_val = att_mask_val * sfval

       !---------------------------------!
       !- Add_offset is red or set to 0 -!
       !---------------------------------!
       CALL sub_read_netcdf_att_val( &
            ncid                   , &
            varid                  , &
            'add_offset'           , &
            aoval                  , &
            0._rprec               )

       att_mask_val = att_mask_val + aoval

        !- Mask 4d array values with 0 -!
        IF (PRESENT(mask)) THEN
           WHERE(var(:,:,:, start_time) == att_mask_val) &
              mask(:,:,:,1)                   = 0_iprec
        ENDIF
        
        WHERE(var(:,:,:, start_time:end_time) == att_mask_val) &
             var(:,:,:, start_time:end_time) = 0._rprec

      ENDIF

      !-------------------------------------------------!
      !- Test if we have to read EIV array values and  -!
      !- if yes,  add it to current                    -!
      !-------------------------------------------------!
      IF ( nceiv /= c_none ) THEN

        !- Read variable dimensions -!
        CALL sub_select_var_dims(ncid, nceiv, tmp_var_id, dimsorder, &
             dimx, dimy, dimz, dimt)

        !- Dynamic allocation of tmp_array -!
        !!NG: 10/07/2009 modifications
        IF (.NOT.ALLOCATED(tmp_array)) THEN
           ALLOCATE(tmp_array(1:dimx, 1:dimy, 1:dimz, 1:dimt))
           CALL sub_memory(size(tmp_array),'r','tmp_array','sub_input_data_main')
        ELSEIF ( &
              (SIZE(tmp_array, dim=1) /= dimx).OR. &
              (SIZE(tmp_array, dim=2) /= dimy).OR. &
              (SIZE(tmp_array, dim=3) /= dimz).OR. &
              (SIZE(tmp_array, dim=4) /= dimt)) THEN
           CALL sub_memory(-size(tmp_array),'r','tmp_array', &
                'sub_input_data_main')
           DEALLOCATE(tmp_array)
           ALLOCATE(tmp_array(1:dimx, 1:dimy, 1:dimz, 1:dimt))
           CALL sub_memory(size(tmp_array),'r','tmp_array','sub_input_data_main')
        ENDIF
        !!NG: 10/07/2009 

        !- Read EIV values in tmp_array -!
        CALL sub_read_netcdf_var4d(ncid, tmp_var_id, tmp_array(:,:,:,:), &
             dimsorder, dims_reg)

        !-----------------------------------!
        !- Test if mask value must be read -!
        !-----------------------------------!
        IF ( ncmask /= c_none) THEN
          !- Read mask value in netcdf file -!

          CALL sub_read_netcdf_att_val( &
               ncid                   , &
               varid                  , &
               ncmask                 , &
               att_mask_val           , &
               0._rprec               )

          !- Mask 4d array values with 0 -!

          WHERE(tmp_array(:,:,:,:) == att_mask_val) &
               tmp_array(:,:,:,:) = 0._rprec
        ENDIF

        !- Add EIV component to current -! 
        var(:,:,:,start_time:end_time) = var(:,:,:,start_time:end_time) + &
             tmp_array(:,:,:,:)

        !- Deallocate memory of tmp_array -!
        !!NG: 10/07/2009 CALL sub_memory(-size(tmp_array),'r','tmp_array','sub_input_data_main')
        !!NG: 10/07/2009 DEALLOCATE(tmp_array)


      ENDIF

      !- Comments -!
      WRITE(lun_standard,*)' - ', TRIM(ncvar),': max ', MAXVAL(var), ' min ', MINVAL(var)

      !- Close netcdf file -!
      CALL sub_close_netcdf_file(ncid)

      !- Exit DoLoop if l_esc == .TRUE.
      IF (l_esc) EXIT

      !----------------------------------------------------------------!
      !- Start_time pointer is switched to stored correctly in the 4d -!
      !- array val the next time series                               -!
      !----------------------------------------------------------------!
      start_time = end_time + 1

    ENDDO

    !- We stop the program if the time series in not completed -!
    IF ( end_time < lmt ) STOP 'mod_input_data: mod_input_data_main: end_time < lmt!'

  END SUBROUTINE sub_input_data_main
  !!***
  !=========================================================================
  !!****f* mod_input_data_seq/sub_input_data_seq_main()
  !! NAME
  !!   sub_input_data_seq_main()
  !!
  !! FUNCTION
  !!   These subroutine reads 3d data time series in different netcdf 
  !!   files.
  !!   We only assume that netcdf file_name are build on this form:
  !!                  file_name = prefix_XXXX_suffix
  !!   Where XXXX, is an integer from "ind0" to "indn" on "ndigits" (here 4).
  !!   
  !!   The number of time step in each netcdf file could be different.
  !!
  !!   If time series is not completed, program stop.
  !!
  !!   WARNING: Today, we assume that the first data time step correspond
  !!            to the first record in the netcdf file. (start_time = 1)
  !!   
  !!   * Begin loop on netcdf files for a variable (current or tracer):
  !!       - build file name,
  !!       - open netcdf file,
  !!       - read dimensions,
  !!       - verify if time dimension is reached,
  !!       - read data,
  !!       - If needed, read EIV current components,
  !!       - if needed, masked values are set to zero,
  !!       - close netcdf file
  !!   * END loop.
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   Except "var", all these arguments are red from a namelist file.
  !!   * cdir   : netcdf file data directory
  !!   * prefix : prefix file name
  !!   * ind0   : first indice
  !!   * indn   : last indice
  !!   * ndigits: number of digits of indices
  !!   * suffix : suffix file name
  !!   * ncvar  : netcdf variable name
  !!   * nceiv  : netcdf EIV variable name for currents (or NONE)
  !!   * ncmask : netcdf mask name                      (or NONE)
  !!
  !!   * var    : 4d array where netcdf values will be affected
  !!
  !! TODO
  !!   * add the possibilty to start time step not at the first netcdf
  !!     file record.
  !!
  !! USES
  !!   * sub_build_filename      (mod_netcdf.f90)
  !!   * sub_open_netcdf_file    (mod_netcdf.f90)
  !!   * sub_select_var_dims     (mod_netcdf.f90)
  !!   * sub_read_netcdf_var4d   (mod_netcdf.f90)
  !!   * sub_read_netcdf_att_val (mod_netcdf.f90)
  !!   * sub_close_netcdf_file   (mod_netcdf.f90)
  !!
  !! USED BY
  !!   * sub_input_data_seq (mod_input_data_seq.f90)
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_seq_main(      & 24,16
       ncid, varid                       , &
       new_file, ind_file                , &
       ind_time, ind_time_size           , &
       dimsorder                         , &
       cdir, prefix, ndigits, suffix     , &
       ncvar, nceiv, ncmask, var, mask, key_ijt)

    !-------------!
    ! Declaration !
    !-------------!
    !- arguments -!
    INTEGER(kind = iprec) , INTENT(inout) :: ncid
    INTEGER(kind = iprec) , INTENT(inout) :: varid
    LOGICAL               , INTENT(inout) :: new_file
    INTEGER(kind = iprec) , INTENT(inout) :: ind_file
    INTEGER(kind = iprec) , INTENT(inout) :: ind_time
    INTEGER(kind = iprec) , INTENT(out)   :: ind_time_size
    INTEGER(kind  = iprec), DIMENSION(:,:,:,:), INTENT(out), OPTIONAL :: mask
    LOGICAL, INTENT(in), OPTIONAL :: key_ijt

    ! dimsorder: array where netcdf variable dimensions order are stored.
    !            This, because sometime dimensions order could be different 
    !            than (x, y, z, t). It could be (x, y, t, nothing) or 
    !            (z, t, nothing, nothing).
    !            More information: see sub_select_var_dims and 
    !            sub_read_netcdf_var4d in mod_netcdf.f90 file.
    INTEGER(kind = iprec), DIMENSION(4), INTENT(out) :: dimsorder 

    INTEGER(kind = iprec), DIMENSION(4) :: tmp_dimsorder 

    CHARACTER(len = *)    , INTENT(in)    :: cdir, prefix, suffix, ncvar, nceiv, ncmask
    INTEGER(kind  = iprec), INTENT(in)    :: ndigits
    REAL(kind     = rprec), DIMENSION(:,:,:,:), INTENT(out) :: var

    REAL(kind = rprec) :: sfval
    REAL(kind = rprec) :: aoval

    !- local variables -!
    CHARACTER(len = 4), PARAMETER :: c_none='NONE' ! 
    INTEGER(kind = iprec) :: dimx         ! dimension in x (i)
    INTEGER(kind = iprec) :: dimy         ! dimension in y (j)
    INTEGER(kind = iprec) :: dimz         ! dimension in z (k)
    INTEGER(kind = iprec) :: dimt         ! dimension in t (l)
    INTEGER(kind = iprec) :: tmp_var_id   ! temporaly variable ID
    REAL(kind    = rprec) :: att_mask_val ! attribute mask value
    CHARACTER(len = 128)  :: c_filename   ! file name

    !!NG:10/07/2009 ! tmp_array: temporaly array
    !!NG:10/07/2009 REAL(kind    = rprec), DIMENSION(:,:,:,:), ALLOCATABLE :: tmp_array

    !-------------!
    ! Code begins !
    !-------------!
    IF (PRESENT(mask)) THEN
       mask(:,:,:,:)=1_iprec
    ENDIF

    IF (new_file) THEN
       !- Build file name -!
       CALL sub_build_filename(ind_file, ndigits, prefix, suffix, c_filename)

       !- Open Netcdf file -!
       CALL sub_open_netcdf_file(cdir, c_filename, ncid)
       new_file = .FALSE.

       !- Read variable dimensions -!
       IF(PRESENT(key_ijt)) THEN
          CALL sub_select_var_dims(ncid, ncvar, varid, dimsorder, &
               dimx, dimy, dimz, dimt, key_ijt)
       ELSE
          CALL sub_select_var_dims(ncid, ncvar, varid, dimsorder, &
               dimx, dimy, dimz, dimt)
       ENDIF

       ind_time_size = dimt
       IF (TRIM(forback) == 'backward'.AND.(ind_time<1)) ind_time = dimt

    ENDIF

    !-----------------------------------------------------------------!
    !- Read netcdf 3d variable and store it in our 4d array variable -!
    !- Fourth dimension is time -!------------------------------------!
    !----------------------------!
    !! NG: 19/01/2009 write(lun_standard,*)'dimsorder ',dimsorder
    !! NG: 19/01/2009 write(lun_standard,*)'dims_reg ',dims_reg
    CALL sub_read_netcdf_var4d(ncid, varid, var(:,:,:,:), &
         dimsorder, dims_reg(:,:), ind_time)

    !- Test if mask value must be read -!
    IF ( ncmask /= c_none) THEN

       !- Read mask value in netcdf file -!
       CALL sub_read_netcdf_att_val( &
            ncid                   , &
            varid                  , &
            ncmask                 , &
            att_mask_val           , &
            0._rprec               )

       !-----------------------------------!
       !- Scale_factor is red or set to 1 -!
       !-----------------------------------!
       CALL sub_read_netcdf_att_val( &
            ncid                   , &
            varid                  , &
            'scale_factor'         , &
            sfval                  , &
            1._rprec               )

       att_mask_val = att_mask_val * sfval

       !---------------------------------!
       !- Add_offset is red or set to 0 -!
       !---------------------------------!
       CALL sub_read_netcdf_att_val( &
            ncid                   , &
            varid                  , &
            'add_offset'           , &
            aoval                  , &
            0._rprec               )

       att_mask_val = att_mask_val + aoval

       !- Mask 4d array values with 0 -!
       IF (PRESENT(mask)) THEN
          WHERE(var(:,:,:,1) == att_mask_val) &
               mask(:,:,:,1)                   = 0_iprec
       ENDIF

       !- Mask 4d array values with 0 -!
       WHERE(var(:,:,:,:) == att_mask_val) &
            var(:,:,:,:) = 0._rprec
    ENDIF

    !-------------------------------------------------!
    !- Test if we have to read EIV array values and  -!
    !- if yes,  add it to current                    -!
    !-------------------------------------------------!
    IF ( nceiv /= c_none ) THEN

       !- Read variable dimensions -!
       CALL sub_select_var_dims(ncid, nceiv, tmp_var_id, tmp_dimsorder, &
            dimx, dimy, dimz, dimt)

       !!NG: 10/07/2009 !- Dynamic allocation of tmp_array -!
       !!NG: 10/07/2009 ALLOCATE(tmp_array(1:dimx, 1:dimy, 1:dimz, 1:it_ind))
       !!NG: 10/07/2009 CALL sub_memory(size(tmp_array),'r','tmp_array','sub_input_data_seq_main')

       !- Dynamic allocation of tmp_array -!
       !!NG: 10/07/2009 modifications
       IF (.NOT.ALLOCATED(tmp_array)) THEN
          ALLOCATE(tmp_array(1:dimx, 1:dimy, 1:dimz, 1:it_ind))
          CALL sub_memory(size(tmp_array),'r','tmp_array','sub_input_data_seq_main')
       ELSEIF ( &
            (SIZE(tmp_array, dim=1) /= dimx).OR. &
            (SIZE(tmp_array, dim=2) /= dimy).OR. &
            (SIZE(tmp_array, dim=3) /= dimz).OR. &
            (SIZE(tmp_array, dim=4) /= it_ind)) THEN
          CALL sub_memory(-size(tmp_array),'r','tmp_array', &
               'sub_input_data_main')
          DEALLOCATE(tmp_array)
          ALLOCATE(tmp_array(1:dimx, 1:dimy, 1:dimz, 1:it_ind))
          CALL sub_memory(size(tmp_array),'r','tmp_array','sub_input_data_seq_main')
       ENDIF
       !!NG: 10/07/2009 

       !- Read EIV values in tmp_array -!
       CALL sub_read_netcdf_var4d(ncid, tmp_var_id, tmp_array(:,:,:,:), &
            tmp_dimsorder, dims_reg, ind_time)

       !-----------------------------------!
       !- Test if mask value must be read -!
       !-----------------------------------!
       IF ( ncmask /= c_none) THEN
          !- Read mask value in netcdf file -!

          CALL sub_read_netcdf_att_val( &
               ncid                   , &
               tmp_var_id             , &
               ncmask                 , &
               att_mask_val           , &
               0._rprec               )

          !- Mask 4d array values with 0 -!

          WHERE(tmp_array(:,:,:,:) == att_mask_val) &
               tmp_array(:,:,:,:) = 0._rprec
       ENDIF

       !- Add EIV component to current -! 
       var(:,:,:,:) = var(:,:,:,:) + tmp_array(:,:,:,:)

       !!NG: 10/07/2009 !- Deallocate memory of tmp_array -!
       !!NG: 10/07/2009 CALL sub_memory(-size(tmp_array),'r','tmp_array','sub_input_data_seq_main')
       !!NG: 10/07/2009 DEALLOCATE(tmp_array)

    ENDIF

    !- Comments -!
    IF (id_comments) THEN
       WRITE(lun_standard,*)' - ', TRIM(ncvar),': max ', MAXVAL(var), ' min ', MINVAL(var)
    ENDIF

    IF (TRIM(forback) == 'forward' ) THEN
       ind_time = ind_time + 1

       IF (ind_time > ind_time_size) THEN
          !- Close netcdf file -!
          CALL sub_close_netcdf_file(ncid)
          new_file = .TRUE.
          ind_file = ind_file + 1
          ind_time = 1
       ENDIF

    ELSE !! Backward

       ind_time = ind_time - 1

       IF (ind_time < 1) THEN
          !- Close netcdf file -!
          CALL sub_close_netcdf_file(ncid)
          new_file = .TRUE.
          ind_file = ind_file - 1
       ENDIF

    ENDIF

  END SUBROUTINE sub_input_data_seq_main
  !!***


  !!****f* mod_input_data/sub_transp_alloc()
  !! NAME
  !!   sub_transp_alloc()
  !!
  !! FUNCTION
  !!   Dynamic allocation of the transport arrays (uu, vv, ww)
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * dim1: dimension 1, or x, or i
  !!   * dim1: dimension 2, or y, or j
  !!   * dim1: dimension 3, or z, or k
  !!   * dim1: dimension 4, or t, or l
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * mod_input_data (Private)
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_transp_alloc(dim1, dim2, dim3, dim4) 2,4

    !-------------!
    ! Declaration !
    !-------------!
    !- arguments -!
    INTEGER(kind=iprec), INTENT(in) :: dim1, dim2, dim3, dim4

    !-------------!
    ! Code begins !
    !-------------!
    !- Dynamic allocation -!
    ALLOCATE(uu(dim1, dim2, dim3, dim4))
    CALL sub_memory(size(uu),'r','uu','sub_transp_alloc')
    ALLOCATE(vv(dim1, dim2, dim3, dim4))
    CALL sub_memory(size(vv),'r','vv','sub_transp_alloc')
    ALLOCATE(ww(dim1, dim2, dim3, dim4))
    CALL sub_memory(size(ww),'r','ww','sub_transp_alloc')

    IF (TRIM(w_surf_option) == 'E-P-R') THEN
        ALLOCATE(epr(dim1, dim2, 1, dim4))
        CALL sub_memory(size(epr),'r','epr','sub_transp_alloc')
    ENDIF

  END SUBROUTINE sub_transp_alloc
  !!***
  !!****f* mod_input_data/sub_transp_dealloc()
  !! NAME
  !!   sub_transp_dealloc()
  !!
  !! FUNCTION
  !!   Deallocate memory of the transport arrays (uu, vv, ww)
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * No argument
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_transp_dealloc() 3,4

    !-------------!
    ! Code begins !
    !-------------!
    !- Deallocate arrays memory -!
    IF (ALLOCATED(uu)) THEN
       CALL sub_memory(-size(uu),'r','uu','sub_transp_dealloc')
       DEALLOCATE(uu)
    ENDIF
    IF (ALLOCATED(vv)) THEN
       CALL sub_memory(-size(vv),'r','vv','sub_transp_dealloc')
       DEALLOCATE(vv)
    ENDIF
    IF (ALLOCATED(ww)) THEN
       CALL sub_memory(-size(ww),'r','ww','sub_transp_dealloc')
       DEALLOCATE(ww)
    ENDIF

    IF (ALLOCATED(epr)) THEN
       CALL sub_memory(-size(epr),'r','epr','sub_transp_dealloc')
       DEALLOCATE(epr)
    ENDIF

  END SUBROUTINE sub_transp_dealloc
  !!***
  !=========================================================================
  !!****f* mod_input_data/sub_tracer_alloc()
  !! NAME
  !!   sub_tracer_alloc()
  !!
  !! FUNCTION
  !!   Dynamic allocation of the tracer arrays (tt, ss, rr)
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * dim1: dimension 1, or x, or i
  !!   * dim1: dimension 2, or y, or j
  !!   * dim1: dimension 3, or z, or k
  !!   * dim1: dimension 4, or t, or l
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * mod_input_data (Private)
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_tracer_alloc(dim1, dim2, dim3, dim4) 2,3

    !-------------!
    ! Declaration !
    !-------------!
    !- arguments -!
    INTEGER(kind=iprec), INTENT(in) :: dim1, dim2, dim3, dim4

    !-------------!
    ! Code begins !
    !-------------!
    !- Dynamic allocation -!
    ALLOCATE(tt(dim1, dim2, dim3, dim4))
    CALL sub_memory(size(tt),'r','tt','sub_tracer_alloc')

    ALLOCATE(ss(dim1, dim2, dim3, dim4))
    CALL sub_memory(size(ss),'r','ss','sub_tracer_alloc')

    ALLOCATE(rr(dim1, dim2, dim3, dim4))
    CALL sub_memory(size(rr),'r','rr','sub_tracer_alloc')

  END SUBROUTINE sub_tracer_alloc
  !!***

  !!****f* mod_input_data/sub_tracer_dealloc()
  !! NAME
  !!   sub_tracer_dealloc()
  !!
  !! FUNCTION
  !!   Deallocate memory of the tracer arrays (tt, ss, rr)
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * No argument
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_tracer_dealloc() 3,3

    !-------------!
    ! Code begins !
    !-------------!
    !- Deallocate arrays -!
    IF (ALLOCATED(tt)) THEN
       CALL sub_memory(-size(tt),'r','tt','sub_tracer_dealloc')
       DEALLOCATE(tt)
    ENDIF
    IF (ALLOCATED(ss)) THEN
       CALL sub_memory(-size(ss),'r','ss','sub_tracer_dealloc')
       DEALLOCATE(ss)
    ENDIF
    IF (ALLOCATED(rr)) THEN
       CALL sub_memory(-size(rr),'r','rr','sub_tracer_dealloc')
       DEALLOCATE(rr)
    ENDIF

  END SUBROUTINE sub_tracer_dealloc
  !!***
  !!****f* mod_input_data/sub_tracer_dealloc()
  !! NAME
  !!   sub_tracer_dealloc()
  !!
  !! FUNCTION
  !!   Deallocate memory of the tracer arrays (tt, ss, rr)
  !!
  !! AUTHOR
  !!   * Origin  : Nicolas Grima (April-May 2005)
  !! 
  !! CREATION DATE
  !!   * April-May 2005
  !!
  !! HISTORY
  !!   Date (dd/mm/yyyy/) - Modification(s)
  !!
  !! ARGUMENTS
  !!   * No argument
  !!
  !! TODO
  !!   
  !!
  !! USED BY
  !!   * trajec
  !!
  !! SOURCE
  !!=======================================================================

  SUBROUTINE sub_input_data_dealloc_mem() 2,5

    IF (ALLOCATED(zeta_a)) THEN
       CALL sub_memory(-size(zeta_a),'r','zeta_a','sub_input_data_dealloc_mem')
       DEALLOCATE(zeta_a)
    END IF

    IF (ALLOCATED(zw0)) THEN
       CALL sub_memory(-size(zw0),'r','zw0','sub_input_data_dealloc_mem')
       DEALLOCATE(zw0)
    END IF

    IF (ALLOCATED(e3u)) THEN
       CALL sub_memory(-size(e3u),'r','e3u','sub_input_data_dealloc_mem')
       DEALLOCATE(e3u)
    END IF

    IF (ALLOCATED(e3v)) THEN
       CALL sub_memory(-size(e3v),'r','e3v','sub_input_data_dealloc_mem')
       DEALLOCATE(e3v)
    END IF

    IF (ALLOCATED(tmp_array)) THEN
       CALL sub_memory(-size(tmp_array),'r','tmp_array','sub_input_data_dealloc_mem')
       DEALLOCATE(tmp_array)
    END IF

  END SUBROUTINE sub_input_data_dealloc_mem
  

END MODULE mod_input_data