!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! - 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.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!****f* ariane/posini()
!! NAME
!! posini() (posini.f90)
!!
!! USAGE
!!
!!
!! FUNCTION
!!
!!
!! AUTHOR
!! * Origin : Bruno Blanke (1992)
!! * F77toF90: Nicolas Grima (April-May 2005)
!!
!! CREATION DATE
!! * April-May 2005
!!
!! HISTORY
!! Date (dd/mm/yyyy/) - Modification(s)
!!
!! RESULT
!!
!!
!! EXAMPLES
!!
!!
!! NOTES
!! ROBODoc header style.
!!
!! TODO
!!
!!
!! PORTABILITY
!! Machine-OS - Fortran90/95 compiler
!! * i686-pc-linux-gnu - ifort
!! * i686-pc-linux-gnu - g95
!!
!! SEE ALSO
!!
!! TODO
!! * Should be transform to a module.
!! * Should be migrate to F90 !!!
!! * 3 similar loops, try to find a common base to reduce writing.
!! * more modularity
!! * more comments.
!!
!! USES
!! * USE mod_precision
!! * USE mod_namelist
!! * USE mod_input_grid
!! * USE mod_input_data
!! * USE mod_posin
!! * USE mod_flags
!! * USE mod_stats
!! * USE mod_stati
!! * USE mod_txt
!! * USE mod_fx
!! * USE mod_fy
!! * USE mod_fz
!! * USE mod_criter0
!! * USE mod_sigma
!! * USE mod_zinter
!!
!! USED BY
!! * trajec
!!
!! SOURCE
!!=========================================================================
SUBROUTINE posini(trmin, nb_traj,nb_sect,trans_total,pos_nfnt) 1,98
!-----------------------------------------------------------------------
! automatic positioning: particles are gathered where the transport is
! the largest, giving an equivalent individual tranport to all
! particles (not to exceed max_transport)
!------------------------------------------------------------------------
!------------------!
! USE ASSOCIAITION !
!------------------!
USE mod_precision
USE mod_namelist
USE mod_input_grid
USE mod_input_data
USE mod_posin
USE mod_flags
USE mod_stats
USE mod_stati
USE mod_txt
USE mod_fx
USE mod_fy
USE mod_fz
USE mod_criter0
USE mod_sigma
USE mod_zinter
USE mod_reducmem
USE mod_lun
USE mod_save_netcdf
!-------------!
! Declaration !
!-------------!
IMPLICIT NONE
!- arguments -!
INTEGER(kind=iprec), INTENT(out) :: nb_traj, nb_sect, pos_nfnt
REAL(kind=rprec) , INTENT(in) :: trmin
REAL(kind=rprec) , INTENT(out) :: trans_total
!- local variables -!
INTEGER(kind=iprec), DIMENSION(:), ALLOCATABLE :: &
it1, & !
it2, & !
jt1, & !
jt2, & !
kt1, & !
kt2, & !
num, & !
ndir !
INTEGER(kind=iprec) :: n !
INTEGER(kind=iprec) :: i,j,k
INTEGER(kind=iprec) :: is,js,ks
INTEGER(kind=iprec) :: it0s,jt0s,kt0s
INTEGER(kind=iprec) :: iu0s,jv0s,kw0s
INTEGER(kind=iprec) :: l !
INTEGER(kind=iprec) :: iu0 !
INTEGER(kind=iprec) :: it0 !
INTEGER(kind=iprec) :: it00 !
INTEGER(kind=iprec) :: n1 !
INTEGER(kind=iprec) :: n2 !
INTEGER(kind=iprec) :: i1 !
INTEGER(kind=iprec) :: j1 !
INTEGER(kind=iprec) :: k1 !
INTEGER(kind=iprec) :: l1 !
INTEGER(kind=iprec) :: jt0 !
INTEGER(kind=iprec) :: jt00 !
INTEGER(kind=iprec) :: ilim1 !
INTEGER(kind=iprec) :: ilim2 !
INTEGER(kind=iprec) :: jv0 !
INTEGER(kind=iprec) :: kt0 !
INTEGER(kind=iprec) :: kt00 !
INTEGER(kind=iprec) :: kw0 !
REAL(kind=rprec) :: u ! zonal current
REAL(kind=rprec) :: v ! meridional current
REAL(kind=rprec) :: w ! vertical current
REAL(kind=rprec) :: init_trans
REAL(kind=rprec) :: prof ! W depth
!- Comments -!
WRITE(lun_standard,*)''
WRITE(lun_standard,*)' ---------------------'
WRITE(lun_standard,*)' = ENTER >>>> POSINI ='
WRITE(lun_standard,*)' ---------------------'
!-------------!
! Code begins !
!-------------!
!- Dynamic Allocations -!
CALL sub_txt_alloc
()
CALL sub_txt_init
()
CALL sub_flags_alloc
(imt, jmt, kmt)
CALL sub_flags_init
(0)
CALL sub_stati_alloc
()
CALL sub_stati_init
(0)
CALL sub_stats_alloc
()
CALL sub_stats_init
()
ALLOCATE(it1(maxsegm))
CALL sub_memory
(size(it1),'i','it1','posini')
ALLOCATE(it2(maxsegm))
CALL sub_memory
(size(it2),'i','it2','posini')
ALLOCATE(jt1(maxsegm))
CALL sub_memory
(size(jt1),'i','jt1','posini')
ALLOCATE(jt2(maxsegm))
CALL sub_memory
(size(jt2),'i','jt2','posini')
ALLOCATE(kt1(maxsegm))
CALL sub_memory
(size(kt1),'i','kt1','posini')
ALLOCATE(kt2(maxsegm))
CALL sub_memory
(size(kt2),'i','kt2','posini')
ALLOCATE(num(maxsegm))
CALL sub_memory
(size(num),'i','num','posini')
ALLOCATE(ndir(maxsegm))
CALL sub_memory
(size(ndir),'i','ndir','posini')
!- Initializations -!
nb_traj = 0
pos_nfnt = 0
nb_sect = 0
trans_total = 0._rprec
!- Open sections.txt file -!
OPEN(unit=lun_dummy, file='sections.txt', action="read", &
POSITION='rewind')
!- DoLoop on segment number -!
DO n = 1, maxsegm
READ(lun_dummy,*) num(n),it1(n),it2(n),jt1(n),jt2(n),kt1(n),kt2(n),sname
ndir(n) = 0
WRITE(lun_output,*)num(n),' ',it1(n),' ',it2(n),' ',jt1(n),' ',jt2(n), &
' ',kt1(n),' ',kt2(n),' ',sname
IF (num(n) < 0) THEN
WRITE(lun_standard,*)''
WRITE(lun_standard,*)'transparent section associated with segment #',n
pos_nfnt = max0(pos_nfnt,-num(n))
ENDIF
IF (ABS(num(n)) > maxsect) THEN
WRITE(lun_error,*)'Problem with segment #' , n
WRITE(lun_error,*)'Maximum number of sections is: ', maxsect
STOP
ENDIF
IF ((it1(n).NE.it2(n)).AND. &
(jt1(n).NE.jt2(n)).AND. &
(kt1(n).NE.kt2(n))) THEN
WRITE(lun_error,*)'Problem with segment #',n
WRITE(lun_error,*)'One pair of indices must be identical.'
STOP
ENDIF
IF ( ( (it1(n) == it2(n)).AND.(jt1(n) == jt2(n)) ).OR. &
( (it1(n) == it2(n)).AND.(kt1(n) == kt2(n)) ).OR. &
( (jt1(n) == jt2(n)).AND.(kt1(n) == kt2(n)) ) ) THEN
WRITE(lun_error,*)'Problem with segment #',n
WRITE(lun_error,*)'Only one pair of indices can be identical'
STOP
ENDIF
IF (num(n) > 0) THEN
secname(num(n)) = sname
IF (num(n) > nb_sect) nb_sect = num(n)
ELSE
fntname(-num(n)) = sname
ENDIF
!-----------------------!
!- Latitudinal section -!
!-----------------------!
IF (it1(n) == it2(n)) THEN
ndir(n) = 1
i = it1(n)
IF (i < 0) THEN
i = -i
ndir(n) = -1
ENDIF
DO k = kt1(n), kt2(n)
DO j = jt1(n), jt2(n)
mtfin(i,j,k) = num(n)
END DO
END DO
ENDIF
IF (it1(n) < 0) it1(n) = -it1(n)
IF (it2(n) < 0) it2(n) = -it2(n)
!------------------------!
!- Longitudinal section -!
!------------------------!
IF (jt1(n) == jt2(n)) THEN
ndir(n) = 2
j = jt1(n)
IF (j < 0) THEN
j = -j
ndir(n) = -2
ENDIF
DO k = kt1(n), kt2(n)
DO i = it1(n), it2(n)
mtfin(i,j,k) = num(n)
END DO
END DO
ENDIF
IF (jt1(n) < 0) jt1(n) = -jt1(n)
IF (jt2(n) < 0) jt2(n) = -jt2(n)
!--------------------!
!- Vertical section -!
!--------------------!
IF (kt1(n) == kt2(n)) THEN
ndir(n) = 3
k = kt1(n)
IF (k < 0) THEN
k = -k
ndir(n) = -3
ENDIF
DO j = jt1(n), jt2(n)
DO i = it1(n), it2(n)
mtfin(i,j,k) = num(n)
END DO
END DO
ENDIF
IF (kt1(n) < 0) kt1(n) = -kt1(n)
IF (kt2(n) < 0) kt2(n) = -kt2(n)
END DO
secname(0)='meanders'
!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
!!===============================================================================!!
!!_______________________________________________________________________________!!
IF (TRIM(bin) == 'nobin') THEN
WRITE(lun_standard,*)'-- initial positions are automatically computed --'
WRITE(lun_standard,*)' bin = ', TRIM(bin)
!----------------!
!- Loop on time -!
!----------------!
DO l = lmin, lmax
!--------------------!
!- Loop on segments -!
!--------------------!
DO n = 1, maxsegm
!------------------------------------- --!
!- Test if the segment number = 1 -!
!- Segment where particules will start. -!
!------------------------------------- --!
IF (num(n) == 1) THEN
!!======================================================================
IF ((ndir(n) == 1).OR.(ndir(n) == -1)) THEN
it0 = it1(n)
iu0 = it0
it00 = it0 + 1
IF (ndir(n) < 0) THEN
iu0 = it0 - 1
it00 = it0 - 1
ENDIF
DO k = kt1(n), kt2(n)
DO j = jt1(n), jt2(n)
CALL sub_reducmem_shift_or_not_ind
(it0, j, k, it0s, js, ks)
IF ( (tmask(it0s,js,ks,1) > 0.5).AND. &
(mtfin(it00,j,k) <= 0 ) ) THEN
!!NG: bug 14/04/2009 IF (key_nointerpolstats) THEN
!!NG: bug 14/04/2009 IF(key_alltracers) THEN
!!NG: bug 14/04/2009 init_temp(nb_traj) = tt(it0s,js,ks,1)
!!NG: bug 14/04/2009 init_salt(nb_traj) = ss(it0s,js,ks,1)
!!NG: bug 14/04/2009 init_dens(nb_traj) = rr(it0s,js,ks,1)
!!NG: bug 14/04/2009 ENDIF ! key_alltracers
!!NG: bug 14/04/2009 ENDIF ! key_nointerpolstats
CALL sub_reducmem_shift_or_not_ind
(iu0, j, k, iu0s, js, ks)
u = uu(iu0s, js, ks, l)
!----------------------------------
! test on the sign of the velocity
!----------------------------------
IF ((ABS(u) > trmin).AND.(REAL(ndir(n), kind = rprec)*u > 0._rprec)) THEN
IF (lmt == 1) THEN
n1 = 1 + INT( SQRT(ABS(u)/max_transport) )
n2 = 1
ELSE
n1 = 1 + INT( (ABS(u)/max_transport)**(1._rprec/3._rprec) )
n2 = n1
ENDIF
init_trans = ABS(u)/REAL(n1,kind=rprec)/REAL(n1,kind=rprec)/REAL(n2,kind=rprec)
DO j1 = 1, n1
DO k1 = 1, n1
DO l1 = 1, n2
nb_traj = nb_traj + 1
IF (nb_traj > nmax) THEN
WRITE(lun_error,*)'too many particles, maximum is: ', nmax
STOP
ENDIF
tfi(nb_traj) = REAL(iu0, kind=rprec)
tfj(nb_traj) = REAL(j, kind=rprec) - 1._rprec - &
.5_rprec/REAL(n1, kind=rprec) + &
REAL(j1, kind=rprec)/REAL(n1, kind=rprec)
tfk(nb_traj) = REAL(k, kind=rprec) - &
.5_rprec/REAL(n1, kind=rprec) + &
REAL(k1, kind=rprec)/REAL(n1, kind=rprec)
tfl(nb_traj) = REAL(l, kind=rprec) - .5_rprec - &
.5_rprec/REAL(n2, kind=rprec) + &
REAL(l1, kind=rprec)/REAL(n2, kind=rprec)
ttr(nb_traj) = init_trans
IF (criter0
(tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
tfl(nb_traj),it0,j,k,l)) THEN
trans_total = trans_total + init_trans
isn(-1) = isn(-1) + 1
sn(-1) = sn(-1) + init_trans
IF (key_alltracers) THEN
IF (.NOT.key_nointerpolstats) THEN
init_temp(nb_traj) = zinter
(tt, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
init_salt(nb_traj) = zinter
(ss, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
IF (key_approximatesigma) THEN
init_dens(nb_traj) = zinter
(rr, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
ELSE
init_dens(nb_traj) = sigma
(zsigma, init_salt(nb_traj), init_temp(nb_traj))
ENDIF
ELSE !!NG: bug correction 14/04/2009
init_temp(nb_traj) = tt(it0s,js,ks,1) !!NG: bug correction 14/04/2009
init_salt(nb_traj) = ss(it0s,js,ks,1) !!NG: bug correction 14/04/2009
init_dens(nb_traj) = rr(it0s,js,ks,1) !!NG: bug correction 14/04/2009
ENDIF
IF (init_temp(nb_traj) <= s0min(5,-1)) s0min(5,-1) = init_temp(nb_traj)
IF (init_salt(nb_traj) <= s0min(6,-1)) s0min(6,-1) = init_salt(nb_traj)
IF (init_dens(nb_traj) <= s0min(7,-1)) s0min(7,-1) = init_dens(nb_traj)
IF (init_temp(nb_traj) >= s0max(5,-1)) s0max(5,-1) = init_temp(nb_traj)
IF (init_salt(nb_traj) >= s0max(6,-1)) s0max(6,-1) = init_salt(nb_traj)
IF (init_dens(nb_traj) >= s0max(7,-1)) s0max(7,-1) = init_dens(nb_traj)
s0x(5,-1) = s0x(5,-1) + init_temp(nb_traj) * init_trans
s0x(6,-1) = s0x(6,-1) + init_salt(nb_traj) * init_trans
s0x(7,-1) = s0x(7,-1) + init_dens(nb_traj) * init_trans
s0x2(5,-1) = s0x2(5,-1) + init_temp(nb_traj) * init_temp(nb_traj) * init_trans
s0x2(6,-1) = s0x2(6,-1) + init_salt(nb_traj) * init_salt(nb_traj) * init_trans
s0x2(7,-1) = s0x2(7,-1) + init_dens(nb_traj) * init_dens(nb_traj) * init_trans
ENDIF
IF (key_roms.OR.key_symphonie) THEN
prof = &
fz
(gi=tfi(nb_traj), gj=tfj(nb_traj), gk=-tfk(nb_traj), il=l)
ELSE
prof = fz
(gk=-tfk(nb_traj))
ENDIF
IF (fx(tfi(nb_traj),tfj(nb_traj)) <= s0min(1,-1)) &
s0min(1,-1) = fx
(tfi(nb_traj),tfj(nb_traj))
IF (fy(tfi(nb_traj),tfj(nb_traj)) <= s0min(2,-1)) &
s0min(2,-1) = fy
(tfi(nb_traj),tfj(nb_traj))
IF (-prof <= s0min(3,-1)) &
s0min(3,-1) = -prof
IF (fx(tfi(nb_traj),tfj(nb_traj)) >= s0max(1,-1)) &
s0max(1,-1) = fx
(tfi(nb_traj),tfj(nb_traj))
IF (fy(tfi(nb_traj),tfj(nb_traj)) >= s0max(2,-1)) &
s0max(2,-1) = fy
(tfi(nb_traj),tfj(nb_traj))
IF (-prof >= s0max(3,-1)) s0max(3,-1) = -prof
s0x(1,-1) = s0x(1,-1) + fx
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x(2,-1) = s0x(2,-1) + fy
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x(3,-1) = s0x(3,-1) - prof * init_trans
s0x2(1,-1) = s0x2(1,-1) + fx
(tfi(nb_traj),tfj(nb_traj)) * &
fx
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x2(2,-1) = s0x2(2,-1)+ fy
(tfi(nb_traj),tfj(nb_traj)) * &
fy
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x2(3,-1) = s0x2(3,-1)+ prof * prof * init_trans
ELSE
nb_traj=nb_traj-1
ENDIF
END DO
END DO
END DO
!
! fin de test sur signe de la vitesse
ENDIF
!
! fin de test sur masque/critere
ENDIF
!
! fins de boucle sur j,k,l
END DO
END DO
!!NG: END DO
!
! fin test it1=it2
ENDIF
!!======================================================================
IF ((ndir(n) == 2).OR.(ndir(n) == -2)) THEN
jt0 = jt1(n)
jv0 = jt0
jt00 = jt0 + 1
IF (ndir(n) < 0) THEN
jv0 = jt0 - 1
jt00 = jt0 - 1
ENDIF
DO k = kt1(n),kt2(n)
IF (key_periodic) THEN
ilim1 = MAX(2,it1(n))
ilim2 = MIN(imt-1,it2(n))
ELSE
ilim1 = it1(n)
ilim2 = it2(n)
ENDIF
DO i = ilim1, ilim2
CALL sub_reducmem_shift_or_not_ind
(i, jt0, k, is, jt0s, ks)
IF ( (tmask(is,jt0s,ks,1) > 0.5).AND.(mtfin(i,jt00,k) <= 0) ) THEN
!!NG: bug 14/04/2009 IF (key_nointerpolstats) THEN
!!NG: bug 14/04/2009 IF (key_alltracers) THEN
!!NG: bug 14/04/2009 init_temp(nb_traj) = tt(is,jt0s,ks,1)
!!NG: bug 14/04/2009 init_salt(nb_traj) = ss(is,jt0s,ks,1)
!!NG: bug 14/04/2009 init_dens(nb_traj) = rr(is,jt0s,ks,1)
!!NG: bug 14/04/2009 ENDIF
!!NG: bug 14/04/2009 ENDIF
CALL sub_reducmem_shift_or_not_ind
(i, jv0, k, is, jv0s, ks)
v = vv(is,jv0s,ks,l)
!----------------------------------
! test on the sign of the velocity
!----------------------------------
IF ((ABS(v) > trmin).AND.(REAL(ndir(n),kind=rprec)*v > 0._rprec)) THEN
IF (lmt == 1) THEN
n1 = 1 + INT( SQRT(ABS(v)/max_transport) )
n2 = 1
ELSE
n1 = 1 + INT( (ABS(v)/max_transport)**(1._rprec/3._rprec) )
n2 = n1
ENDIF
init_trans = ABS(v)/REAL(n1,kind=rprec)/REAL(n1,kind=rprec)/REAL(n2,kind=rprec)
DO i1 = 1,n1
DO k1 = 1,n1
DO l1 = 1,n2
nb_traj = nb_traj + 1
IF (nb_traj > nmax) THEN
WRITE(lun_error,*)'too many particles, maximum is: ',nmax
STOP
ENDIF
tfi(nb_traj) = REAL(i,kind=rprec) -1._rprec - &
.5_rprec/REAL(n1,kind=rprec) + &
REAL(i1,kind=rprec)/REAL(n1,kind=rprec)
tfj(nb_traj) = REAL(jv0,kind=rprec)
tfk(nb_traj) = REAL(k,kind=rprec) - &
.5_rprec/REAL(n1,kind=rprec) + &
REAL(k1,kind=rprec)/REAL(n1,kind=rprec)
tfl(nb_traj) = REAL(l,kind=rprec) -.5_rprec - &
.5_rprec/REAL(n2,kind=rprec) + &
REAL(l1,kind=rprec)/REAL(n2,kind=rprec)
ttr(nb_traj) = init_trans
IF (criter0
(tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
tfl(nb_traj),i,jt0,k,l)) THEN
trans_total = trans_total + init_trans
isn(-1) = isn(-1) + 1
sn(-1) = sn(-1) + init_trans
IF (key_alltracers) THEN
IF (.NOT.key_nointerpolstats) THEN
init_temp(nb_traj) = zinter
(tt, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
init_salt(nb_traj) = zinter
(ss, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
IF (key_approximatesigma) THEN
init_dens(nb_traj) = zinter
(rr, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
ELSE
init_dens(nb_traj) = sigma
(zsigma, init_salt(nb_traj), init_temp(nb_traj))
ENDIF
ELSE !!NG: bug correction 14/04/2009
init_temp(nb_traj) = tt(is,jt0s,ks,1) !!NG: bug correction 14/04/2009
init_salt(nb_traj) = ss(is,jt0s,ks,1) !!NG: bug correction 14/04/2009
init_dens(nb_traj) = rr(is,jt0s,ks,1) !!NG: bug correction 14/04/2009
ENDIF
IF (init_temp(nb_traj) <= s0min(5,-1)) s0min(5,-1) = init_temp(nb_traj)
IF (init_salt(nb_traj) <= s0min(6,-1)) s0min(6,-1) = init_salt(nb_traj)
IF (init_dens(nb_traj) <= s0min(7,-1)) s0min(7,-1) = init_dens(nb_traj)
IF (init_temp(nb_traj) >= s0max(5,-1)) s0max(5,-1) = init_temp(nb_traj)
IF (init_salt(nb_traj) >= s0max(6,-1)) s0max(6,-1) = init_salt(nb_traj)
IF (init_dens(nb_traj) >= s0max(7,-1)) s0max(7,-1) = init_dens(nb_traj)
s0x(5,-1) = s0x(5,-1) + init_temp(nb_traj) * init_trans
s0x(6,-1) = s0x(6,-1) + init_salt(nb_traj) * init_trans
s0x(7,-1) = s0x(7,-1) + init_dens(nb_traj) * init_trans
s0x2(5,-1) = s0x2(5,-1) + init_temp(nb_traj) * init_temp(nb_traj) * init_trans
s0x2(6,-1) = s0x2(6,-1) + init_salt(nb_traj) * init_salt(nb_traj) * init_trans
s0x2(7,-1) = s0x2(7,-1) + init_dens(nb_traj) * init_dens(nb_traj) * init_trans
ENDIF
IF (key_roms.OR.key_symphonie) THEN
prof = &
fz
(gi = tfi(nb_traj), &
& gj = tfj(nb_traj), &
& gk = -tfk(nb_traj), &
& il = l )
ELSE
prof = fz
(gk=-tfk(nb_traj))
ENDIF
IF (fx(tfi(nb_traj),tfj(nb_traj)) <= s0min(1,-1)) &
s0min(1,-1) = fx
(tfi(nb_traj),tfj(nb_traj))
IF (fy(tfi(nb_traj),tfj(nb_traj)) <= s0min(2,-1)) &
s0min(2,-1) = fy
(tfi(nb_traj),tfj(nb_traj))
IF (-prof <= s0min(3,-1)) s0min(3,-1) = -prof
IF (fx(tfi(nb_traj),tfj(nb_traj)) >= s0max(1,-1)) &
s0max(1,-1) = fx
(tfi(nb_traj),tfj(nb_traj))
IF (fy(tfi(nb_traj),tfj(nb_traj)) >= s0max(2,-1)) &
s0max(2,-1) = fy
(tfi(nb_traj),tfj(nb_traj))
IF (-prof >= s0max(3,-1)) s0max(3,-1) = -prof
s0x(1,-1) = s0x(1,-1) + fx
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x(2,-1) = s0x(2,-1) + fy
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x(3,-1) = s0x(3,-1) - prof * init_trans
s0x2(1,-1) = s0x2(1,-1) + fx
(tfi(nb_traj),tfj(nb_traj)) * &
fx
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x2(2,-1) = s0x2(2,-1) + fy
(tfi(nb_traj),tfj(nb_traj)) * &
fy
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x2(3,-1) = s0x2(3,-1) + prof * prof *init_trans
ELSE
nb_traj = nb_traj - 1
ENDIF
END DO
END DO
END DO
!
! fin de test sur signe de la vitesse
ENDIF
!
! fin de test sur masque/critere
ENDIF
!
! fins de boucle sur i,k,l
END DO
END DO
!
! fin test jt1=jt2
ENDIF
!!======================================================================
IF ((ndir(n) == 3).OR.(ndir(n) == -3)) THEN
kt0 = kt1(n)
kw0 = kt0 + 1
kt00 = kt0 + 1
IF (ndir(n) < 0) THEN
kw0 = kt0
kt00 = kt0 - 1
ENDIF
IF (kt0 == 0) THEN
kt0 = 1
ENDIF
DO j = jt1(n), jt2(n)
DO i = it1(n), it2(n)
CALL sub_reducmem_shift_or_not_ind
(i, j, kt0, is, js, kt0s)
IF ( (tmask(is,js,kt0s,1) > 0.5).AND.(mtfin(i,j,kt00) <= 0) ) THEN
!!NG: bug 14/04/2009 IF (key_nointerpolstats) THEN
!!NG: bug 14/04/2009 IF (key_alltracers) THEN
!!NG: bug 14/04/2009 init_temp(nb_traj) = tt(is,js,kt0s,1)
!!NG: bug 14/04/2009 init_salt(nb_traj) = ss(is,js,kt0s,1)
!!NG: bug 14/04/2009 init_dens(nb_traj) = rr(is,js,kt0s,1)
!!NG: bug 14/04/2009 ENDIF
!!NG: bug 14/04/2009 ENDIF
CALL sub_reducmem_shift_or_not_ind
(i, j, kw0, is, js, kw0s)
w = ww(is,js,kw0s,l)
!----------------------------------
! test on the sign of the velocity
!----------------------------------
IF ((ABS(w) > trmin).AND.(REAL(ndir(n),kind=rprec)*w < 0._rprec)) THEN
IF (lmt == 1) THEN
n1 = 1 + INT( SQRT(ABS(w)/max_transport) )
n2 = 1
ELSE
n1 = 1 + INT( (ABS(w)/max_transport)**(1._rprec/3._rprec) )
n2 = n1
ENDIF
init_trans = ABS(w)/REAL(n1,kind=rprec)/REAL(n1,kind=rprec)/REAL(n2,kind=rprec)
DO j1 = 1,n1
DO i1 = 1,n1
DO l1 = 1,n2
nb_traj = nb_traj+1
IF (nb_traj > nmax) THEN
WRITE(lun_error,*)'too many particles, maximum is: ',nmax
STOP
ENDIF
tfi(nb_traj) = REAL(i,kind=rprec) - 1._rprec - &
.5_rprec/REAL(n1,kind=rprec) + &
REAL(i1,kind=rprec)/REAL(n1,kind=rprec)
tfj(nb_traj) = REAL(j,kind=rprec) - 1._rprec - &
.5_rprec/REAL(n1,kind=rprec) + &
REAL(j1,kind=rprec)/REAL(n1,kind=rprec)
tfk(nb_traj) = REAL(kw0,kind=rprec)
tfl(nb_traj) = REAL(l,kind=rprec) - .5_rprec - &
.5_rprec/REAL(n2,kind=rprec) + &
REAL(l1,kind=rprec)/REAL(n2,kind=rprec)
ttr(nb_traj) = init_trans
IF (criter0
(tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
tfl(nb_traj),i,j,kt0,l)) THEN
trans_total = trans_total + init_trans
isn(-1) = isn(-1) + 1
sn(-1) = sn(-1) + init_trans
IF (key_alltracers) THEN
IF (.NOT.key_nointerpolstats) THEN
init_temp(nb_traj) = zinter
(tt, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
init_salt(nb_traj) = zinter
(ss, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
IF (key_approximatesigma) THEN
init_dens(nb_traj) = zinter
(rr, &
tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj),tfl(nb_traj))
ELSE
init_dens(nb_traj) = sigma
(zsigma, init_salt(nb_traj), init_temp(nb_traj))
ENDIF
ELSE !!NG: bug correction 14/04/2009
init_temp(nb_traj) = tt(is,js,kt0s,1) !!NG: bug correction 14/04/2009
init_salt(nb_traj) = ss(is,js,kt0s,1) !!NG: bug correction 14/04/2009
init_dens(nb_traj) = rr(is,js,kt0s,1) !!NG: bug correction 14/04/2009
ENDIF
IF (init_temp(nb_traj) <= s0min(5,-1)) &
s0min(5,-1) = init_temp(nb_traj)
IF (init_salt(nb_traj) <= s0min(6,-1)) &
s0min(6,-1) = init_salt(nb_traj)
IF (init_dens(nb_traj) <= s0min(7,-1)) &
s0min(7,-1) = init_dens(nb_traj)
IF (init_temp(nb_traj) >= s0max(5,-1)) &
s0max(5,-1) = init_temp(nb_traj)
IF (init_salt(nb_traj) >= s0max(6,-1)) &
s0max(6,-1) = init_salt(nb_traj)
IF (init_dens(nb_traj) >= s0max(7,-1)) &
s0max(7,-1) = init_dens(nb_traj)
s0x(5,-1) = s0x(5,-1) + init_temp(nb_traj) * init_trans
s0x(6,-1) = s0x(6,-1) + init_salt(nb_traj) * init_trans
s0x(7,-1) = s0x(7,-1) + init_dens(nb_traj) * init_trans
s0x2(5,-1) = s0x2(5,-1) + init_temp(nb_traj) * init_temp(nb_traj) * init_trans
s0x2(6,-1) = s0x2(6,-1) + init_salt(nb_traj) * init_salt(nb_traj) * init_trans
s0x2(7,-1) = s0x2(7,-1) + init_dens(nb_traj) * init_dens(nb_traj) * init_trans
ENDIF
IF (key_roms.OR.key_symphonie) THEN
prof= &
fz
(gi=tfi(nb_traj), gj=tfj(nb_traj), gk=-tfk(nb_traj), il=l)
ELSE
prof=fz
(gk=-tfk(nb_traj))
ENDIF
IF (fx(tfi(nb_traj),tfj(nb_traj)) <= s0min(1,-1)) &
s0min(1,-1) = fx
(tfi(nb_traj),tfj(nb_traj))
IF (fy(tfi(nb_traj),tfj(nb_traj)) <= s0min(2,-1)) &
s0min(2,-1) = fy
(tfi(nb_traj),tfj(nb_traj))
IF (-prof <= s0min(3,-1)) s0min(3,-1) = -prof
IF (fx(tfi(nb_traj),tfj(nb_traj)) >= s0max(1,-1)) &
s0max(1,-1) = fx
(tfi(nb_traj),tfj(nb_traj))
IF (fy(tfi(nb_traj),tfj(nb_traj)) >= s0max(2,-1)) &
s0max(2,-1) = fy
(tfi(nb_traj),tfj(nb_traj))
IF (-prof>= s0max(3,-1)) s0max(3,-1) = -prof
s0x(1,-1) = s0x(1,-1) + fx
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x(2,-1) = s0x(2,-1) + fy
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x(3,-1) = s0x(3,-1) - prof * init_trans
s0x2(1,-1) = s0x2(1,-1) + &
fx
(tfi(nb_traj),tfj(nb_traj)) * &
fx
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x2(2,-1) = s0x2(2,-1) + &
fy
(tfi(nb_traj),tfj(nb_traj)) * &
fy
(tfi(nb_traj),tfj(nb_traj)) * init_trans
s0x2(3,-1) = s0x2(3,-1) + prof * prof * init_trans
ELSE
nb_traj = nb_traj - 1
ENDIF
END DO
END DO
END DO
!
! fin de test sur signe de la vitesse
ENDIF
!
! fin de test sur masque/critere
ENDIF
!
! fins de boucle sur i,j,l
END DO
END DO
!
! fin test kt1=kt2
ENDIF
!
! fin test num(n)
ENDIF
!
! fin du positionnement
ENDDO
ENDDO
ELSE
WRITE(lun_standard,*) '-- posini: bin = ', TRIM(bin), ' --'
ENDIF
!- Deallocate dynamic memory array
CALL sub_memory
(-size(it1),'i','it1','posini')
DEALLOCATE(it1)
CALL sub_memory
(-size(it2),'i','it2','posini')
DEALLOCATE(it2)
CALL sub_memory
(-size(jt1),'i','jt1','posini')
DEALLOCATE(jt1)
CALL sub_memory
(-size(jt2),'i','jt2','posini')
DEALLOCATE(jt2)
CALL sub_memory
(-size(kt1),'i','kt1','posini')
DEALLOCATE(kt1)
CALL sub_memory
(-size(kt2),'i','kt2','posini')
DEALLOCATE(kt2)
CALL sub_memory
(-size(num),'i','num','posini')
DEALLOCATE(num)
CALL sub_memory
(-size(ndir),'i','ndir','posini')
DEALLOCATE(ndir)
WRITE(lun_standard,*)''
WRITE(lun_standard,*)' ----------------------'
WRITE(lun_standard,*)' = POSINI >>>> EXIT ='
WRITE(lun_standard,*)' ----------------------'
RETURN
END SUBROUTINE posini
!!***