!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! - 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/trajec
!! NAME
!! trajec (trajec.f90 - Main program - Fortran 90)
!!
!! USAGE
!!
!!
!! FUNCTION
!! --- Particle trajectories from the analytical computation of 3D
!! streamlines in a 3D velocity fiels.
!!
!! --- Qualitative or quantitative diagnostics
!!
!! --- Computations done over periods equal to the available sampling
!! for the velocity field, and over which the velocity is assumed
!! constant
!!
!! --- The program works with the OPA tensorial formalism
!!
!! --- input: simplified meshmask, 3D TRANSPORT field (velocity x
!! section), parameters (on various files)
!!
!! --- output: individual trajectories or quantitative outputs
!!
!! --- no-vectorized version
!!
!! --- some comments (yes !) in English
!!
!! --- copyright: Bruno Blanke (Bruno.Blanke@univ-brest.fr)
!! first version: Summer 1992
!! this version: March 2007
!! http://www.ifremer.fr/lpo/blanke/ARIANE
!!
!! 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
!!
!! (see mod_lun)
!! current input files:
!! -------------------
!! current output files:
!! ---------------------
!!
!! TODO
!!
!!
!! PORTABILITY
!! Machine-OS - Fortran90/95 compiler
!! * i686-pc-linux-gnu - ifort
!!
!! SEE ALSO
!!
!! SOURCE
!!=========================================================================
SUBROUTINE trajec 1,162
USE mod_precision
USE mod_cst
USE mod_namelist
USE mod_input_grid
USE mod_input_data
USE mod_output_data
USE mod_trajec_subs
USE mod_init_particules
USE mod_posin
USE mod_flags
USE mod_stats
USE mod_stati
USE mod_txt
USE mod_lun
USE mod_fx
USE mod_fy
USE mod_fz
USE mod_zinter
USE mod_sigma
USE mod_quant
USE mod_reducmem
USE mod_orca
USE mod_save_netcdf
!--------------!
! DECLARATIONS !
!--------------!
IMPLICIT NONE
INTEGER(kind=iprec) :: &
nsect2 = 0 , &
ind_traj = 0 , &
nfinold = 0 , & !
ibuoy = 0 , & !
iswap = 0 , & !
iperfl = 0 , & !
jperfl = 0 , & !
i_mod = 1 , & ! ntraj / 10
i, j,n , & !
is, js, ks , & !
nis, njs, nks , & !
i1, i2 , & !
i1s, i2s , & !
j1, j2 , & !
j1s, j2s , & !
k1, k2 , & !
k1s, k2s , & !
lu
INTEGER(kind=iprec), DIMENSION(:), ALLOCATABLE :: &
i0 , & !
j0 , & !
k0 , & !
l0 , & !
nfin
REAL(kind = qprec) :: ttt
REAL(kind=rprec) :: & !
subtime = -1._rprec , & !
tfic = 0._rprec , & !
tlap = 0._rprec , & !
tfin = 0._rprec , & !
ft , & !
rrr , & !
prof , & !
age , & !
tnext , & !
tswap , & !
x0, y0, z0 , & !
u, v, w , & !
du, dv, dw , & !
tx, ty, tz , & !
tfac , & !
t, s, r , & !
xt, yt, zt, st , & !
trtot1 , & !
r_lmt , & !
dumm
REAL(kind=rprec), DIMENSION(:), ALLOCATABLE :: & !
trueage, & !
agenew , & !
gi , & !
gj , & !
gk , & !
gl , & !
hi , & !
hj , & !
hk , & !
hl , & !
tf , & !
sf , & !
rf , & !
init_depth , & ! !! NG: 15_09_2008
final_depth !! NG: 15_09_2008
REAL(kind=rprec), DIMENSION(:,:,:), ALLOCATABLE :: & !
array_qual_traj
!-------------------------------------------!
!- DYNAMIC ALLOCATIONS AND INITIALIZATIONS -!
!-------------------------------------------!
CALL sub_posin_alloc
(nmax)
ALLOCATE(i0(1)) ; i0(:) = 0
CALL sub_memory
(size(i0),'i','i0','trajec')
ALLOCATE(j0(1)) ; j0(:) = 0
CALL sub_memory
(size(j0),'i','j0','trajec')
ALLOCATE(k0(1)) ; k0(:) = 0
CALL sub_memory
(size(k0),'i','k0','trajec')
ALLOCATE(l0(1)) ; l0(:) = 0
CALL sub_memory
(size(l0),'i','l0','trajec')
ALLOCATE(nfin(nmax)) ; nfin(:) = -1
CALL sub_memory
(size(nfin),'i','nfin','trajec')
ALLOCATE(trueage(nmax)); trueage(:)= 0._rprec
CALL sub_memory
(size(trueage),'r','trueage','trajec')
ALLOCATE(agenew(nmax)); agenew(:)= 0._rprec
CALL sub_memory
(size(agenew),'r','agenew','trajec')
tage(:) = 0._rprec ! define in mod_posin.f90
ALLOCATE(gi(1)) ; gi(:) = 0._rprec
CALL sub_memory
(size(gi),'r','gi','trajec')
ALLOCATE(gj(1)) ; gj(:) = 0._rprec
CALL sub_memory
(size(gj),'r','gj','trajec')
ALLOCATE(gk(1)) ; gk(:) = 0._rprec
CALL sub_memory
(size(gk),'r','gk','trajec')
ALLOCATE(gl(1)) ; gl(:) = 0._rprec
CALL sub_memory
(size(gl),'r','gl','trajec')
ALLOCATE(hi(nmax)) ; hi(:) = 0._rprec
CALL sub_memory
(size(hi),'r','hi','trajec')
ALLOCATE(hj(nmax)) ; hj(:) = 0._rprec
CALL sub_memory
(size(hj),'r','hj','trajec')
ALLOCATE(hk(nmax)) ; hk(:) = 0._rprec
CALL sub_memory
(size(hk),'r','hk','trajec')
ALLOCATE(hl(nmax)) ; hl(:) = 0._rprec
CALL sub_memory
(size(hl),'r','hl','trajec')
ALLOCATE(tf(nmax)) ; tf(:) = 0._rprec
CALL sub_memory
(size(tf),'r','tf','trajec')
ALLOCATE(sf(nmax)) ; sf(:) = 0._rprec
CALL sub_memory
(size(sf),'r','sf','trajec')
ALLOCATE(rf(nmax)) ; rf(:) = 0._rprec
CALL sub_memory
(size(rf),'r','rf','trajec')
!! NG: 15_09_2008
ALLOCATE(init_depth(nmax)) ; init_depth(:) = 0._rprec
CALL sub_memory
(size(init_depth),'i','init_depth','trajec')
ALLOCATE(final_depth(nmax)) ; final_depth(:) = 0._rprec
CALL sub_memory
(size(final_depth),'i','final_depth','trajec')
!======================!
!- Read region limits -!
!======================!
CALL sub_reducmem_read_reg_lim
()
IF (key_ascii_outputs) THEN
!==========================!
! Open ASCII output files -!
!==========================!
CALL sub_output_data_open_files
()
ENDIF
!-------------------!
!- READ INPUT MESH -!
!-------------------!---------------------------------------------------!
!-- Allocate and read coordinates and scale factors
!-- xx_tt, xx_uu, xx_vv, xx_ff
!-- yy_tt, yy_uu, yy_vv, yy_ff
!-- zz_tt, zz_ww
!-- e2u
!-- e1v
!-- e1t, e2t, e3t
!-- tmask
!-----------------------------------------------------------------------!
CALL sub_input_grid
()
!==========================
!- QUALITATIVE VARIABLES -!
!=================================================================
! sampling time (in seconds) for the available transport field
!=================================================================
tfic = tunit * REAL(ntfic, kind=rprec)
!=================================================================
! sampling time (in seconds) between two successive output positions
!=================================================================
tlap = ABS(delta_t) * REAL(frequency, kind=rprec)
!=================================================================
! maximum integration time (im seconds) for individual trajectories
!=================================================================
tfin = tlap * REAL(ABS(nb_output), kind=rprec)
!==============================================
! mask is written (for graphical trajectories)
!==============================================
IF ((TRIM(mode) == 'qualitative').AND.(mask)) THEN
DO j = 1, jmt
DO i = 1, imt
IF (tmask(i,j,1,1) == 0._rprec) THEN
IF (key_alltracers) THEN
WRITE(lun_traj,7055)0,xx_tt(i,j,1,1),yy_tt(i,j,1,1),0.,0.,0.,0.,0.
ELSE
WRITE(lun_traj,7055)0,xx_tt(i,j,1,1),yy_tt(i,j,1,1),0.,0.
ENDIF
ENDIF
END DO
END DO
ENDIF
!!NG: BEGIN DIRECT DATA ACCESS
!-------------------!
!- READ INPUT DATA -!
!-------------------!---------------------------------------------------!
!-- Allocate and read (or compute) uu, vv, ww, [tt], [ss] and [rr].
!-- and store them in the memory (RAM) of the computer.
!-----------------------------------------------------------------------!
CALL sub_input_data
()
!=================================================================
! opposite transports (-1) if negative time step
!=================================================================
IF (TRIM(forback) == 'backward') THEN
delta_t = - delta_t
uu(:,:,:,:) = -uu(:,:,:,:)
vv(:,:,:,:) = -vv(:,:,:,:)
ww(:,:,:,:) = -ww(:,:,:,:)
ENDIF
!!NG: END DIRECT DATA ACCESS
!=================================================================
! we initialize the *new* particle position
!=================================================================
hi(:) = 0._rprec
hj(:) = 0._rprec
hk(:) = 0._rprec
hl(:) = 0._rprec
!=======================================!
!- Read or Compute particule positions -!
!=======================================!
CALL sub_init_particules_positions
()
!=======================!
!- NetCDF data outputs -!
!=======================!
!---------------------------!
!- Open Netcdf output file -!
!---------------------------!
CALL sub_save_netcdf_data_init
()
!------------------------------------------------!
!- Save Initial Positions in netcdf output file -!
!------------------------------------------------!
CALL sub_save_netcdf_init_pos
()
!- Deallocate init_temp, init_salt, init_dens (mod_posin) -!
CALL sub_posin_dealloc_init_tracers
()
IF (TRIM(mode) == 'quantitative') THEN
CALL sub_quant_alloc
(imt, jmt, kmt)
CALL sub_quant_init
() ! initialize a few arrays
!-------------------------------------------------------------------------
! we do not yet know whether a final criterion (criter1) will be satisfied
!-------------------------------------------------------------------------
icrit1 = 0
secname(nsect+1) = 'Criter1'
ELSEIF (TRIM(mode) == 'qualitative') THEN
ALLOCATE(array_qual_traj(ntraj,nb_output+1,nstat))
CALL sub_memory
(size(array_qual_traj),'r','array_qual_traj','trajec')
array_qual_traj(:,:,:)=mask_value
ELSE
STOP
ENDIF
!==================================================================!
!==================================================================!
! General loop on the number of trajectories. !
!==================================================================!
!==================================================================!
WRITE(lun_standard,*)''
WRITE(lun_standard,*)'============================================='
WRITE(lun_standard,1122) ntraj
WRITE(lun_standard,*)'============================================='
1122 FORMAT (' = ARIANE IS COMPUTING ',I0.10,' PARTICLES =')
IF ( ntraj < 10 ) THEN
i_mod = 1
ELSEIF ( ntraj < 100) THEN
i_mod = 10
ELSEIF ( ntraj < 1000) THEN
i_mod = 100
ELSEIF ( ntraj < 10000) THEN
i_mod = 1000
ELSEIF ( ntraj < 100000) THEN
i_mod = 10000
ELSEIF ( ntraj < 1000000) THEN
i_mod = 100000
ELSEIF ( ntraj < 10000000) THEN
i_mod = 1000000
ELSEIF ( ntraj < 100000000) THEN
i_mod = 10000000
ELSEIF ( ntraj < 1000000000) THEN
i_mod = 100000000
ELSE
i_mod = 1000000000
ENDIF
!=====================================================================!
! DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 !
!=====================================================================!
DL0:DO n = 1, ntraj
!------------!
!- Comments -!
!------------!
!!NG IF(n == ntraj) THEN
!!NG WRITE(lun_standard,'(I0,A)', advance='yes') n, '.'
!!NG ELSEIF ((n == 1).OR.(MOD(n, i_mod) == 0)) THEN
!!NG WRITE(lun_standard,'(I0,3A)', advance='no') n, ' - '
!!NG ENDIF
IF ((n == 1).OR.(MOD(n, i_mod) == 0).OR.(n == ntraj)) THEN
WRITE(lun_standard,'(10X,I0)') n
ENDIF
!====================================================================
! Few initializations (transport arrays) [QUANT]
!====================================================================
IF (TRIM(mode) == 'quantitative') THEN
nfinold = 0
IF (.NOT.key_eco) THEN
!----------------------!
!- Take a lot of time -!
!----------------------!
CALL sub_quant_inittmp
() ! Inititalize arrays
ENDIF
ENDIF
!====================================================================
! time index must be within [0.5, lmt+.5[; otherwise time translation
!NG time index must be within ]0.5, lmt+.5]; otherwise time translation
!====================================================================
!NG !!! To work "lmt" must be greater than zero !!! (else infinit loop)
!NG !!! "lmt" must be tested before to use it. !!!
!NG !!! If "lmt" is tested, GOTO could be removed. !!!
100 IF (tfl(n) < 0.5_rprec) THEN
tfl(n) = tfl(n) + REAL(lmt, kind=rprec)
GOTO 100
ENDIF
IF (tfl(n) >= (REAL(lmt,kind=rprec)+0.5_rprec)) THEN
tfl(n) = tfl(n) - REAL(lmt, kind=rprec)
GOTO 100
ENDIF
!====================================================================
! the program is able to deal with constant-depth particles
! (qualitative experiments only) [QUALI]
!====================================================================
ibuoy = 0
IF (tfk(n) <= 0._rprec) THEN
IF (TRIM(mode) == 'quantitative') THEN
WRITE(lun_error,*)''
WRITE(lun_error,*)'--- ERROR ---'
WRITE(lun_error,*)'Isobaric particles allowed in QUALITATIVE experiments ONLY '
WRITE(lun_error,*)'Particule number n = ', n
WRITE(lun_error,*)' tfk(n) = ', tfk(n)
WRITE(lun_error,*)'--- STOP ARIANE ---'
STOP
ENDIF
ibuoy = 1
tfk(n) = -tfk(n)
ENDIF
!====================================================================
! time (tfl(n)) and position (VELOCITY grid: fi, tfj(n), tfk(n))
! indices used to detect the relevant embedding temperature
! gridcell (i, j, k, l).
!====================================================================
i0(1) = INT(tfi(n)) + 1
j0(1) = INT(tfj(n)) + 1
k0(1) = INT(tfk(n))
l0(1) = NINT(tfl(n))
!===========================!
!- Periodicity in OPA-ORCA -!
!===========================!
IF (.NOT.(key_roms.OR.key_symphonie)) THEN
CALL sub_orca_north_pole
(i0(1), j0(1))
CALL sub_orca_east_west_periodic
(i0(1))
ENDIF
!====================================================================
! test of the initial position
!====================================================================
CALL sub_reducmem_shift_or_not_ind
(i0(1),j0(1),k0(1),is,js,ks)
IF (tmask(is,js,ks,1) <= .5_rprec) THEN
WRITE(lun_error,7000)' false start', &
n,i0(1),j0(1),k0(1),l0(1),tfi(n),tfj(n),tfk(n),tfl(n)
CYCLE DL0
ENDIF
!====================================================================
! tfi(n), tfj(n), tfk(n) and tfl(n) refer to the initial position
! gi, gj, gk and gl refer to the current position
! BEWARE: we use a reverse vertical axis for easier computations
!====================================================================
gi(1) = tfi(n)
gj(1) = tfj(n)
gk(1) = -tfk(n)
gl(1) = tfl(n)
!====================================================================
! age initialization
!====================================================================
age = tage(n)
!====================================================================
! next output instant (for positions in qualitative experiments)
!====================================================================
tnext = tlap
!====================================================================
! next time swap instant (between two samples of velocity field)
!====================================================================
!NG "l" is an integer => REAL(l, kind=rprec)
IF (TRIM(forback) == 'forward') THEN
tswap = tfic * ABS( tfl(n)-(l0(1)+.5_rprec) ) + tage(n)
ELSE
tswap = tfic * ABS( tfl(n)-(l0(1)-.5_rprec) ) + tage(n)
ENDIF
!!BBL - march 2007
!!NG: option to follow an experience interrupted. Age is stored
!!NG: in the final position file.
!!NG: Dont USE this option in backward mode !!!
IF (key_read_age) THEN
dumm = tswap / tfic
tswap = NINT(tswap-NINT(dumm)*tfic)+NINT(dumm)*tfic
ENDIF
!!NG: to be compatible to sequential version - 09/01/2007
!!NG: we had also the key_nointerpolstats
!===================================================================
! Compute diagnostics for initial tracers values
!===================================================================
CALL sub_reducmem_shift_or_not_ind
(i0(1),j0(1),k0(1),is,js,ks)
IF (key_alltracers) THEN
IF (key_nointerpolstats) THEN
tf(n) = tt(is,js,ks,l0(1))
sf(n) = ss(is,js,ks,l0(1))
rf(n) = rr(is,js,ks,l0(1))
ELSE
tf(n)=zinter
(tt,tfi(n),tfj(n),-tfk(n),tfl(n))
sf(n)=zinter
(ss,tfi(n),tfj(n),-tfk(n),tfl(n))
IF (key_approximatesigma) THEN
rrr=zinter
(rr,tfi(n),tfj(n),-tfk(n),tfl(n))
ELSE
rrr=sigma
(zsigma,sf(n),tf(n))
ENDIF
ENDIF
ENDIF
!!NG : 09/01/2007
!====================================================================
! the transport on the initial section is stored in the relevant
! transport arrays (quantitative experiment) [QUANT]
!====================================================================
IF (TRIM(mode) == 'quantitative') THEN
IF (key_roms.OR.key_symphonie) THEN
prof=fz
(gi=tfi(n), gj=tfj(n), gk=-tfk(n), il=l0(1))
ELSE
prof=fz
(gk=-tfk(n))
ENDIF
!! NG: 15_09_2008
!! NG: Here we stored the depth of the particle at the intial position.
!! NG: Because with ROMS or symphonie the zz_ww in mod_fz is different at
!! NG: time step.
init_depth(n)=prof
CALL sub_reducmem_shift_or_not_ind
(NINT(tfi(n)),NINT(tfj(n)),NINT(tfk(n)),nis,njs,nks)
IF ( (ABS(ANINT(tfi(n))-tfi(n)) <= 1.e-5_rprec) .AND. (uu(nis,js,ks,l0(1)) > 0._rprec) ) THEN
IF (key_alltracers) THEN
CALL sub_quant_initside_u
(NINT(tfi(n)), j0(1), k0(1), l0(1), &
prof, ttr(n), tf(n), sf(n), rrr)
ELSE
CALL sub_quant_initside_u
(NINT(tfi(n)), j0(1), k0(1), l0(1), &
prof, ttr(n))
ENDIF
ENDIF
IF ((ABS(ANINT(tfj(n))-tfj(n)) <= 1.e-5_rprec).AND.(vv(is,njs,ks,l0(1)) > 0._rprec)) THEN
IF (key_alltracers) THEN
CALL sub_quant_initside_v
(i0(1), NINT(tfj(n)), k0(1), l0(1), &
prof, ttr(n), tf(n), sf(n), rrr)
ELSE
CALL sub_quant_initside_v
(i0(1), NINT(tfj(n)), k0(1), l0(1), &
prof, ttr(n))
ENDIF
ENDIF
IF ((ABS(ANINT(tfk(n))-tfk(n)) <= 1.e-5_rprec).AND.(ww(is,js,nks,l0(1)) < 0._rprec)) THEN
CALL sub_quant_initside_w
(i0(1), j0(1), NINT(tfk(n)), l0(1), ttr(n))
ENDIF
!====================================================================
! the initial position is written on "traj.qt" [QUALI]
!====================================================================
ELSEIF (TRIM(mode) == 'qualitative') THEN
x0 = fx
(gi(1),gj(1))
y0 = fy
(gi(1),gj(1))
IF (key_roms.OR.key_symphonie) THEN
z0 = fz
(gi=gi(1), gj=gj(1), gk=gk(1), il=l0(1))
ELSE
z0 = fz
(gk=gk(1))
ENDIF
ind_traj = 1
array_qual_traj(n,ind_traj,1) = x0
array_qual_traj(n,ind_traj,2) = y0
array_qual_traj(n,ind_traj,3) = z0
array_qual_traj(n,ind_traj,4) = 0._rprec
IF (key_alltracers) THEN
!! WRITE(lun_traj,7055) n, x0, y0, z0, 0., tf(n), sf(n), rrr
array_qual_traj(n,ind_traj,5) = tf(n)
array_qual_traj(n,ind_traj,6) = sf(n)
array_qual_traj(n,ind_traj,7) = rrr
!! ELSE
!! WRITE(lun_traj,7055) n, x0, y0, z0, 0.
ENDIF
ELSE
STOP
ENDIF
iswap = 0 !!NG-bug: -- BUG FIXED by this initialization --
!! In qualitatif mode it was not possible to compute
!! more than one particule.
!====================================================================
! entry point (new mesh)
!====================================================================
!=====================================================================!
! DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 !
!=====================================================================!
DL1: DO WHILE( (TRIM(mode) == 'quantitative').OR.((age < tfin).AND.(iswap /= 2)) )
!====================================================================
! linear interpolation for the local velocity
!====================================================================
CALL sub_reducmem_shift_or_not_ind
(i0(1),j0(1),k0(1),is,js,ks)
u = uu(is-1,js ,ks,l0(1)) + &
(gi(1) - REAL(i0(1)-1, kind=rprec)) * (uu(is,js,ks ,l0(1)) - uu(is-1,js ,ks,l0(1)))
v = vv(is ,js-1,ks,l0(1)) + &
(gj(1) - REAL(j0(1)-1, kind=rprec)) * (vv(is,js,ks ,l0(1)) - vv(is ,js-1,ks,l0(1)))
w = ww(is ,js ,ks,l0(1)) - &
(gk(1) + REAL(k0(1) , kind=rprec)) * (ww(is,js,ks+1,l0(1)) - ww(is ,js ,ks,l0(1)))
!====================================================================
! entry and exit indices (depending on the sign of the velocity)
!====================================================================
i1 = i0(1) - 1 !
i2 = i0(1) !
IF (u < 0._rprec) THEN ! i1 = i - (1 + SIGN(1,u)) * 1/2
i2 = i0(1) - 1 ! i2 = i - (1 - SIGN(1,u)) * 1/2
i1 = i0(1) ! ????? u=0. ?????
ENDIF !
!
j1 = j0(1) - 1 !
j2 = j0(1) !
IF (v < 0._rprec) THEN ! j1 = j - (1 + SIGN(1,v)) * 1/2
j2 = j0(1) - 1 ! j2 = j - (1 - SIGN(1,v)) * 1/2
j1 = j0(1) ! ????? v=0. ?????
ENDIF !
!
k1 = k0(1) + 1 !
k2 = k0(1) !
IF (w < 0._rprec) THEN ! k1 = k + (1 + SIGN(1,w)) * 1/2
k2 = k0(1) + 1 ! k2 = k + (1 - SIGN(1,w)) * 1/2
k1 = k0(1) ! ????? w=0. ?????
ENDIF !
!====================================================================
! times to get to each edge of the (i, j, k) grid cell:
! - infinite, if Uexit * Ulocal < 0
! - linear relationship, if Uexit = Ulocal
! - logarithmi! formulation, in all other cases
! BEWARE: use of log(u)-log(v) instead of log (u/v),
! better for numerical results
!====================================================================
CALL sub_reducmem_shift_or_not_ind
(i0(1) ,j0(1) ,k0(1) ,is ,js ,ks )
CALL sub_reducmem_shift_or_not_ind
(i1,j1,k1,i1s,j1s,k1s)
CALL sub_reducmem_shift_or_not_ind
(i2,j2,k2,i2s,j2s,k2s)
du = ( uu(i2s,js ,ks ,l0(1)) - uu(i1s,js ,ks ,l0(1)) ) * REAL(i2-i1, kind=rprec)
dv = ( vv(is ,j2s,ks ,l0(1)) - vv(is ,j1s,ks ,l0(1)) ) * REAL(j2-j1, kind=rprec)
dw = ( ww(is ,js ,k2s,l0(1)) - ww(is ,js ,k1s,l0(1)) ) * REAL(k1-k2, kind=rprec)
!===============
! Scale factors
!===============
IF (key_roms.OR.key_symphonie) THEN
!- ROMS case => E3T is dependent of time -!
tfac = e3t(is,js,ks,l0(1)) * e1t(is,js,1,1) * e2t(is,js,1,1)
ELSE
!- OPA case => E3T is not dependent of time -!
IF (key_partialsteps) THEN
tfac = e3t(is,js,ks,1) * e1t(is,js,1,1) * e2t(is,js,1,1)
ELSE
tfac = e3t(1,1,ks,1) * e1t(is,js,1,1) * e2t(is,js,1,1)
ENDIF
ENDIF
!=============================
! Compute time to exit a cell
!=============================
IF (( u * uu(i2s,js,ks,l0(1))) <= 0._rprec) THEN
tx = 1.e35_rprec
ELSEIF (du == 0._rprec) THEN
tx = ( REAL(i2, kind=rprec) - gi(1) ) / u
ELSE
tx = ( LOG(ABS(uu(i2s,js,ks,l0(1)))) - LOG(ABS(u)) ) / du
ENDIF
IF ((v * vv(is,j2s,ks,l0(1))) <= 0._rprec) THEN
ty = 1.e35_rprec
ELSEIF (dv == 0._rprec) THEN
ty = ( REAL(j2, kind=rprec) - gj(1) ) / v
ELSE
ty = ( LOG(ABS(vv(is,j2s,ks,l0(1)))) - LOG(ABS(v)) ) / dv
ENDIF
!!NG: pour BBL
!!OLD: IF ((ibuoy /= 0).OR.((w*ww(is,js,k2s,l0(1))) <= 0.)) THEN
IF ((key_2dquant).OR.(ibuoy /= 0).OR.((w*ww(is,js,k2s,l0(1))) <= 0._rprec)) THEN
tz = 1.e35_rprec
ELSEIF (dw == 0._rprec) THEN
tz = -( REAL(k2, kind=rprec) + gk(1) ) / w
ELSE
tz = ( LOG(ABS(ww(is,js,k2s,l0(1)))) - LOG(ABS(w)) ) / dw
ENDIF
!====================================================================
! the true exit time is given by the shortest times (of tx, ty, tz)
! entry and exit indices are then updated
!====================================================================
ttt = MIN(tx,ty,tz)
!!NG: pour BBL
IF (.NOT.key_2dquant) THEN
IF ( (ttt >= 1.e34_qprec) .AND. (ibuoy == 0) ) THEN
WRITE(lun_error,7000)' dead end',n,i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
CYCLE DL0
ENDIF
ENDIF
IF (ttt < 0._qprec) THEN
IF (ABS(tfac*ttt) > 1._rprec) THEN
WRITE(lun_error,7000)' negative time',n,i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
CYCLE DL0
ENDIF
ttt = 0._qprec
ENDIF
!====================================================================
! do we need to update the velocity field WITHIN the gridcell ?
!====================================================================
iswap = 0
!====================================================================
! update needed
!====================================================================
IF ((lmt > 1).AND.((age+ttt*tfac) > tswap)) THEN
ttt = ( tswap - age ) / tfac
iswap = 1
ENDIF
!====================================================================
! end of the integration [Qualitative]
!====================================================================
IF ((TRIM(mode) == 'qualitative').AND.((age+ttt*tfac) > tfin)) THEN
ttt = ( tfin - age ) / tfac
iswap = 2
ENDIF
IF (TRIM(forback) == 'forward') THEN
hl(n) = gl(1) + ttt * tfac / tfic
ELSE
hl(n) = gl(1) - ttt * tfac / tfic
ENDIF
!====================================================================
! case 1: we do not reach the side of the cell (for a given direction)
! diagnostic of the final position
! - linear formulation, if Uexit = Ulocal
! - exponential formulation, in other cases
!====================================================================
IF (tx > ttt) THEN
CALL sub_dont_reachside
(hi(n), du, gi(1), u, ttt, i1, i2, np=n)
ENDIF
IF (ty > ttt) THEN
CALL sub_dont_reachside
(hj(n), dv, gj(1), v, ttt, j1, j2)
ENDIF
!!NG: BBL
IF ((ibuoy == 0).AND.(.NOT.key_2dquant)) THEN
IF (tz > ttt) THEN
CALL sub_dont_reachside
(hk(n), dw, gk(1), w, ttt, -k1, -k2)
ENDIF
ELSEIF((ibuoy /= 0).OR.(key_2dquant)) THEN
hk(n) = gk(1) ! IF (ibuoy /= 0) hk(n) = gk(1)
ENDIF
!!NG IF (ibuoy /= 0) hk(n) = gk(1)
!====================================================================
! case 2: we DO reach the side of the cell (for a given direction)
! no need for time update within the gridcell
! - we record the exiting transport
! - we update the positions index (hi, hj(n), hk(n))
!====================================================================
IF (tx <= ttt) THEN
hi(n)=REAL(i2,kind=rprec)
IF (TRIM(mode) == 'quantitative') &
CALL sub_quant_reachside_u
(i2, j0(1), k0(1), l0(1), hi(n), hj(n), hk(n), hl(n), ttr(n))
IF (i2 > i1) THEN
i1 = i2
i2 = i2 + 1
ELSE
i1 = i2
i2 = i2 - 1
ENDIF
ENDIF
IF (ty <= ttt) THEN
hj(n) = REAL(j2,kind=rprec)
IF (TRIM(mode) == 'quantitative') THEN
CALL sub_quant_reachside_v
( &
i0(1), j2, k0(1), l0(1), &
hi(n), hj(n), hk(n), hl(n), &
ttr(n))
ENDIF
IF (j2 > j1) THEN
j1 = j2
j2 = j2 + 1
ELSE
j1 = j2
j2 = j2 - 1
ENDIF
ENDIF
IF (tz <= ttt) THEN
IF (TRIM(mode) == 'quantitative') THEN
CALL sub_quant_reachside_w
(i0(1), j0(1), k2, l0(1), ttr(n))
ENDIF
!!NG: BBL
IF (.NOT.key_2dquant) THEN
IF (ibuoy == 0) THEN
hk(n) = -REAL(k2,kind=rprec)
IF (k2 > k1) THEN
k1 = k2
k2 = k2 + 1
ELSE
k1 = k2
k2 = k2 - 1
ENDIF
ENDIF
ENDIF
ENDIF
!====================================================================
! update of the T-grid index of the new embedding gridcell
!====================================================================
i0(1) = max0(i1, i2)
j0(1) = max0(j1, j2)
k0(1) = min0(k1, k2)
!!NG-bug: IF k == 0 => free-surface bug
!===========================!
!- Periodicity in OPA-ORCA -!
!===========================!
IF (.NOT.(key_roms.OR.key_symphonie)) THEN
CALL sub_orca_north_pole
(i1 , j1 , j0(1)+1)
CALL sub_orca_north_pole
(i2 , j2 , j0(1)+1)
CALL sub_orca_north_pole
(hi(n), hj(n), j0(1)+1)
CALL sub_orca_north_pole
(i0(1), j0(1), j0(1)+1, jperfl)
CALL sub_orca_east_west_periodic
(i1 , min0(i1,i2)+1, i0(1))
CALL sub_orca_east_west_periodic
(i2 , min0(i1,i2)+1, i0(1))
CALL sub_orca_east_west_periodic
(hi(n), min0(i1,i2)+1, i0(1))
CALL sub_orca_east_west_periodic
(i0(1), min0(i1,i2)+1, i0(1), iperfl)
ENDIF
!====================================================================
! age update
!====================================================================
agenew(n) = age + ttt * tfac
IF (iswap == 1) THEN
hl(n) = NINT(hl(n)+.5_rprec) - .5_rprec
agenew(n) = tswap
ENDIF
!====================================================================
! we stop the integration if we reach a limit of the domain (it may
! happen in QUALITATIVE experiments) [QUALI]
!====================================================================
IF (TRIM(mode) == 'qualitative') THEN
IF (k0(1) == 0) THEN
WRITE(lun_error,7000)'out of zdomain',n,i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
CYCLE DL0
ENDIF
IF ((j0(1) == dims_reg(2,1)).OR.(j0(1) == dims_reg(2,2))) THEN
WRITE(lun_error,7000)'out of ydomain',n,i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
CYCLE DL0
ENDIF
IF (.NOT.key_periodic) THEN
IF ((i0(1) == dims_reg(1,1)).OR.(i0(1) == dims_reg(1,2))) THEN
WRITE(lun_error,7000)'out of xdomain',n,i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
CYCLE DL0
ENDIF
ENDIF
ENDIF
!====================================================================
! quantitative analysis
!
! in QUANTITATIVE experiments, the 2D projections of the transport
! field associated with 1 particle is summed depending on the final
! section [QUANT]
!====================================================================
IF (TRIM(mode) == 'quantitative') THEN
!------------------------------------------------------------------
! mtfin informs about the sections used in quantitative experiments
!------------------------------------------------------------------
nfin(n) = mtfin(i0(1),j0(1),k0(1))
!------------------------------------------------
! test of a *transparent* hydrological criterion
!------------------------------------------------
IF (criter2
(tfi(n),tfj(n),tfk(n),tfl(n), &
hi(n),hj(n),hk(n),hl(n),&
i0(1),j0(1),k0(1),l0(1))) THEN
trueage(n) = agenew(n) / tcyc
IF (key_nointerpolstats) THEN
IF (key_alltracers) THEN
CALL sub_reducmem_shift_or_not_ind
(i0(1),j0(1),k0(1),is,js,ks)
t=tt(is,js,ks,l0(1))
s=ss(is,js,ks,l0(1))
r=rr(is,js,ks,l0(1))
ENDIF
ELSE
IF (key_alltracers) THEN
t=zinter
(tt,hi(n),hj(n),hk(n),hl(n))
s=zinter
(ss,hi(n),hj(n),hk(n),hl(n))
IF (key_approximatesigma) THEN
r=zinter
(rr,hi(n),hj(n),hk(n),hl(n))
ELSE
r=sigma
(zsigma,s,t)
ENDIF
ENDIF
ENDIF
IF (key_alltracers) THEN
WRITE(lun_trans,1002)n,hi(n),hj(n),-hk(n),hl(n),ttr(n),trueage(n),t,s,r
ELSE
WRITE(lun_trans,1002)n,hi(n),hj(n),-hk(n),hl(n),ttr(n),trueage(n)
ENDIF
1002 FORMAT(i0,3(1x,f0.3),1x,f0.1,1x,f0.3,1x,f0.2, 3(1x,f0.3))
ENDIF
!===============================================================
! a final criterion may define a "final section" (criter1)
!===============================================================
IF ((nfin(n) > 0).OR.(criter1(tfi(n),tfj(n),tfk(n),tfl(n), &
hi(n),hj(n),hk(n),hl(n), &
i0(1),j0(1),k0(1),l0(1)))) THEN
EXIT DL1
!============================================================
! are we on a transparent section ?
!============================================================
ELSEIF (nfin(n) < 0) THEN
IF (nfinold /= nfin(n)) THEN
nfinold = nfin(n)
lu = 90 - nfin(n)
idiru = 'E'
idirv = 'N'
idirw = '^'
IF (u < 0) idiru = 'W'
IF (v < 0) idirv = 'S'
IF (w < 0) idirw = 'v'
trueage(n) = agenew(n) / tcyc
IF (key_nointerpolstats) THEN
IF (key_alltracers) THEN
CALL sub_reducmem_shift_or_not_ind
(i0(1),j0(1),k0(1),is,js,ks)
t=tt(is,js,ks,l0(1))
s=ss(is,js,ks,l0(1))
r=rr(is,js,ks,l0(1))
ENDIF
ELSE
IF (key_alltracers) THEN
t=zinter
(tt,hi(n),hj(n),hk(n),hl(n))
s=zinter
(ss,hi(n),hj(n),hk(n),hl(n))
IF (key_approximatesigma) THEN
r=zinter
(rr,hi(n),hj(n),hk(n),hl(n))
ELSE
r=sigma
(zsigma,s,t)
ENDIF
ENDIF
ENDIF
IF (key_alltracers) THEN
WRITE(lu,1001)n,hi(n),hj(n),-hk(n),hl(n),ttr(n),trueage(n),idiru, &
idirv,idirw, t,s,r
ELSE
WRITE(lu,1001)n,hi(n),hj(n),-hk(n),hl(n),ttr(n),trueage(n),idiru, &
idirv,idirw
ENDIF
1001 FORMAT(i0,3(1x,f0.3),1x,f0.1,1x,f0.3,1x,f0.2,3(1x,a2), &
3(1x,f0.3))
ENDIF
!============================================================
! IF (nfin(n) == 0)
!============================================================
ELSE
nfinold=0
ENDIF
ENDIF !=> IF (TRIM(mode) == 'quantitative') THEN
!------------------------------
! end of quantitative analysis
!------------------------------
!!NG: Here a Bug in Roms mode when iswap /= 0, because l was modified
!!NG: and in qualitative mode l is used to compute depth (z0)
!!NG: in the fz module....
!!NG: The "time index update of the velocity field if needed" has been
!!NG: moved after the "qualitative" test...
!====================================================================
! THIS SHOULD NEVER HAPPEN !
!====================================================================
CALL sub_reducmem_shift_or_not_ind
(i0(1),j0(1),k0(1),is,js,ks)
IF (tmask(is,js,ks,1) == 0._rprec) THEN
WRITE(lun_error,7000)' coast crash',n,i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
WRITE(lun_coast_crash,*)' '
WRITE(lun_coast_crash,*)'coast crash for particle # ',n
WRITE(lun_coast_crash,*)'i j k: ',i0(1),' ',j0(1),' ',k0(1)
WRITE(lun_coast_crash,*)'is js ks: ',is,' ',js,' ',ks
WRITE(lun_coast_crash,*)'gi gj gk:',gi(1),' ',gj(1),' ',gk(1)
WRITE(lun_coast_crash,*)'hi(n) hj(n) hk(n):',hi(n),' ',hj(n),' ',hk(n)
WRITE(lun_coast_crash,*)'u: ',uu(is-1,js ,ks ,l0(1)),uu(is,js,ks,l0(1))
WRITE(lun_coast_crash,*)'v: ',vv(is ,js-1,ks ,l0(1)),vv(is,js,ks,l0(1))
WRITE(lun_coast_crash,*)'w: ',ww(is ,js ,ks+1,l0(1)),ww(is,js,ks,l0(1))
WRITE(lun_coast_crash,*)'tx ty tz ttt: ',tx,ty,tz,ttt
STOP
ENDIF
!====================================================================
! do we need an output of the position (qualitative experiment) ?
! (linear interpolation)
! - age: previous age
! - agenew: current age
! - tnext: output time [QUALI]
!====================================================================
IF (TRIM(mode) == 'qualitative') THEN
DL2:DO WHILE (agenew(n) >= tnext)
ft = (tnext-age) / (agenew(n)-age)
xt = gi(1) + (hi(n)-gi(1)) * ft
yt = gj(1) + (hj(n)-gj(1)) * ft
zt = gk(1) + (hk(n)-gk(1)) * ft
IF (TRIM(forback) == 'forward') THEN
st = gl(1) + (tnext-age) / tfic
ELSE
st = gl(1) - (tnext-age) / tfic
ENDIF
IF (key_periodic) THEN
IF (iperfl == 1) THEN
xt = gi(1) + (hi(n)-REAL(imt-2,kind=rprec)-gi(1)) * ft
ENDIF
IF (iperfl == 2) THEN
xt = gi(1) + (hi(n)+REAL(imt-2,kind=rprec)-gi(1)) * ft
ENDIF
ENDIF
IF (key_jfold) THEN
!-----------------------------------
! extrapolation not yet implemented
!-----------------------------------
IF (jperfl == 1) THEN
xt = gi(1)
yt = gj(1)
zt = gk(1)
ENDIF
ENDIF
!-------------------------
! position output on file
!-------------------------
x0 = fx
(xt,yt)
y0 = fy
(xt,yt)
IF (key_roms.OR.key_symphonie) THEN
z0 = fz
(gi=xt, gj=yt, gk=zt, il=l0(1))
ELSE
z0 = fz
(gk=zt)
ENDIF
IF (tnext <= tfin) THEN
ind_traj = ind_traj + 1
array_qual_traj(n,ind_traj,1) = x0
array_qual_traj(n,ind_traj,2) = y0
array_qual_traj(n,ind_traj,3) = z0
array_qual_traj(n,ind_traj,4) = tnext/tcyc
IF (key_alltracers) THEN
tf(n) = zinter
(tt,xt,yt,zt,st)
sf(n) = zinter
(ss,xt,yt,zt,st)
IF (key_approximatesigma) THEN
rrr = zinter
(rr,xt,yt,zt,st)
ELSE
rrr = sigma
(zsigma,sf(n),tf(n))
ENDIF
array_qual_traj(n,ind_traj,5) = tf(n)
array_qual_traj(n,ind_traj,6) = sf(n)
array_qual_traj(n,ind_traj,7) = rrr
!! WRITE(lun_traj,7055)n,x0,y0,z0,tnext/tcyc,tf(n),sf(n),rrr
!! ELSE
!! WRITE(lun_traj,7055)n,x0,y0,z0,tnext/tcyc
ENDIF
ENDIF
tnext = tnext + tlap
ENDDO DL2
ENDIF
!====================================================================
! time index update of the velocity field if needed
!====================================================================
IF (iswap /= 0) THEN
IF (TRIM(forback) == 'forward') THEN
l0(1) = l0(1) + 1
ELSE
l0(1) = l0(1) - 1
ENDIF
IF (l0(1) < 1) l0(1) = l0(1) + lmt
IF (l0(1) > lmt) l0(1) = l0(1) - lmt
tswap = tswap + tfic
ENDIF
!====================================================================
! swap for space and time index positions
!====================================================================
gi(1) = hi(n)
gj(1) = hj(n)
gk(1) = hk(n)
gl(1) = hl(n)
age = agenew(n)
!=============!
! game over ? !
!=============!
ENDDO DL1
!=====================================================================!
! DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 DL1 !
!=====================================================================!
IF (TRIM(mode) == 'quantitative') THEN
!==================================================================
! a final criterion uses "1 + nsect" as the section index
! we flag the existence of such a criterion
!==================================================================
IF (nfin(n) <= 0) THEN
nfin(n) = nsect + 1
icrit1 = 1
ENDIF
!==================================================================
! we distinguish meandering and recirculating particles (criter0)
!==================================================================
IF ((TRIM(bin) == 'nobin').AND.(nfin(n) == 1).AND. &
criter0
(hi(n),hj(n),hk(n),hl(n),i0(1),j0(1),k0(1),l0(1))) nfin(n)=0
!==================================================================
! hereafter,
! nfin=0, if the initial section (with criter0) has been reached
! nfin=1, if the initial section (no criter0) or a final section has
! been reached
! nfin=1+nsect, if a final criterion (criter1) is satisfied
!==================================================================
sectrans(nfin(n)) = sectrans(nfin(n)) + ttr(n)
!==================================================================
! there IS NOW storage of transport projections for "criter1" cases
! (but never try to compute 2D streamfunctions with these fields !)
!==================================================================
IF (.NOT.key_eco) THEN
CALL sub_quant_store_transp
(nfin(n))
ENDIF
!==================================================================
! statistics...
!==================================================================
CALL sub_quant_statistics
(nfin(n), ttr(n), agenew(n), tcyc, &
hi(n), hj(n), hk(n), hl(n), tfi(n), tfj(n), tfk(n), tfl(n), &
i0(1), j0(1), k0(1), l0(1), tf(n), sf(n), rf(n), trueage(n))
!! NG: 15_09_2008
!! NG: Here we store the final depth which depend of time for ROMS
!! NG: and symphonie.
IF (key_roms.OR.key_symphonie) THEN
prof=fz
(gi=hi(n), gj=hj(n), gk=hk(n), il=1)
ELSE
prof=fz
(gk=hk(n))
ENDIF
final_depth(n)=prof
WRITE(lun_output,7000)secname(nfin(n)),n, &
i0(1),j0(1),k0(1),l0(1), &
tfi(n),tfj(n),tfk(n),tfl(n)
7000 FORMAT(a14,', #',i0,' en:', 1x,i0,1x,i0,1x,i0,1x,i0,', init.=',4(1x,f0.2),': ',f0.2)
ENDIF
!================
! shoot again...
!================
ENDDO DL0 !=> DO n=1,ntraj
!=====================================================================!
! DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 DL0 !
!=====================================================================!
!======================================================================!
!END END END - General loop on the number of trajectories - END END END!
!======================================================================!
WRITE(lun_standard,*)''
WRITE(lun_standard,*)'-----------------------'
WRITE(lun_standard,*)'= Writing Output Data ='
WRITE(lun_standard,*)'-----------------------'
!!===============================================================================!!
!!============================== QUALITATIVE MODE ===============================!!
!!===============================================================================!!
IF (TRIM(mode) == 'qualitative' ) THEN
!----------------!
!- Trajectories -!
!----------------!
!- Netcdf -!
!----------!
IF (key_alltracers) THEN
CALL sub_save_netcdf_trajectories
( &
array_qual_traj(:,:,1), array_qual_traj(:,:,2), &
array_qual_traj(:,:,3), array_qual_traj(:,:,4), &
array_qual_traj(:,:,5), array_qual_traj(:,:,6), &
array_qual_traj(:,:,7) )
ELSE
CALL sub_save_netcdf_trajectories
( &
array_qual_traj(:,:,1), array_qual_traj(:,:,2), &
array_qual_traj(:,:,3), array_qual_traj(:,:,4) )
ENDIF
!---------!
!- ASCII -!
!---------!
IF (key_ascii_outputs) THEN
DO n = 1, ntraj
DO i = 1, nb_output+1
WRITE(lun_traj,7055)n,array_qual_traj(n,i,1:nstat)
ENDDO
ENDDO
ENDIF
7055 FORMAT (i0,3(1x,f0.5),1x,f0.5,1x,3(1x,f0.5))
!!================================================================================!!
!!============================== QUANTITATIVE MODE ===============================!!
!!================================================================================!!
ELSE
!-------------------!
!- Final positions -!
!-------------------!
!---------!
!- ASCII -!
!---------!
IF (key_ascii_outputs) THEN
DO n = 1, ntraj
IF (key_alltracers) THEN
WRITE(lun_fin_pos,7059)hi(n),hj(n),-hk(n),hl(n), &
ttr(n),trueage(n),nfin(n), tf(n),sf(n),rf(n)
ELSE
WRITE(lun_fin_pos,7059)hi(n),hj(n),-hk(n),hl(n), &
ttr(n),trueage(n),nfin(n)
ENDIF
ENDDO
ENDIF
7059 FORMAT(3(1x,f0.3),1x,f0.3,1x,f0.3,1x,f0.3,1x,i0,3(1x,f0.3))
!----------------------------!
!- statistics file in ASCII -!
!----------------------------!
WRITE(lun_stats,*)' '
IF (key_unitm3) THEN
WRITE(lun_stats,*)'total transport (in m3/s): ', trtot/REAL(lmt,kind=rprec), &
' ( x lmt =' ,trtot, ')'
WRITE(lun_stats,*)'max_transport (in m3/s) : ', max_transport
ELSE
WRITE(lun_stats,*)'total transport (in Sv) : ', trtot / 1.e6_rprec
WRITE(lun_stats,*)'max_transport (in Sv) : ', max_transport / 1.e6_rprec
ENDIF
WRITE(lun_stats,*)'# particles : ' , ntraj
IF ((TRIM(bin) == 'nobin').AND.(isn(-1) > 0)) THEN
!------------------------------------------
! we mask uninitialized min. and max. ages
!------------------------------------------
DO i = 1, nstat
IF (ABS(s1min(i,-1)) > 1.e10_rprec) s1min(i,-1) = 0._rprec
IF (ABS(s1max(i,-1)) > 1.e10_rprec) s1max(i,-1) = 0._rprec
IF (ABS(s0min(i,-1)) > 1.e10_rprec) s0min(i,-1) = 0._rprec
IF (ABS(s0max(i,-1)) > 1.e10_rprec) s0max(i,-1) = 0._rprec
END DO
WRITE(lun_stats,*) ' '
WRITE(lun_stats,7057)'initial state ',' ','#',isn(-1)
WRITE(lun_stats,7157)' stats. for: ',(chpstat(i),i=1,nstat)
WRITE(lun_stats,7257)' min: ' ,(s0min(i,-1),i=1,nstat)
WRITE(lun_stats,7257)' max: ' ,(s0max(i,-1),i=1,nstat)
WRITE(lun_stats,7257)'mean: ' ,(s0x(i,-1)/sn(-1),i=1,nstat)
WRITE(lun_stats,7257)'std. dev.: ' ,&
(SQRT(ABS((s0x2(i,-1)-s0x(i,-1)*s0x(i,-1)/sn(-1))/(sn(-1)-1._rprec))),i=1,nstat)
ENDIF
IF (TRIM(bin) /= 'nobin') THEN
WRITE(lun_stats,*)' '
WRITE(lun_stats,*)'no stats for initial state available since'
WRITE(lun_stats,*)' automatic positioning was by-passed'
ENDIF
WRITE(lun_stats,*)' '
!================================
! 2D flow projection and z stats
!================================
subtime = REAL(lmt, kind=rprec) / REAL(1+lmax-lmin, kind =rprec)
!NG Performances (sub_quant_store_transp => * r_lmt)
r_lmt = 1._rprec / REAL(lmt, kind=rprec)
!!NG IF (key_eco) THEN
!!NG WRITE(lun_xy_zonal)((subtime*uxy(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_merid)((subtime*vxy(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
!!NG WRITE(lun_yz_merid)((subtime*vyz(j,k,1)*r_lmt,j=1,jmt),k=1,kmt)
!!NG WRITE(lun_yz_verti)((subtime*wyz(j,k,1)*r_lmt,j=1,jmt),k=1,kmt)
!!NG WRITE(lun_xz_zonal)((subtime*uxz(i,k,1)*r_lmt,i=1,imt),k=1,kmt)
!!NG WRITE(lun_xz_verti)((subtime*wxz(i,k,1)*r_lmt,i=1,imt),k=1,kmt)
!!NG WRITE(lun_xy_stats)((uh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((vh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((zuh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((zvh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((z2uh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((z2vh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG
!!NG IF (key_alltracers) THEN
!!NG WRITE(lun_xy_stats)((tuh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((tvh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((t2uh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((t2vh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((suh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((svh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((s2uh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((s2vh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((ruh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((rvh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((r2uh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG WRITE(lun_xy_stats)((r2vh(i,j,1)*r_lmt ,i=1,imt),j=1,jmt)
!!NG ENDIF
!!NG
!!NG ELSE
!!NG
!!NG WRITE(lun_xy_zonal)(((subtime*uxy(i,j,n)*r_lmt,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_merid)(((subtime*vxy(i,j,n)*r_lmt,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_yz_merid)(((subtime*vyz(j,k,n)*r_lmt,j=1,jmt),k=1,kmt),n=0,nsect2)
!!NG WRITE(lun_yz_verti)(((subtime*wyz(j,k,n)*r_lmt,j=1,jmt),k=1,kmt),n=0,nsect2)
!!NG WRITE(lun_xz_zonal)(((subtime*uxz(i,k,n)*r_lmt,i=1,imt),k=1,kmt),n=0,nsect2)
!!NG WRITE(lun_xz_verti)(((subtime*wxz(i,k,n)*r_lmt,i=1,imt),k=1,kmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((uh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((vh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((zuh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((zvh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((z2uh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((z2vh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG
!!NG IF (key_alltracers) THEN
!!NG WRITE(lun_xy_stats)(((tuh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((tvh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((t2uh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((t2vh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((suh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((svh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((s2uh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((s2vh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((ruh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((rvh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((r2uh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG WRITE(lun_xy_stats)(((r2vh(i,j,n)*r_lmt ,i=1,imt),j=1,jmt),n=0,nsect2)
!!NG ENDIF
!!NG
!!NG ENDIF
uxy(:,:,:) = uxy(:,:,:) * subtime * r_lmt
vxy(:,:,:) = vxy(:,:,:) * subtime * r_lmt
vyz(:,:,:) = vyz(:,:,:) * subtime * r_lmt
wyz(:,:,:) = wyz(:,:,:) * subtime * r_lmt
uxz(:,:,:) = uxz(:,:,:) * subtime * r_lmt
wxz(:,:,:) = wxz(:,:,:) * subtime * r_lmt
uh(:,:,:) = uh(:,:,:) * r_lmt
vh(:,:,:) = vh(:,:,:) * r_lmt
zuh(:,:,:) = zuh(:,:,:) * r_lmt
zvh(:,:,:) = zvh(:,:,:) * r_lmt
z2uh(:,:,:) = z2uh(:,:,:) * r_lmt
z2vh(:,:,:) = z2vh(:,:,:) * r_lmt
IF (key_alltracers) THEN
tuh(:,:,:) = tuh(:,:,:) * r_lmt
tvh(:,:,:) = tvh(:,:,:) * r_lmt
t2uh(:,:,:) = t2uh(:,:,:) * r_lmt
t2vh(:,:,:) = t2vh(:,:,:) * r_lmt
suh(:,:,:) = suh(:,:,:) * r_lmt
svh(:,:,:) = svh(:,:,:) * r_lmt
s2uh(:,:,:) = s2uh(:,:,:) * r_lmt
s2vh(:,:,:) = s2vh(:,:,:) * r_lmt
ruh(:,:,:) = ruh(:,:,:) * r_lmt
rvh(:,:,:) = rvh(:,:,:) * r_lmt
r2uh(:,:,:) = r2uh(:,:,:) * r_lmt
r2vh(:,:,:) = r2vh(:,:,:) * r_lmt
ENDIF
!!---------------------------!
!!= Quantitative statistics =!
!!---------------------------!
trtot1 = 0._rprec
nsect2 = nsect
IF (icrit1 == 1) THEN
nsect2 = 1 + nsect
ENDIF
DO n = 0, nsect2
IF (key_unitm3) THEN
WRITE(lun_stats,7357)secname(n), subtime*sectrans(n)*r_lmt, n
ELSE
WRITE(lun_stats,7357)secname(n), subtime*sectrans(n)/1.e6_rprec*r_lmt, n
ENDIF
trtot1 = trtot1 + sectrans(n)
END DO
IF (key_unitm3) THEN
WRITE(lun_stats,7357)'total' ,subtime*trtot*r_lmt
WRITE(lun_stats,7357)'except mnds',subtime*(trtot-sectrans(0))*r_lmt
WRITE(lun_stats,7357)'lost' ,subtime*(trtot-trtot1)*r_lmt
ELSE
WRITE(lun_stats,7357)'total' ,subtime*trtot/1.e6_rprec*r_lmt
WRITE(lun_stats,7357)'except mnds',subtime*(trtot-sectrans(0))/1.e6_rprec*r_lmt
WRITE(lun_stats,7357)'lost' ,subtime*(trtot-trtot1)/1.e6_rprec*r_lmt
ENDIF
7357 FORMAT(a15,1x,f0.4,1x,i0)
DO n=0,nsect2
IF (isn(n) > 0) THEN
WRITE(lun_stats,*)' '
WRITE(lun_stats,7057)'final state ',TRIM(secname(n)),' # ',isn(n)
WRITE(lun_stats,7157)'stats. ini: ',(chpstat(i),i=1,nstat)
WRITE(lun_stats,7257)' min: ',(s0min(i,n),i=1,nstat)
WRITE(lun_stats,7257)' max: ',(s0max(i,n),i=1,nstat)
WRITE(lun_stats,7257)'mean: ',(s0x(i,n)/sn(n),i=1,nstat)
WRITE(lun_stats,7257)'std. dev.: ',( SQRT( &
ABS((s0x2(i,n)-s0x(i,n)*s0x(i,n)/sn(n))/(sn(n)-1._rprec))),i=1,nstat)
WRITE(lun_stats,7157)'stats. fin: ',(chpstat(i),i=1,nstat)
WRITE(lun_stats,7257)' min: ',(s1min(i,n),i=1,nstat)
WRITE(lun_stats,7257)' max: ',(s1max(i,n),i=1,nstat)
WRITE(lun_stats,7257)'mean: ',(s1x(i,n)/sn(n),i=1,nstat)
WRITE(lun_stats,7257)'std. dev.: ',( SQRT( &
ABS((s1x2(i,n)-s1x(i,n)*s1x(i,n)/sn(n))/(sn(n)-1._rprec))),i=1,nstat)
ENDIF
END DO
!----------------------------------------------------------------------!
!- Define statistic variables and their attributes in quantative mode -!
!----------------------------------------------------------------------!
CALL sub_save_netcdf_data_init_stats
()
!------------------------------------!
!- Save statistics in Netcdf Format -!
!------------------------------------!
Call sub_save_netcdf_stats
()
!----------------------------!
!- Close Netcdf output file -!
!----------------------------!
WRITE(lun_standard,*)''
WRITE(lun_standard,*)'--------------------------------'
WRITE(lun_standard,*)'= Close Statistics NetCDF File ='
WRITE(lun_standard,*)'--------------------------------'
CALL sub_save_netcdf_data_close
(ncid_data_stats)
ENDIF
7057 FORMAT(a13,a15,a2,i7)
7157 FORMAT(a13,7a10)
7257 FORMAT(a13,7f10.3)
!=====================================!
!- Binary storage of final positions -!
!=====================================!
IF (key_ascii_outputs) THEN
DO n = 1, ntraj
WRITE(lun_final)hi(n),hj(n),-hk(n),hl(n),ttr(n),agenew(n)
WRITE(lun_init)tfi(n),tfj(n),tfk(n),tfl(n),ttr(n),tage(n)
ENDDO
ENDIF
!! NG: 15_09_2008
IF (TRIM(mode) == 'quantitative') THEN
!-----------------!
!- Initial Depth -!
!-----------------!
!- Netcdf -!
!----------!
CALL sub_save_netcdf_init_depth
(init_depth(:))
ENDIF
!-------------------!
!- Final positions -!
!-------------------!
!- Netcdf -!
!----------!
IF (key_alltracers) THEN
CALL sub_save_netcdf_final_pos
( &
hi(:),hj(:),-hk(:),hl(:) , &
agenew(:), nfin(:) , &
final_depth(:) , &
tf(:), sf(:), rf(:) )
ELSE
CALL sub_save_netcdf_final_pos
( &
hi(:),hj(:),-hk(:),hl(:) , &
agenew(:), nfin(:),final_depth(1:ntraj) )
ENDIF
!-----------------------------!
!- Close Netcdf output files -!
!-----------------------------!
WRITE(lun_standard,*)''
WRITE(lun_standard,*)'-------------------------------'
WRITE(lun_standard,*)'= Close Positions NetCDF File ='
WRITE(lun_standard,*)'-------------------------------'
CALL sub_save_netcdf_data_close
(ncid_data_pos)
!=====================!
!- DEALLOCATE memory -!
!=====================!
WRITE(lun_standard,*)''
WRITE(lun_standard,*)'---------------------'
WRITE(lun_standard,*)'= Deallocate Memory ='
WRITE(lun_standard,*)'---------------------'
IF (ALLOCATED(i0)) THEN
CALL sub_memory
(-size(i0),'i','i0','trajec_seq')
DEALLOCATE(i0)
END IF
IF (ALLOCATED(j0)) THEN
CALL sub_memory
(-size(j0),'i','j0','trajec_seq')
DEALLOCATE(j0)
END IF
IF (ALLOCATED(k0)) THEN
CALL sub_memory
(-size(k0),'i','k0','trajec_seq')
DEALLOCATE(k0)
END IF
IF (ALLOCATED(l0)) THEN
CALL sub_memory
(-size(l0),'i','l0','trajec_seq')
DEALLOCATE(l0)
END IF
IF (ALLOCATED(nfin)) THEN
CALL sub_memory
(-size(nfin),'i','nfin','trajec_seq')
DEALLOCATE(nfin)
END IF
IF (ALLOCATED(agenew)) THEN
CALL sub_memory
(-size(agenew),'r','agenew','trajec_seq')
DEALLOCATE(agenew)
END IF
IF (ALLOCATED(trueage)) THEN
CALL sub_memory
(-size(trueage),'r','truage','trajec_seq')
DEALLOCATE(trueage)
END IF
IF (ALLOCATED(gi)) THEN
CALL sub_memory
(-size(gi),'r','gi','trajec_seq')
DEALLOCATE(gi)
END IF
IF (ALLOCATED(gj)) THEN
CALL sub_memory
(-size(gj),'r','gj','trajec_seq')
DEALLOCATE(gj)
END IF
IF (ALLOCATED(gk)) THEN
CALL sub_memory
(-size(gk),'r','gk','trajec_seq')
DEALLOCATE(gk)
END IF
IF (ALLOCATED(gl)) THEN
CALL sub_memory
(-size(gl),'r','gl','trajec_seq')
DEALLOCATE(gl)
END IF
IF (ALLOCATED(hi)) THEN
CALL sub_memory
(-size(hi),'r','hi','trajec_seq')
DEALLOCATE(hi)
END IF
IF (ALLOCATED(hj)) THEN
CALL sub_memory
(-size(hj),'r','hj','trajec_seq')
DEALLOCATE(hj)
END IF
IF (ALLOCATED(hk)) THEN
CALL sub_memory
(-size(hk),'r','hk','trajec_seq')
DEALLOCATE(hk)
END IF
IF (ALLOCATED(hl)) THEN
CALL sub_memory
(-size(hl),'r','hl','trajec_seq')
DEALLOCATE(hl)
END IF
IF (ALLOCATED(tf)) THEN
CALL sub_memory
(-size(tf),'r','tf','trajec_seq')
DEALLOCATE(tf)
END IF
IF (ALLOCATED(sf)) THEN
CALL sub_memory
(-size(sf),'r','sf','trajec_seq')
DEALLOCATE(sf)
END IF
IF (ALLOCATED(rf)) THEN
CALL sub_memory
(-size(rf),'r','rf','trajec_seq')
DEALLOCATE(rf)
END IF
!! NG: 15_09_2008
IF (ALLOCATED(init_depth)) THEN
CALL sub_memory
(-size(init_depth),'i','init_depth','trajec_seq')
DEALLOCATE(init_depth)
END IF
IF (ALLOCATED(final_depth)) THEN
CALL sub_memory
(-size(final_depth),'i','final_depth','trajec_seq')
DEALLOCATE(final_depth)
END IF
IF (ALLOCATED(array_qual_traj)) THEN
CALL sub_memory
(-size(array_qual_traj),'r','array_qual_traj','trajec_seq')
DEALLOCATE(array_qual_traj)
END IF
CALL sub_netcdf_dealloc_mem
()
CALL sub_input_data_dealloc_mem
()
CALL sub_coord_dealloc
()
CALL sub_scalef_dealloc
()
IF (key_roms) CALL sub_h_roms_dealloc
()
CALL sub_transp_dealloc
()
IF (key_alltracers) THEN
CALL sub_tracer_dealloc
()
ENDIF
CALL sub_posin_dealloc
()
CALL sub_tmask_dealloc
()
IF (TRIM(mode) == 'quantitative') THEN ! [QUANT]
CALL sub_txt_dealloc
()
CALL sub_flags_dealloc
()
CALL sub_stats_dealloc
()
CALL sub_stati_dealloc
()
CALL sub_quant_dealloc
()
ENDIF
!=================
! MAIN PART - END
!=================
END SUBROUTINE trajec
!!***