!!****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(nstat, trmin, maxsect,t0,nsegm,ntraj,nsect,trtot,nfnt,ibin) 2,172
!-----------------------------------------------------------------------
! automatic positioning: particles are gathered where the transport is
! the largest, giving an equivalent individual tranport to all
! particles (not to exceed t0)
!------------------------------------------------------------------------
!------------------!
! 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
!-------------!
! Declaration !
!-------------!
IMPLICIT NONE
!- arguments -!
INTEGER(kind=iprec), INTENT(in) :: ibin, nstat
INTEGER(kind=iprec), INTENT(out) :: maxsect,nsegm, ntraj, nsect, nfnt
REAL(kind=rprec) , INTENT(in) :: t0, trmin
REAL(kind=rprec) , INTENT(out) :: trtot
!- local variables -!
INTEGER(kind=iprec), DIMENSION(:), ALLOCATABLE :: &
it1, & !
it2, & !
jt1, & !
jt2, & !
kt1, & !
kt2, & !
num, & !
ndir !
INTEGER(kind=iprec) :: n !
INTEGER(kind=iprec) :: maxsegm, idum !
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) :: t ! temperature
REAL(kind=rprec) :: s ! salinity
REAL(kind=rprec) :: r ! density
REAL(kind=rprec) :: u ! zonal current
REAL(kind=rprec) :: v ! meridional current
REAL(kind=rprec) :: w ! vertical current
REAL(kind=rprec) :: tr ! trajectory
REAL(kind=rprec) :: prof ! W depth
!-------------!
! Code begins !
!-------------!
!- Open the file param0 and compute the number of segments -!
!- and the number of sections. -!
!- (remember that a section could be composed by more than -!
!- one segment) -!
maxsegm = 1
maxsect = 1
OPEN(unit=lun_dummy, file='param0', action="read")
DO
READ(lun_dummy,*, END=7777, ERR=7777) idum
maxsegm = maxsegm + 1
IF (idum > maxsect) THEN
maxsect = maxsect + 1
ENDIF
ENDDO
7777 continue
CLOSE(lun_dummy)
!- Dynamic Allocations -!
CALL sub_txt_alloc
(maxsect, nstat)
CALL sub_txt_init
()
CALL sub_flags_alloc
(imt, jmt, kmt)
CALL sub_flags_init
(0)
CALL sub_stati_alloc
(maxsect)
CALL sub_stati_init
(0)
CALL sub_stats_alloc
(nstat, maxsect)
CALL sub_stats_init
()
ALLOCATE(it1(maxsegm))
ALLOCATE(it2(maxsegm))
ALLOCATE(jt1(maxsegm))
ALLOCATE(jt2(maxsegm))
ALLOCATE(kt1(maxsegm))
ALLOCATE(kt2(maxsegm))
ALLOCATE(num(maxsegm))
ALLOCATE(ndir(maxsegm))
!- Initializations -!
ntraj = 0
nfnt = 0
nsect = 0
trtot = 0.
OPEN(unit=lun_dummy, file='param0', action="read")
!- DoLoop on segment number -!
DO n=1,maxsegm+1
READ(lun_dummy,*,END=7054,ERR=7054) &
num(n),it1(n),it2(n),jt1(n),jt2(n),kt1(n),kt2(n),sname
ndir(n)=0
WRITE(58,*)num(n),' ',it1(n),' ',it2(n),' ',jt1(n),' ',jt2(n), &
' ',kt1(n),' ',kt2(n),' ',sname
IF (num(n) < 0) THEN
WRITE(58,*)'transparent section associated with segment #',n
nfnt=max0(nfnt,-num(n))
ENDIF
IF (ABS(num(n)) > maxsect) THEN
WRITE(0,*)'problem with segment #',n
WRITE(0,*)' 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(0,*)'problem with segment #',n
WRITE(0,*)' 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(0,*)'problem with segment #',n
WRITE(0,*)' only one pair of indices can be identical'
STOP
ENDIF
IF (num(n) > 0) THEN
secname(num(n))=sname
IF (num(n) > nsect) nsect=num(n)
ELSE
fntname(-num(n))=sname
ENDIF
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)
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)
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
WRITE(0,*)'# of segments to be read exceeds maxsegm'
STOP
7054 nsegm=n-1
WRITE(0,*)'# of segments that could be read: ',nsegm
secname(0)='meanders'
WRITE(58,*)'nsegm= ',nsegm
IF (ibin < 1) THEN
WRITE(0,*)'initial positions are automatically computed'
DO n = 1, nsegm
IF (num(n) == 1) THEN
IF ((ndir(n) == 1).OR.(ndir(n) == -1)) THEN
it0=it1(n)
iu0=it0 !NG Performance !
it00=it0+1 !
IF (ndir(n) < 0) THEN ! iu0 = it0 + INT( ( 1.- REAL(ndir(n), kind=rprec) ) * 0.5)
iu0=it0-1 ! it00 = it0 + SIGN(1, ndir(n))
it00=it0-1 !
ENDIF !NG
DO l=1,lmt
DO k=kt1(n),kt2(n)
DO j=jt1(n),jt2(n)
CALL sub_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
IF (key_nointerpolstats) THEN
IF (key_alltracers) THEN
t=tt(it0s,js,ks,l)
s=ss(it0s,js,ks,l)
r=rr(it0s,js,ks,l)
ENDIF
ENDIF
CALL sub_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.(float(ndir(n))*u > 0)) THEN
IF (lmt == 1) THEN
n1=1+INT( SQRT(ABS(u)/t0) )
n2=1
ELSE
n1=1+INT( (ABS(u)/t0)**(1./3.) )
n2=n1
ENDIF
tr=ABS(u)/n1/n1/n2
DO j1=1,n1
DO k1=1,n1
DO l1=1,n2
ntraj=ntraj+1
IF (ntraj > nmax) THEN
WRITE(0,*)'too many particles, maximum is: ',nmax
STOP
ENDIF
tfi(ntraj)=iu0
tfj(ntraj)=j-1.-.5/n1+float(j1)/n1
tfk(ntraj)=k-.5/n1+float(k1)/n1
tfl(ntraj)=l-.5-.5/n2+float(l1)/n2
ttr(ntraj)=tr
IF (criter0
(tfi(ntraj),tfj(ntraj),-tfk(ntraj), &
tfl(ntraj),it0,j,k,l)) THEN
trtot=trtot+tr
isn(-1)=isn(-1)+1
sn(-1)=sn(-1)+tr
IF (key_alltracers) THEN
IF (.NOT.key_nointerpolstats) THEN
t=zinter &
(tt,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
s=zinter &
(ss,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
IF (key_approximatesigma) THEN
r=zinter &
(rr,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
ELSE
r=sigma
(zsigma,s,t)
ENDIF
ENDIF
IF (t <= s0min(5,-1)) s0min(5,-1)=t
IF (s <= s0min(6,-1)) s0min(6,-1)=s
IF (r <= s0min(7,-1)) s0min(7,-1)=r
IF (t >= s0max(5,-1)) s0max(5,-1)=t
IF (s >= s0max(6,-1)) s0max(6,-1)=s
IF (r >= s0max(7,-1)) s0max(7,-1)=r
s0x(5,-1)=s0x(5,-1)+t*tr
s0x(6,-1)=s0x(6,-1)+s*tr
s0x(7,-1)=s0x(7,-1)+r*tr
s0x2(5,-1)=s0x2(5,-1)+t*t*tr
s0x2(6,-1)=s0x2(6,-1)+s*s*tr
s0x2(7,-1)=s0x2(7,-1)+r*r*tr
!!NG IF (key_outgridt) THEN
!!NG WRITE(50,7050) &
!!NG tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr, &
!!NG t,s,r
!!NG ELSE
WRITE(50,7050) &
tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr, &
t,s,r
!!NG ENDIF
ELSE
!!NG IF (key_outgridt) THEN
!!NG WRITE(50,7050) &
!!NG tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr
!!NG ELSE
WRITE(50,7050) &
tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr
!!NG ENDIF
ENDIF
IF (key_roms) THEN
prof=fz
(gi=tfi(ntraj), gj=tfj(ntraj), gk=-tfk(ntraj), il=l)
ELSE
prof=fz
(gk=-tfk(ntraj))
ENDIF
IF (fx(tfi(ntraj),tfj(ntraj)) <= s0min(1,-1)) &
s0min(1,-1)=fx
(tfi(ntraj),tfj(ntraj))
IF (fy(tfi(ntraj),tfj(ntraj)) <= s0min(2,-1)) &
s0min(2,-1)=fy
(tfi(ntraj),tfj(ntraj))
IF (-prof <= s0min(3,-1)) &
s0min(3,-1)=-prof
IF (fx(tfi(ntraj),tfj(ntraj)) >= s0max(1,-1)) &
s0max(1,-1)=fx
(tfi(ntraj),tfj(ntraj))
IF (fy(tfi(ntraj),tfj(ntraj)) >= s0max(2,-1)) &
s0max(2,-1)=fy
(tfi(ntraj),tfj(ntraj))
IF (-prof >= s0max(3,-1)) &
s0max(3,-1)=-prof
s0x(1,-1)=s0x(1,-1)+fx
(tfi(ntraj),tfj(ntraj))*tr
s0x(2,-1)=s0x(2,-1)+fy
(tfi(ntraj),tfj(ntraj))*tr
s0x(3,-1)=s0x(3,-1)- prof *tr
s0x2(1,-1)=s0x2(1,-1)+fx
(tfi(ntraj),tfj(ntraj))* &
fx
(tfi(ntraj),tfj(ntraj))*tr
s0x2(2,-1)=s0x2(2,-1)+fy
(tfi(ntraj),tfj(ntraj))* &
fy
(tfi(ntraj),tfj(ntraj))*tr
s0x2(3,-1)=s0x2(3,-1)+ prof * &
prof *tr
ELSE
ntraj=ntraj-1
ENDIF
END DO
END DO
END DO
7050 FORMAT(4(1x,f7.3),1x,f12.3,3(1x,f6.3))
!
! 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
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 l=1,lmt
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_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
IF (key_nointerpolstats) THEN
IF (key_alltracers) THEN
t=tt(is,jt0s,ks,l)
s=ss(is,jt0s,ks,l)
r=rr(is,jt0s,ks,l)
ENDIF
ENDIF
CALL sub_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.(float(ndir(n))*v > 0)) THEN
IF (lmt == 1) THEN
n1=1+INT( SQRT(ABS(v)/t0) )
n2=1
ELSE
n1=1+INT( (ABS(v)/t0)**(1./3.) )
n2=n1
ENDIF
tr=ABS(v)/n1/n1/n2
DO i1=1,n1
DO k1=1,n1
DO l1=1,n2
ntraj=ntraj+1
IF (ntraj > nmax) THEN
WRITE(0,*)'too many particles, maximum is: ',nmax
STOP
ENDIF
tfi(ntraj)=i-1.-.5/n1+float(i1)/n1
tfj(ntraj)=jv0
tfk(ntraj)=k-.5/n1+float(k1)/n1
tfl(ntraj)=l-.5-.5/n2+float(l1)/n2
ttr(ntraj)=tr
IF (criter0
(tfi(ntraj),tfj(ntraj),-tfk(ntraj), &
tfl(ntraj),i,jt0,k,l)) THEN
trtot=trtot+tr
isn(-1)=isn(-1)+1
sn(-1)=sn(-1)+tr
IF (key_alltracers) THEN
IF (.NOT.key_nointerpolstats) THEN
t=zinter &
(tt,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
s=zinter &
(ss,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
IF (key_approximatesigma) THEN
r=zinter &
(rr,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
ELSE
r=sigma
(zsigma,s,t)
ENDIF
ENDIF
IF (t <= s0min(5,-1)) s0min(5,-1)=t
IF (s <= s0min(6,-1)) s0min(6,-1)=s
IF (r <= s0min(7,-1)) s0min(7,-1)=r
IF (t >= s0max(5,-1)) s0max(5,-1)=t
IF (s >= s0max(6,-1)) s0max(6,-1)=s
IF (r >= s0max(7,-1)) s0max(7,-1)=r
s0x(5,-1)=s0x(5,-1)+t*tr
s0x(6,-1)=s0x(6,-1)+s*tr
s0x(7,-1)=s0x(7,-1)+r*tr
s0x2(5,-1)=s0x2(5,-1)+t*t*tr
s0x2(6,-1)=s0x2(6,-1)+s*s*tr
s0x2(7,-1)=s0x2(7,-1)+r*r*tr
!!NG IF (key_outgridt) THEN
!!NG WRITE(50,7050) &
!!NG tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr, &
!!NG t,s,r
!!NG ELSE
WRITE(50,7050) &
tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr, &
t,s,r
!!NG ENDIF
ELSE
!!NG IF (key_outgridt) THEN
!!NG WRITE(50,7050) &
!!NG tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr
!!NG ELSE
WRITE(50,7050) &
tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr
!!NG ENDIF
ENDIF
IF (key_roms) THEN
prof=fz
(gi=tfi(ntraj), gj=tfj(ntraj), gk=-tfk(ntraj), il=l)
ELSE
prof=fz
(gk=-tfk(ntraj))
ENDIF
IF (fx(tfi(ntraj),tfj(ntraj)) <= s0min(1,-1)) &
s0min(1,-1)=fx
(tfi(ntraj),tfj(ntraj))
IF (fy(tfi(ntraj),tfj(ntraj)) <= s0min(2,-1)) &
s0min(2,-1)=fy
(tfi(ntraj),tfj(ntraj))
IF (-prof <= s0min(3,-1)) &
s0min(3,-1)=-prof
IF (fx(tfi(ntraj),tfj(ntraj)) >= s0max(1,-1)) &
s0max(1,-1)=fx
(tfi(ntraj),tfj(ntraj))
IF (fy(tfi(ntraj),tfj(ntraj)) >= s0max(2,-1)) &
s0max(2,-1)=fy
(tfi(ntraj),tfj(ntraj))
IF (-prof >= s0max(3,-1)) &
s0max(3,-1)=-prof
s0x(1,-1)=s0x(1,-1)+fx
(tfi(ntraj),tfj(ntraj))*tr
s0x(2,-1)=s0x(2,-1)+fy
(tfi(ntraj),tfj(ntraj))*tr
s0x(3,-1)=s0x(3,-1)- prof *tr
s0x2(1,-1)=s0x2(1,-1)+fx
(tfi(ntraj),tfj(ntraj))* &
fx
(tfi(ntraj),tfj(ntraj))*tr
s0x2(2,-1)=s0x2(2,-1)+fy
(tfi(ntraj),tfj(ntraj))* &
fy
(tfi(ntraj),tfj(ntraj))*tr
s0x2(3,-1)=s0x2(3,-1)+ prof * &
prof *tr
ELSE
ntraj=ntraj-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
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 l=1,lmt
DO j=jt1(n),jt2(n)
DO i=it1(n),it2(n)
CALL sub_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
IF (key_nointerpolstats) THEN
IF (key_alltracers) THEN
t=tt(is,js,kt0s,l)
s=ss(is,js,kt0s,l)
r=rr(is,js,kt0s,l)
ENDIF
ENDIF
CALL sub_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.(float(ndir(n))*w < 0)) THEN
IF (lmt == 1) THEN
n1=1+INT( SQRT(ABS(w)/t0) )
n2=1
ELSE
n1=1+INT( (ABS(w)/t0)**(1./3.) )
n2=n1
ENDIF
tr=ABS(w)/n1/n1/n2
DO j1=1,n1
DO i1=1,n1
DO l1=1,n2
ntraj=ntraj+1
IF (ntraj > nmax) THEN
WRITE(0,*)'too many particles, maximum is: ',nmax
STOP
ENDIF
tfi(ntraj)=i-1.-.5/n1+float(i1)/n1
tfj(ntraj)=j-1.-.5/n1+float(j1)/n1
tfk(ntraj)=kw0
tfl(ntraj)=l-.5-.5/n2+float(l1)/n2
ttr(ntraj)=tr
IF (criter0
(tfi(ntraj),tfj(ntraj),-tfk(ntraj), &
tfl(ntraj),i,j,kt0,l)) THEN
trtot=trtot+tr
isn(-1)=isn(-1)+1
sn(-1)=sn(-1)+tr
IF (key_alltracers) THEN
IF (.NOT.key_nointerpolstats) THEN
t=zinter &
(tt,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
s=zinter &
(ss,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
IF (key_approximatesigma) THEN
r=zinter &
(rr,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
ELSE
r=sigma
(zsigma,s,t)
ENDIF
ENDIF
IF (t <= s0min(5,-1)) s0min(5,-1)=t
IF (s <= s0min(6,-1)) s0min(6,-1)=s
IF (r <= s0min(7,-1)) s0min(7,-1)=r
IF (t >= s0max(5,-1)) s0max(5,-1)=t
IF (s >= s0max(6,-1)) s0max(6,-1)=s
IF (r >= s0max(7,-1)) s0max(7,-1)=r
s0x(5,-1)=s0x(5,-1)+t*tr
s0x(6,-1)=s0x(6,-1)+s*tr
s0x(7,-1)=s0x(7,-1)+r*tr
s0x2(5,-1)=s0x2(5,-1)+t*t*tr
s0x2(6,-1)=s0x2(6,-1)+s*s*tr
s0x2(7,-1)=s0x2(7,-1)+r*r*tr
!!NG
!!NG IF (key_outgridt) THEN
!!NG WRITE(50,7050) &
!!NG tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr, &
!!NG t,s,r
!!NG ELSE
WRITE(50,7050) &
tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr, &
t,s,r
!!NG ENDIF
ELSE
!!NG IF (key_outgridt) THEN
!!NG WRITE(50,7050) &
!!NG tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr
!!NG ELSE
WRITE(50,7050) &
tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr
!!NG ENDIF
ENDIF
IF (key_roms) THEN
prof=fz
(gi=tfi(ntraj), gj=tfj(ntraj), gk=-tfk(ntraj), il=l)
ELSE
prof=fz
(gk=-tfk(ntraj))
ENDIF
IF (fx(tfi(ntraj),tfj(ntraj)) <= s0min(1,-1)) &
s0min(1,-1)=fx
(tfi(ntraj),tfj(ntraj))
IF (fy(tfi(ntraj),tfj(ntraj)) <= s0min(2,-1)) &
s0min(2,-1)=fy
(tfi(ntraj),tfj(ntraj))
IF (-prof <= s0min(3,-1)) &
s0min(3,-1)=-prof
IF (fx(tfi(ntraj),tfj(ntraj)) >= s0max(1,-1)) &
s0max(1,-1)=fx
(tfi(ntraj),tfj(ntraj))
IF (fy(tfi(ntraj),tfj(ntraj)) >= s0max(2,-1)) &
s0max(2,-1)=fy
(tfi(ntraj),tfj(ntraj))
IF (-prof>= s0max(3,-1)) &
s0max(3,-1)=-prof
s0x(1,-1)=s0x(1,-1)+fx
(tfi(ntraj),tfj(ntraj))*tr
s0x(2,-1)=s0x(2,-1)+fy
(tfi(ntraj),tfj(ntraj))*tr
s0x(3,-1)=s0x(3,-1)- prof *tr
s0x2(1,-1)=s0x2(1,-1)+fx
(tfi(ntraj),tfj(ntraj))* &
fx
(tfi(ntraj),tfj(ntraj))*tr
s0x2(2,-1)=s0x2(2,-1)+fy
(tfi(ntraj),tfj(ntraj))* &
fy
(tfi(ntraj),tfj(ntraj))*tr
s0x2(3,-1)=s0x2(3,-1)+ prof * &
prof *tr
ELSE
ntraj=ntraj-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
END DO
!
! fin test kt1=kt2
ENDIF
!
! fin test num(n)
ENDIF
!
! fin du positionnement
END DO
ELSE
WRITE(0,*) ' posini: ibin =', ibin
ENDIF
!- Deallocate dynamic memory array
DEALLOCATE(it1)
DEALLOCATE(it2)
DEALLOCATE(jt1)
DEALLOCATE(jt2)
DEALLOCATE(kt1)
DEALLOCATE(kt2)
DEALLOCATE(num)
DEALLOCATE(ndir)
RETURN
END SUBROUTINE posini
!!***