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

Module mod_init_particules 9,6

  Use mod_precision
  USE mod_memory
  Use mod_namelist
  Use mod_lun
  Use mod_input_data
  Use mod_posin
  Use mod_seq

  Implicit None

  Integer( kind = iprec) :: & !
       maxsect = 0        , & ! quantitative mode = max of sections
       maxsegm = 0        , & ! quantitative mode = max of segments
       ntraj   = 0        , & ! number of trajectories
       nsect   = 0        , & ! number of sections
       nfnt    = 0        , & !
       icrit1  = 0

  Real(kind = rprec)     :: & !
       trtot   = 0._rprec  

Contains

  !=========================================================================
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !========================================================================
  ! four options to initialize the positions:
  ! - qualitative: on ASCII file "initial_positions.txt"
  ! - qualitative: on BINARY file "initial.bin"
  ! - quantitative: auto positionning (from file "segments")
  ! - quantitative: from file "initial.bin"
  !=========================================================================
  !=========================================================================


  Subroutine sub_init_particules_positions() 3,28

    Integer(kind = iprec) :: nb_dims
    Integer(kind = iprec), Dimension(1) :: dims

    !================!
    !- Declarations -!
    !================!
    Integer(kind=iprec) :: & !
         n               , & !
         ind             , & !
         nb_subset             , & !
         alloc_size      , & !
         ncid            , & !
         varid
    Integer(kind=iprec) :: is_loop


    Real(kind=rprec) ::     &  !
         trmin   =    1._rprec ! A particle is set if the transport is
                               ! higher thant this value. 
                               ! Reduce this value if your transport is
                               ! very small.


    Integer(kind=iprec), Dimension(:), Allocatable :: ind_subset


    Write(lun_standard,*)''
    Write(lun_standard,*)'==============================='
    Write(lun_standard,*)'= INITIAL PARTICULE POSITIONS ='
    Write(lun_standard,*)'==============================='
    Write(lun_standard,*)'       (read or computed)'

    !=====================!
    !- Quantitative MODE -!
    !=====================!
    If (Trim(mode) == 'quantitative') Then

       !-------------------------------------------------!
       !- Compute in sections.txt : maxsect and maxsegm -!
       !-------------------------------------------------!
       Call sub_init_part_maxsect_maxsegm()

       !---------------------!
       !- Auto positioning - !
       !---------------------!
       If (key_sequential) Then

          id_comments = .True. !- Write input data reading information -!

          Call posini_seq(trmin,ntraj,nsect,trtot,nfnt)

       Else

          Call posini(trmin, ntraj, nsect, trtot, nfnt)

       Endif

       !----------------------------------------
       ! initialization of transparent sections
       !----------------------------------------
       !!NG      IF (nfnt > 0) THEN
       !!NG        DO lu = 91, 90 + nfnt
       !!NG          WRITE(fname,'(a12,a15)')'transparent.',fntname(lu-90)
       !!NG          OPEN(lu,form='FORMATTED',file=fname)
       !!NG        END DO
       !!NG      ENDIF

    Endif

    !=====================================================!
    !- indices read on ASCII file (QUALITATIVE analysis) -!
    !=====================================================!
    If ((Trim(mode) == 'qualitative').And.(Trim(bin) == 'nobin')) Then

       Write(lun_output,*)
       Write(lun_output,*)'======================================================'
       Write(lun_output,*)'ASCII positions are read on file initial_positions.txt'
       Write(lun_output,*)'======================================================'

       Open(unit=lun_dummy, file='initial_positions.txt', action="read")

       Do n = 1, nmax + 1

          Read(lun_dummy,*,End=7054,ERR=7054)tfi(n),tfj(n),tfk(n),tfl(n),ttr(n)

       End Do

       Write(lun_error,*)
       Write(lun_error,*)'# of positions in initial.positions.txt exceeds nmax'
       Write(lun_error,*)'--- STOP ARIANE ---'
       Stop

7054   ntraj = n - 1

       Write(lun_output,*)
       Write(lun_output,*)'# of particles that could be read: ',ntraj

       If ( ntraj == 0) Then
          Write(lun_error,*)
          Write(lun_error,*)')-: The number of particule is 0 :-( '
          Write(lun_error,*)'Please verify that you have correctly set the file'
          Write(lun_error,*)'initial.positions.txt (param0 is obsoleted).'
          Write(lun_error,*)'--- STOP ARIANE ---'
       Endif

       If (key_sequential) 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 ((key_roms).OR.(key_symphonie))THEN
                alloc_size = 3
             ELSE
                IF (TRIM(w_surf_option) == 'E-P-R') THEN
                   alloc_size = 4
                ELSE
                   alloc_size = 3
                ENDIF
             ENDIF
          ENDIF

!!$          If (key_alltracers) Then
!!$             IF ((key_roms).OR.(key_symphonie))THEN
!!$                alloc_size = 7
!!$             ELSE
!!$                alloc_size = 6
!!$             ENDIF
!!$          Else
!!$             alloc_size = 3
!!$          Endif

          Call sub_seq_alloc(alloc_size)

       Endif

    Endif

    !===========================================!
    !- initial positions in file "initial.bin" -!
    !===========================================!
    If (Trim(bin) /= 'nobin') Then
       !-----------------------------
       ! indices read on BINARY file
       !-----------------------------
       !!NG      WRITE(lun_output,*)
       !!NG      WRITE(lun_output,*)'==============================================='
       !!NG      WRITE(lun_output,*)'BINARY positions are read from file initial.bin'
       !!NG      WRITE(lun_output,*)'==============================================='
       !!NG
       !!NG      OPEN(UNIT=lun_initial,file='initial.bin', form='UNFORMATTED', ACTION='READ')
       !!NG
       !!NG      DO n = 1, nmax + 1
       !!NG
       !!NG        READ(UNIT=lun_initial,END=7053,ERR=7053)tfi(n),tfj(n),tfk(n),tfl(n),ttr(n)
       !!NG
       !!NG        trtot = trtot + ttr(n)
       !!NG
       !!NG      END DO
       !!NG
       !!NG      WRITE(lun_error,*)
       !!NG      WRITE(lun_error,*)'# of available positions in initial.bin exceeds nmax'
       !!NG      WRITE(lun_error,*)'--- STOP ARIANE ---'
       !!NG      STOP
       !!NG
       !!NG7053  ntraj = n - 1

       Write(lun_output,*)
       Write(lun_output,*)'=================================================='
       Write(lun_output,*)'= Positions are read from ariane_initial.nc file ='
       Write(lun_output,*)'=================================================='

       !------------------------!
       !- Open the netcdf file -!
       !------------------------!
       Call sub_open_netcdf_file('.', 'ariane_initial.nc', ncid)
       Write(lun_standard, *)'  --- ariane_initial.nc is opened ---'

       !-------!
       !- TFI -! (rien à voir avec la chaine TF1 ;-)
       !-------!
       Call sub_read_netcdf_varid_ndims(     &
            ncid   = ncid                  , &
            nc_var = Trim(init_final)//'_x', &
            varid  = varid                 , &
            ndims  = nb_dims                 )

       Call sub_read_netcdf_var_dims( &
            ncid = ncid             , &
            varid = varid           , &
            ndims = nb_dims         , &
            dims  = dims(:)           )

       ntraj = dims(1)
       Write(lun_output,*)
       Write(lun_output,*)'# of particles that could be read: ',ntraj

       Call sub_read_netcdf_var1d( ncid, varid, tfi(1:ntraj))
       Write(lun_standard, *)'  - X positions are read from ', &
            &Trim(init_final)//'_x', ' -'

       !-------!
       !- TFJ -!
       !-------!
       Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_y', varid)

       Call sub_read_netcdf_var1d( ncid, varid, tfj(1:ntraj))
       Write(lun_standard, *)'  - Y positions are read from ', &
            &Trim(init_final)//'_y', ' -'

       !-------!
       !- TFK -!
       !-------!
       Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_z', varid)

       Call sub_read_netcdf_var1d( ncid, varid, tfk(1:ntraj))
       Write(lun_standard, *)'  - Z positions are read from ', &
            &Trim(init_final)//'_z', ' -'

       !-------!
       !- TFL -!
       !-------!
       Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_t', varid)

       Call sub_read_netcdf_var1d( ncid, varid, tfl(1:ntraj))
       Write(lun_standard, *)'  - Time positions read from ', &
            &Trim(init_final)//'_t', ' -'

       !-------!
       !- AGE -!
       !-------!
       If (key_read_age) Then
          Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_age', varid)

          Call sub_read_netcdf_var1d( ncid, varid, tage(1:ntraj))

          trtot = Sum(tage(1:ntraj)) / Real(ntraj, kind = rprec)

          Write(lun_standard, *)'  - Age are read from ', &
               &Trim(init_final)//'_age', ' - (mean =',  trtot, ')'
       Endif

       !-------!
       !- TTR -! (transport)
       !-------!
       Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_transp', varid)

       Call sub_read_netcdf_var1d( ncid, varid, ttr(1:ntraj))

       trtot = Sum(ttr(1:ntraj))

       Write(lun_standard, *)'  - Transports are read from ', &
            &Trim(init_final)//'_transp', ' - (total = ', trtot, ')'

       If (key_alltracers.And.(Trim(mode) == 'quantitative')) Then

          !---------------!
          !- TEMPERATURE -!
          !---------------!
          Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_temp', varid)

          Call sub_read_netcdf_var1d( ncid, varid, init_temp(1:ntraj))
          Write(lun_standard, *)'  - Temperature values are read from ', &
               &Trim(init_final)//'_temp', ' -'

          !------------!
          !- SALYNITY -!
          !------------!
          Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_salt', varid)

          Call sub_read_netcdf_var1d( ncid, varid, init_salt(1:ntraj))
          Write(lun_standard, *)'  - Salinity values are read from ', &
               &Trim(init_final)//'_salt', ' -'

          !-----------!
          !- DENSITY -!
          !-----------!
          Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_dens', varid)

          Call sub_read_netcdf_var1d( ncid, varid, init_dens(1:ntraj))
          Write(lun_standard, *)'  - Density values are read from ', &
               &Trim(init_final)//'_dens', ' -'

       Endif

       !-------------------------!
       !- Close the netcdf file -!
       !-------------------------!
       Call sub_close_netcdf_file(ncid)
       Write(lun_standard, *)'  --- ariane_initial.nc is closed ---'

    Endif

    !=========================!
    !- Read a subset of data -!
    !=========================!
    If (Trim(bin) == 'subbin') Then

       Allocate(ind_subset(nmax))
       call sub_memory(size(ind_subset),'i','ind_subset','sub_init_particules_positions')

       ind_subset(:)=0

       Open(lun_subset, file='subset.txt')

       Write(lun_output,*)'A subset of indices is reading on file subset.txt.'

       Do nb_subset = 1, ntraj+1
          Read(lun_subset,*,End=1010,ERR=1010)ind_subset(nb_subset)
       End Do

       Write(lun_error,*)
       Write(lun_error,*)'# of indices in subset exceeds ntraj'
       Write(lun_error,*)'--- STOP ARIANE ---'
       Stop

1010   Continue

       Close (lun_subset)

       nb_subset = nb_subset - 1

       If (Maxval(ind_subset(:)) > ntraj) Then
          Write(lun_error,*)'Error: mod_init_particules.f90'
          Write(lun_error,*)' An indice in the subset.txt file is'
          Write(lun_error,*)' higher than the number of particules !!!'
          Write(lun_error,*)'--- STOP ARIANE ---'
          Stop
       Endif

       trtot     = 0._rprec

       ! Because the indices of the subset particules increase
       ! we can not use an intermediate array.
       ! This means that ind-subset(ind) is always >= at ind.
       ! A little bit dangerous...
       Do ind = 1, nb_subset
          tfi(ind) = tfi(ind_subset(ind))
          tfj(ind) = tfj(ind_subset(ind))
          tfk(ind) = tfk(ind_subset(ind))
          tfl(ind) = tfl(ind_subset(ind))
          ttr(ind) = ttr(ind_subset(ind))
          tage(ind) = tage(ind_subset(ind))
          trtot    = trtot + ttr(ind)
          if (key_alltracers.And.(Trim(mode) == 'quantitative')) THEN
            init_temp(ind) = init_temp(ind_subset(ind))
            init_salt(ind) = init_salt(ind_subset(ind))
            init_dens(ind) = init_dens(ind_subset(ind))
          endif
       End Do


       ntraj = nb_subset

       Write(lun_output,*)
       Write(lun_output,*)'# of particles that were kept as a subset: ', ntraj

       call sub_memory(-size(ind_subset),'i','ind_subset','sub_init_particules_positions')
       Deallocate(ind_subset)


    Endif

    !========================================!
    !- Write initial positions in ASCI file -!
    !========================================!
    If ((Trim(mode) == 'quantitative').And.key_ascii_outputs ) Then

       Write(lun_standard,*)'  - Writing initial positions...'

       Do is_loop = 1, ntraj
          Write(lun_init_pos,7050) &
               tfi(is_loop),tfj(is_loop),tfk(is_loop),tfl(is_loop), &
               ttr(is_loop)
       Enddo

       Close(lun_init_pos)

    Endif

7050 Format(5(1x,f0.3))

    !---------------------------------------------!
    !- Compute lmin and lmax in qualitative mode -!
    !---------------------------------------------!
    If (Trim(mode) == 'qualitative') Then

       If(Trim(forback) == 'forward')  Then

          Write(lun_standard,*) &
               '  - Computing lmin in qualitative mode...'
          lmin = Nint(Minval(tfl, mask = tfl >  0  ))
          Write(lun_standard,*) '    - lmin =', lmin

          IF (lmin == 0) THEN
             Write(lun_error,*) &
                ' lmin = 0, there is a problem in the initial_positions.txt file'
             STOP
          ENDIF

       Elseif(Trim(forback) == 'backward') Then

          Write(lun_standard,*) &
               '  - Computing lmax in qualitative mode...'
          lmax = Nint(Maxval(tfl, mask = tfl < lmt+1))
          Write(lun_standard,*) '    - lmax =', lmax

       Endif

    Endif


    !======================================!
    !- Allocate arrays in sequential mode -!
    !======================================!
    If ((Trim(bin) /= 'nobin')) Then

       If (key_sequential) 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 ((key_roms).OR.(key_symphonie))THEN
                alloc_size = 3
             ELSE
                IF (TRIM(w_surf_option) == 'E-P-R') THEN
                   alloc_size = 4
                ELSE
                   alloc_size = 3
                ENDIF
             ENDIF
          ENDIF
!!$          If (key_alltracers) Then
!!$             IF ((key_roms).OR.(key_symphonie))THEN
!!$                alloc_size = 7
!!$             ELSE
!!$                alloc_size = 6
!!$             ENDIF
!!$          Else
!!$             alloc_size = 3
!!$          Endif

          Call sub_seq_alloc(alloc_size)
       Endif

    Endif

    !========================================================================
    ! I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - 
    ! [QUANT]
    !========================================================================
    If (Trim(mode) == 'quantitative') Then
       Write(lun_output,*)' '
       Write(lun_output,*)'QUANTITATIVE EXPERIMENT'
       Write(lun_output,*)'======================='
       Write(lun_output,*)'    minimum gridcell transport: ', trmin
       Write(lun_output,*)'  maximum individual transport: ', max_transport
       Write(lun_output,*)'      exact number of segments: ', maxsegm
       Write(lun_output,*)'     exact number of particles: ', ntraj
       Write(lun_output,*)'      exact number of sections: ', nsect
       Write(lun_output,*)'    total documented transport: ', trtot
       Write(lun_output,*)'number of transparent sections: ', nfnt
       Write(lun_output,*)' '
    Endif

  End Subroutine sub_init_particules_positions

  !=========================================================================
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !=========================================================================


  Subroutine sub_init_part_maxsect_maxsegm() 1

    Integer(kind=iprec) :: idum

    !-----------------------------------------------------------------!
    !- Open the file sections.txt and compute the number of segments -!
    !- and the number of sections.                                   -!
    !- (remember that a section could be composed by more than       -!
    !-  one segment)                                                 -!
    !-----------------------------------------------------------------!
    Open(unit=lun_dummy, file='sections.txt', action="read")
    Do 
       Read(lun_dummy,*, End=7777, ERR=7777) idum
       maxsegm = maxsegm + 1
       If (idum > maxsect) Then
          maxsect = idum
       Endif
    Enddo

7777 Continue

    Write(lun_output,*)'Number of segments in sections.txt file = ', maxsegm

    Close(lun_dummy)

  End Subroutine sub_init_part_maxsect_maxsegm

End Module mod_init_particules