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