-
Matthias Schmiddunser authored
Self-Scheduling Model for multiple processors.
Matthias Schmiddunser authoredSelf-Scheduling Model for multiple processors.
post_tracking.f90 43.40 KiB
module definitions
type tracking_cloud
integer np,ntracks
real, dimension(:), pointer :: time
real, dimension(:,:),pointer :: x,y,z,press,temp,velo,veloZ
integer,dimension(:),pointer :: step
integer,dimension(:,:),pointer :: tracks
end type tracking_cloud
end module definitions
!==============================================================================!
!==============================================================================!
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!==============================================================================!
!==============================================================================!
! |
! TRACKING CLOUD POST PROCESSING Jan. 2015 |
! written by Matthias Schmiddunser |
! using age calculation form Pecube and Terra |
! |
!==============================================================================!
!==============================================================================!
! |
! This routine reads in the tracking output files, extracts the surfaced |
! tracks, does variable calculation and writes ASCII vtk files. |
! |
! CAUTION: Particle ID (Douar) and Point Number (VTK) differ: DOUAR numbering |
! starts at 1, VTK numbering starts at 0. |
program post_tracking
use definitions
!use m_ageCalculation
implicit none
include 'mpif.h'
type (tracking_cloud) :: trcl
!input parameters from input_tracking_output.txt
character(150) :: input_dir
character(150) :: output_dir
integer :: output_mode
integer :: output_scaled
integer :: output_first
integer :: output_last
integer :: output_step
integer :: keep_nstep
integer :: min_track_length
integer :: input_first
integer :: input_step
real :: pad_age
integer :: output_AHe
integer :: output_ZHe
integer :: output_AFT
integer :: output_ZFT
integer :: output_Muscovite
integer :: output_velo
integer :: output_veloz
integer :: output_tempinteger :: output_press
integer :: output_time
integer :: output_istep
integer :: output_ID
!buffer variables for reading
real :: read_time,lref,veloref,rhoref,tempscale
integer :: read_step,read_ntracks
integer,dimension(:,:),allocatable :: read_tracks
integer,dimension(:),allocatable :: keep_tracks
!mpi and data transfer variables
logical :: flag
integer :: provided,iproc,nproc,ierr,maxtag
integer,dimension(MPI_STATUS_SIZE) :: mpistatus
integer :: sender,tracksent,send_nstep,endstep
real,dimension(:),allocatable :: send_time,send_temp
real,dimension(:),allocatable :: ages_AHe,ages_ZHe,ages_AFT,ages_ZFT,ages_Muscovite
real,dimension(:),allocatable :: allages
!age calculation
real(8) :: apatiteHe, zirconHe
real(8) :: zirconFT, apatiteFT
real(8) :: muscovite_age
!other variables
integer :: i,k,l,err,dummyint,ipoint,itrack
integer :: output_nstep,input_nstep,input_last,out_track_any,out_nages
integer :: write_nstep,write_points,write_tracks
real :: dummyreal
logical :: keepme
character(4) :: cistep
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
provided=0
call mpi_init_thread(MPI_THREAD_MULTIPLE,provided,ierr)
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
if (iproc==0) then
write (*,'(a)') '!------------------------------------------------------------------------------|'
write (*,'(a)') '!------------------------------------------------------------------------------|'
write (*,'(a)') '! |'
write (*,'(a)') '! 8888888b. .d88888b. 888 888 d8888 8888888b. |'
write (*,'(a)') '! 888 "Y88b d88P" "Y88b 888 888 d88888 888 Y88b |'
write (*,'(a)') '! 888 888 888 888 888 888 d88P888 888 888 |'
write (*,'(a)') '! 888 888 888 888 888 888 d88P 888 888 d88P |'
write (*,'(a)') '! 888 888 888 888 888 888 d88P 888 8888888P" |'
write (*,'(a)') '! 888 888 888 888 888 888 d88P 888 888 T88b |'
write (*,'(a)') '! 888 .d88P Y88b. .d88P Y88b. .d88P d8888888888 888 T88b |'
write (*,'(a)') '! 8888888P" "Y88888P" "Y88888P" d88P 888 888 T88b |'
write (*,'(a)') '! |'
write (*,'(a)') '! DOUAR-WSMP post-processor for the tracking cloud |'
write (*,'(a)') '!------------------------------------------------------------------------------|'
write (*,'(a)') '!------------------------------------------------------------------------------|'
write (*,*) ''
end if
!------------------------------------------------------------------------------|
! Input parameters
!------------------------------------------------------------------------------|
!Reading input parameters from post-input-file:
if (iproc==0) then
open(unit=10, file='input_tracking_outputs.txt',status='old')
read (10,*) input_dir
read (10,*) output_dir
read (10,*) output_mode
read (10,*) output_scaled
read (10,*) output_first
read (10,*) output_last
read (10,*) output_step
read (10,*) keep_nstep
read (10,*) min_track_length
!keep_nstep has to be less than output_step
keep_nstep = min(keep_nstep,output_step-1)
read (10,*) input_first
read (10,*) input_step
read (10,*) !here come the age parameters
read (10,*) pad_age
read (10,*) output_AHe
read (10,*) output_ZHe
read (10,*) output_AFT
read (10,*) output_ZFT
read (10,*) output_Muscovite
read (10,*) !here come the track parameters
read (10,*) output_velo
read (10,*) output_veloz
read (10,*) output_temp
read (10,*) output_press
read (10,*) output_time
read (10,*) output_istep
read (10,*) output_ID
close (10)
!calculating and andjusting input and output timesteps
output_nstep = (output_last-output_first)/output_step+1
output_last = output_first+(output_nstep-1)*output_step
input_last = output_last
input_nstep = (input_last-input_first)/input_step+1
input_last = input_first+(input_nstep-1)*input_step
!Checking for slash at end of input and output dir
l = len(trim(input_dir))
if (input_dir(l:l) .ne. '/') then
input_dir = input_dir(1:l)//'/'
end if
l = len(trim(output_dir))
if (output_dir(l:l) .ne. '/') then
output_dir = output_dir(1:l)//'/'
end if
end if
!Distribute output parameters to other processes
call mpi_bcast(output_mode,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_AHe,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_ZHe,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_AFT,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_ZFT,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_Muscovite,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_velo,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_veloz,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_temp,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_press,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_time,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_istep,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(output_ID,1,mpi_integer,0,mpi_comm_world,ierr)
!determine number of ages that will be calculated
out_nages = output_AHe+output_ZHe+output_AFT+output_ZFT+output_Muscovite
allocate (allages(out_nages),stat=err); if (err.ne.0) call stop_run ('Error alloc allages')
write (*,*) iproc, out_nages
!writing output parameters to screen
if (iproc==0) then
write (*,'(a,a)') 'Simulation Directory: ',input_dir
write (*,'(a,i2,a,i3,a,i3,a,i3)') 'Input Timesteps: Reading ',input_nstep,' input files '// &
'from step ',input_first,' to ',input_last,' in steps of ',input_step
write (*,'(a,a)') 'Output Directory: ',output_dir
write (*,'(a,i2,a,i3,a,i3,a,i3)') 'Output Timesteps: Generating ',output_nstep,' output files '// &
'from step ',output_first,' to ',output_last,' in steps of ',output_step
write (*,'(a,i2,a)') 'accumulating ',keep_nstep,' steps for surfaced tracks'
write (*,'(a,i1)') 'Output mode: ',output_mode
write (*,'(a,l)') 'Write natural units: ',(output_scaled==1)
write (*,'(a)') ''
write (*,'(a)') ' Age variables: '
write (*,'(a,l)') 'Write AHe Ages: ',(output_AHe==1)
write (*,'(a,l)') 'Write ZHe Ages: ',(output_ZHe==1)
write (*,'(a,l)') 'Write AFT Ages: ',(output_AFT==1)
write (*,'(a,l)') 'Write ZFT Ages: ',(output_ZFT==1)
write (*,'(a,l)') 'Write Muscovite Ages: ',(output_Muscovite==1)
write (*,'(a)') ''
write (*,'(a)') 'Track variables: '
write (*,'(a,l)') 'Write track velocity: ',(output_velo==1)
write (*,'(a,l)') 'Write track z velocity: ',(output_veloz==1)
write (*,'(a,l)') 'Write track temperature: ',(output_temp==1)
write (*,'(a,l)') 'Write track pressure: ',(output_press==1)
write (*,'(a,l)') 'Write track time: ',(output_time==1)
write (*,'(a,l)') 'Write track step number: ',(output_istep==1)
write (*,'(a,l)') 'Write particle/track ID: ',(output_ID==1)
write (*,'(a)') ''
write (*,'(a)') '!------------------------------------------------------------------------------|'
write (*,'(a)') ''
end if
!------------------------------------------------------------------------------|
! Reading Tracking Cloud Files
!------------------------------------------------------------------------------|
if (iproc==0) then
!reading input file header for array sizes
call int_to_char (cistep,4,input_last)
open (unit=20,file=trim(input_dir)//'TROUT/tracking_t'//cistep//'.bin',status='old',form='unformatted')
read (20) trcl%np,read_ntracks,read_step,read_time
read (20) lref,veloref,rhoref,tempscale
close (20)
!allocating master arrays for input
allocate (trcl%x(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%x')
allocate (trcl%y(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%y')
allocate (trcl%z(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%z')
allocate (trcl%temp(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%temp')
allocate (trcl%press(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%press')
allocate (trcl%step(input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%step')
allocate (trcl%time(input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%time')
allocate (read_tracks(read_ntracks,3),stat=err); if (err.ne.0) call stop_run ('Error alloc read_tracks')
!read all input files and save in collective arrays
k = 1
do i = input_first,input_last,input_step
!read all information for this time step
call int_to_char (cistep,4,i)
open (unit=20,file=trim(input_dir)//'TROUT/tracking_t'//cistep//'.bin',status='old',form='unformatted')
read (20) dummyint,dummyint,trcl%step(k),trcl%time(k)
read (20)
read (20) trcl%x(:,k)
read (20) trcl%y(:,k)
read (20) trcl%z(:,k)
read (20) trcl%temp(:,k) read (20) trcl%press(:,k)
!read track information only from last timestep
if (i==input_last) then
write (*,'(a)') 'All surfaced tracks:'
write (*,'(a)') 'track ID, Particle ID, start, end'
do l=1,read_ntracks
read (20) read_tracks(l,:)
write (*,'(4i6)') l, read_tracks(l,:)
end do
write (*,'(a)') ''
end if
!advance column in master array
k = k+1
close (20)
end do
write (*,'(a)') 'Finished reading input files'
write (*,'(a)') ''
!select only tracks that surface at an output timestep or up to keep_nstep
!timesteps before each output timestep and are at least min_track_lenght
!timesteps long.
!Then adjust track timesteps to actual position in master array.
!The offset k has to be recorded in order to associate the track with the
!correct timestep.
allocate (keep_tracks(read_ntracks),stat=err); if (err.ne.0) call stop_run ('Error alloc keep_tracks')
keep_tracks = 0
trcl%ntracks = 0
do k=0,keep_nstep
do i=1,read_ntracks
keepme = read_tracks(i,3) >= output_first-k
keepme = keepme .AND. (mod((read_tracks(i,3)+k-output_first),output_step)==0)
keepme = keepme .AND. read_tracks(i,3)-read_tracks(i,2)+1 >= min_track_length
if (keepme) then
keep_tracks(i) = 1000+k
trcl%ntracks = trcl%ntracks+1
end if
end do
end do
!allocate tracks array of needed length
if (trcl%ntracks > 0) then
allocate (trcl%tracks(trcl%ntracks,7),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%tracks')
trcl%tracks = 0
!copy relevant tracks into new track array
!adjust track timsteps to position in new track array
!always round up for start and round down for end in order not to include parts of other tracks
itrack=1
write (*,'(a)') 'Selected surfaced tracks:'
write (*,'(a)') '*: in master array'
write (*,'(a)') 'Particle ID, start*, end*, length*, output timestep, output step*, track ID'
do i=1,read_ntracks
if (keep_tracks(i)>0) then
!1st column is particle ID
trcl%tracks(itrack,1) = read_tracks(i,1)
!2nd column is start in master array, round up
trcl%tracks(itrack,2)= max(1,ceiling(1.0+(read_tracks(i,2)-input_first)/(1.0*input_step)))
!3rd columd is end in master array, round down
trcl%tracks(itrack,3)= max(1,1+(read_tracks(i,3)-input_first)/input_step) !integer division
!4th column is track length (number of points in master array)
trcl%tracks(itrack,4) = trcl%tracks(itrack,3)-trcl%tracks(itrack,2)+1
!5th column is Douar timestep of surfacing, upshifted for those tracks
!within keep_nstep
trcl%tracks(itrack,5) = read_tracks(i,3)+(keep_tracks(i)-1000) !6th column is output timestep in master array
trcl%tracks(itrack,6)= max(1,1+(trcl%tracks(itrack,5)-input_first)/input_step) !integer division
!7th column is track ID, i.e. the row of the original track array from input file
trcl%tracks(itrack,7) = i
write (*,'(7i6)') trcl%tracks(itrack,:)
itrack = itrack+1
end if
end do
else
write (*,'(a)') 'No surfaced tracks'
end if
deallocate (read_tracks)
deallocate (keep_tracks)
end if
!call MPI_barrier (mpi_comm_world, ierr)
call mpi_bcast(trcl%ntracks,1,mpi_integer,0,mpi_comm_world,ierr)
!------------------------------------------------------------------------------|
! Calculate Thermochron Ages
!------------------------------------------------------------------------------|
ages: if ((output_mode==1 .or. output_mode==2) .AND. trcl%ntracks>0) then
call MPI_COMM_GET_ATTR (MPI_COMM_WORLD, MPI_TAG_UB, maxtag, flag, ierr)
if (trcl%ntracks > maxtag) call stop_run ('Too many tracks for safe communications')
workshare: if (iproc==0) then
!master distributing work
if (output_AHe==1) allocate (ages_AHe(trcl%ntracks))
if (output_ZHe==1) allocate (ages_ZHe(trcl%ntracks))
if (output_AFT==1) allocate (ages_AFT(trcl%ntracks))
if (output_ZFT==1) allocate (ages_ZFT(trcl%ntracks))
if (output_Muscovite==1) allocate (ages_Muscovite(trcl%ntracks))
tracksent = 0
!wait for answer, distribute new track
workdistr: do i=2-nproc,trcl%ntracks
!for first nproc-1 times, send to processes without receive
!distribute one track to each processor
if (i>0) then
!get result from other process
call MPI_RECV (allages,out_nages,MPI_REAL, MPI_ANY_SOURCE,MPI_ANY_TAG, &
MPI_COMM_WORLD,mpistatus,ierr)
!Identify process and track ID
sender = mpistatus(MPI_SOURCE)
itrack = mpistatus(MPI_TAG)
!write (*,*) iproc, ': received track',itrack,' from proc ',sender, 'length = ',out_nages
!Write into result arrays
k = 1
if (output_AHe==1) then
ages_AHe(itrack) = allages(k)
call correct_number(ages_AHe(itrack))
k = k+1
end if
if (output_ZHe==1) then
ages_ZHe(itrack) = allages(k)
call correct_number(ages_ZHe(itrack))
k = k+1
end if
if (output_AFT==1) then
ages_AFT(itrack) = allages(k)
call correct_number(ages_AFT(itrack))
k = k+1
end if
if (output_ZFT==1) then ages_ZFT(itrack) = allages(k)
call correct_number(ages_ZFT(itrack))
k = k+1
end if
if (output_Muscovite==1) then
ages_Muscovite(itrack) = allages(k)
call correct_number(ages_Muscovite(itrack))
end if
else
!effecively loops from 1 to nproc-1 for sender
sender = i+nproc-1
end if
!Send new information
if (tracksent < trcl%ntracks) then
!Length of time and temperature arrays
!Add two entrys at front for initalisation
!If surfaced earlier than output timestep, add another entry at end
endstep = min(1,trcl%tracks(tracksent+1,6)-trcl%tracks(tracksent+1,3))
send_nstep = 2 + trcl%tracks(tracksent+1,4)+endstep
allocate (send_time(send_nstep))
allocate (send_temp(send_nstep))
!pad send-arrays at beginning, part 1
send_time(3:send_nstep-endstep) = &
trcl%time(trcl%tracks(tracksent+1,2):trcl%tracks(tracksent+1,3))
send_time(1:2) = (/0.0,0.0/)
send_temp(3:send_nstep-endstep) = trcl%temp(trcl%tracks(tracksent+1,1), &
trcl%tracks(tracksent+1,2):trcl%tracks(tracksent+1,3))
send_temp(1) = send_temp(3)
send_temp(2) = send_temp(3)
!pad send arrays at end if necessary
if (endstep>0) then
send_time(send_nstep) = trcl%time(trcl%tracks(tracksent+1,6))
send_temp(send_nstep) = send_temp(send_nstep-1)
end if
!Calculate relative Age in Ma
send_time = send_time(send_nstep) - send_time
send_time = send_time*lref/veloref
!scale temperature
send_temp = send_temp*tempscale
!pad age array at the beginning, part 2
send_time(1) = pad_age
send_time(2) = 0.5*(send_time(1)+send_time(3))
!write (*,*) iproc, ': send track',tracksent+1,' of length',send_nstep,' to proc ',sender
call MPI_SEND (send_time,send_nstep,MPI_REAL,sender,tracksent+1,MPI_COMM_WORLD,ierr)
call MPI_SEND (send_temp,send_nstep,MPI_REAL,sender,tracksent+1,MPI_COMM_WORLD,ierr)
tracksent = tracksent+1
deallocate (send_time)
deallocate (send_temp)
else
!send end signal
call MPI_SEND (0.0,1,MPI_REAL,sender,0,MPI_COMM_WORLD,ierr)
end if
end do workdistr
else workshare
!slaves doing actual calculation
agecalc: do
call MPI_probe (0, MPI_ANY_TAG, MPI_COMM_WORLD, mpistatus, ierr)
itrack = mpistatus(MPI_TAG) !write (*,*) iproc, ': probed track', itrack
!tag 0 is end signal
if (itrack==0) then
!write (*,*) iproc, ': receive end signal'
!empty send buffer - necessary?
call MPI_RECV (dummyreal,1,MPI_REAL,0,MPI_ANY_TAG,MPI_COMM_WORLD,mpistatus,ierr)
exit
end if
call MPI_get_count (mpistatus,MPI_REAL,send_nstep,ierr)
allocate (send_time(send_nstep))
allocate (send_temp(send_nstep))
!write (*,*) iproc, ': receiving track', itrack, ' of length', send_nstep
call MPI_RECV (send_time,send_nstep,MPI_REAL,0,MPI_ANY_TAG,MPI_COMM_WORLD,mpistatus,ierr)
call MPI_RECV (send_temp,send_nstep,MPI_REAL,0,MPI_ANY_TAG,MPI_COMM_WORLD,mpistatus,ierr)
!do age calculation
k=1
if (output_AHe==1 .or. output_ZHe==1) then
!call Mad_He(dble(send_time), dble(send_temp), send_nstep, apatiteHe, zirconHe)
apatiteHe = 1.d0
zirconHe = 2.d0
if (output_AHe==1) then
allages(k) = real(apatiteHe)
k = k+1
end if
if (output_ZHe==1) then
allages(k) = real(zirconHe)
k = k+1
end if
end if
if (output_ZFT==1) then
!call ZFT(dble(send_time), dble(send_temp), send_nstep, zirconFT)
zirconFT = 3.d0
allages(k) = real(zirconFT)
k = k+1
end if
if (output_AFT==1) then
!call get_aft_age(dble(send_time), dble(send_temp), send_nstep, apatiteFT)
apatiteFT = 4.d0
allages(k) = real(apatiteFT)
k = k+1
end if
if (output_Muscovite==1) then
!call Muscovite(dble(send_time), dble(send_temp), send_nstep, muscovite_age)
muscovite_age = 5.d0
allages(k) = real(muscovite_age)
k = k+1
end if
!send results
!write (*,*) iproc, ': send results of track', itrack, 'length = ',out_nages
call MPI_SEND (allages,out_nages,MPI_REAL,0,itrack, MPI_COMM_WORLD,ierr)
deallocate (send_time)
deallocate (send_temp)
end do agecalc
end if workshare
end if ages
!------------------------------------------------------------------------------|
! Scaling Values to natural units!------------------------------------------------------------------------------|
!check what age prediction needs...
!units of scaling parameters:
!lref: km
!veloref: mm/yr = km/Myr
!rhoref: kg/m^3
!tempscale: °C
if (iproc==0) then
if (output_scaled==1) then
do k=1,input_nstep
do i=1,trcl%np
trcl%x(i,k) = trcl%x(i,k)*lref !pos in km
trcl%y(i,k) = trcl%y(i,k)*lref
trcl%z(i,k) = trcl%z(i,k)*lref
trcl%temp(i,k) = trcl%temp(i,k)*tempscale !temp in C
trcl%press(i,k) = trcl%press(i,k)*1.0e-3*9.81*rhoref*lref !pressure in MPa
end do
trcl%time(k) = trcl%time(k)*lref/veloref !time in Myr
end do
end if
end if
!------------------------------------------------------------------------------|
! Calculating Output values
!------------------------------------------------------------------------------|
if (iproc==0) then
!track values
velo: if (output_velo==1) then
allocate (trcl%velo(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%velo')
do i=1,trcl%np
trcl%velo(i,1) = sqrt( &
(trcl%x(i,2)-trcl%x(i,1))**2+(trcl%y(i,2)-trcl%y(i,1))**2+(trcl%z(i,2)-trcl%z(i,1))**2) &
/(trcl%time(2)-trcl%time(1))
end do
do i=1,trcl%np
do k=2,input_nstep
trcl%velo(i,k) = sqrt( &
(trcl%x(i,k)-trcl%x(i,k-1))**2+(trcl%y(i,k)-trcl%y(i,k-1))**2+(trcl%z(i,k)-trcl%z(i,k-1))**2) &
/(trcl%time(k)-trcl%time(k-1))
end do
end do
end if velo
veloz: if (output_veloz==1) then
allocate (trcl%veloz(trcl%np,input_nstep),stat=err); if (err.ne.0) call stop_run ('Error alloc trcl%veloz')
do i=1,trcl%np
trcl%veloz(i,1) = (trcl%z(i,2)-trcl%z(i,1))/(trcl%time(2)-trcl%time(1))
end do
do i=1,trcl%np
do k=2,input_nstep
trcl%veloz(i,k) = (trcl%z(i,k)-trcl%z(i,k-1))/(trcl%time(k)-trcl%time(k-1))
end do
end do
end if veloz
end if
!------------------------------------------------------------------------------|
! Checking for NaNs and variable ranges
!------------------------------------------------------------------------------|
if (iproc==0) then
do k=1,input_nstep
call correct_number(trcl%time(k))
do i=1,trcl%np
call correct_number(trcl%temp(i,k)) call correct_number(trcl%press(i,k))
if (output_velo==1) call correct_number(trcl%velo(i,k))
if (output_veloz==1) call correct_number(trcl%veloz(i,k))
end do
end do
end if
!------------------------------------------------------------------------------|
! Writing Output Files
!------------------------------------------------------------------------------|
if (iproc==0) then
if ((output_mode==1.or.output_mode==2) .and. trcl%ntracks>0) then
writeages: do l=output_first,output_last,output_step
!count number of tracks for this timestep
write_tracks = 0
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write_tracks = write_tracks+1
end if
end do
if (write_tracks>0) then
!create output file age calculation results
call int_to_char (cistep,4,l)
open (unit=40,file=trim(output_dir)//'ages_t'//cistep//'.vtk',status='replace')
!write header
write(40,'(a)') '# vtk DataFile Version 3.1'
write(40, '(a,i5,a,i4,a)') 'DOUAR thermochron ages for timestep',l, &
', accumulating',keep_nstep,' timesteps'
write(40, '(a)') 'ASCII'
write(40, '(a)') 'DATASET UNSTRUCTURED_GRID'
write(40, '(a)') ''
!write point coordinates
write(40, '(a,i10,a)') 'POINTS ', write_tracks, ' FLOAT'
do i=1,trcl%ntracks
!go through all tracks
if (trcl%tracks(i,5)==l) then
!if track surfaces at this timestep: write last point of track (at surface)
write(40,'(3ES12.5E2)') trcl%x(trcl%tracks(i,1),trcl%tracks(i,3)), &
trcl%y(trcl%tracks(i,1),trcl%tracks(i,3)), &
trcl%z(trcl%tracks(i,1),trcl%tracks(i,3))
end if
end do
write(40, '(a)') ''
!write cell data
write(40, '(a,i10,i10)') 'CELLS ', write_tracks, 2*write_tracks
ipoint = 0
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
!write length of track
write(40, '(i1,i6)') 1, ipoint
!advance point ID
ipoint = ipoint+1
end if
end do
write (40, '(a)') ''
!write cell data type, 1=single vertex
write(40, '(a,i10)') 'CELL_TYPES ',write_tracks
do i = 1,write_tracks
write(40, '(i2)') 1
end do
write (40, '(a)') ''
!write point datasets
if (out_nages > 0) write (40,'(a,i10)') 'POINT_DATA ',write_tracks
if (output_AHe==1) then
write (40,'(a)') 'SCALARS ApatiteHe FLOAT 1'
write (40,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write(40,'(3ES12.5E2)') ages_AHe(i)
end if
end do
write(40, '(a)') ''
end if
if (output_ZHe==1) then
write (40,'(a)') 'SCALARS ZirconHe FLOAT 1'
write (40,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write(40,'(3ES12.5E2)') ages_ZHe(i)
end if
end do
write(40, '(a)') ''
end if
if (output_AFT==1) then
write (40,'(a)') 'SCALARS ApatiteFT FLOAT 1'
write (40,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write(40,'(3ES12.5E2)') ages_AFT(i)
end if
end do
write(40, '(a)') ''
end if
if (output_ZFT==1) then
write (40,'(a)') 'SCALARS ZirconFT FLOAT 1'
write (40,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write(40,'(3ES12.5E2)') ages_ZFT(i)
end if
end do
write(40, '(a)') ''
end if
if (output_Muscovite==1) then
write (40,'(a)') 'SCALARS Muscovite FLOAT 1'
write (40,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write(40,'(3ES12.5E2)') ages_Muscovite(i)
end if
end do
write(40, '(a)') ''
end if
close (40)
end if
end do writeages
end if
out_track_any = output_velo+output_veloz+output_temp+output_press+output_time+output_istep+output_ID
if (output_mode==2 .and. trcl%ntracks>0) then
!output ages and surfaced tracks for each output timestep
do l=output_first,output_last,output_step
!count number of tracks and number of points in those for this timestep
write_tracks = 0
write_points = 0
do i=1,trcl%ntracks if (trcl%tracks(i,5)==l) then
write_tracks = write_tracks+1
write_points = write_points+trcl%tracks(i,4)
end if
end do
if (write_tracks>0) then
call int_to_char (cistep,4,l)
open (unit=30,file=trim(output_dir)//'tracks_t'//cistep//'.vtk',status='replace')
!write header
write(30,'(a)') '# vtk DataFile Version 3.1'
write(30, '(a)') 'DOUAR tracking cloud output of surfaced tracks'
write(30, '(a)') 'ASCII'
write(30, '(a)') 'DATASET UNSTRUCTURED_GRID'
write(30, '(a)') ''
!write point coordinates
write(30, '(a,i10,a)') 'POINTS ', write_points, ' FLOAT'
do i=1,trcl%ntracks
!go through all tracks
if (trcl%tracks(i,5)==l) then
!if track surfaces at this timestep:
do k=trcl%tracks(i,2),trcl%tracks(i,3)
!write all point positions within track extent
write(30,'(3ES12.5E2)') trcl%x(trcl%tracks(i,1),k), &
trcl%y(trcl%tracks(i,1),k),trcl%z(trcl%tracks(i,1),k)
end do
end if
end do
write(30, '(a)') ''
!write cell data
write(30, '(a,i10,i10)') 'CELLS ', write_tracks, write_points+write_tracks
ipoint = 0
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
!write length of track
write(30, '(i4)',advance='no') trcl%tracks(i,4)
do k=trcl%tracks(i,2),trcl%tracks(i,3)
!write point ID
write(30, '(i9)', advance='no') ipoint
!advance point ID
ipoint = ipoint+1
end do
!end of track
write (30, '(a)') ''
end if
end do
write (30, '(a)') ''
!write cell data type, 4=polyline
write(30, '(a,i10)') 'CELL_TYPES ',write_tracks
do i = 1,write_tracks
write(30, '(i2)') 4
end do
write (30, '(a)') ''
!write point datasets
if (out_track_any > 0) write (30,'(a,i10)') 'POINT_DATA ',write_points
if (output_velo==1) then
write (30,'(a)') 'SCALARS Velo FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(3ES12.5E2)') trcl%velo(trcl%tracks(i,1),k)
end do end if
end do
write(30, '(a)') ''
end if
if (output_veloz==1) then
write (30,'(a)') 'SCALARS VeloZ FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(3ES12.5E2)') trcl%veloz(trcl%tracks(i,1),k)
end do
end if
end do
write(30, '(a)') ''
end if
if (output_temp==1) then
write (30,'(a)') 'SCALARS Temperature FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(3ES12.5E2)') trcl%temp(trcl%tracks(i,1),k)
end do
end if
end do
write(30, '(a)') ''
end if
if (output_press==1) then
write (30,'(a)') 'SCALARS Pressure FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(3ES12.5E2)') trcl%press(trcl%tracks(i,1),k)
end do
end if
end do
write(30, '(a)') ''
end if
if (output_time==1) then
write (30,'(a)') 'SCALARS Time FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(3ES12.5E2)') trcl%time(k)
end do
end if
end do
write(30, '(a)') ''
end if
if (output_istep==1) then
write (30,'(a)') 'SCALARS Step INTEGER 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(i5)') trcl%step(k)
end do
end if
end do
write(30, '(a)') ''
end if
if (output_ID==1) then
write (30,'(a)') 'SCALARS Track_ID INTEGER 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
do k=trcl%tracks(i,2),trcl%tracks(i,3)
write(30,'(i5)') trcl%tracks(i,6)
end do
end if
end do
write(30, '(a)') ''
end if
close (30)
else
write (*,'(a,i4)') 'No surfaced tracks for step ', l
end if
end do
end if
if (output_mode==3) then
!output all tracks
do l=output_first,output_last,output_step
call int_to_char (cistep,4,l)
open (unit=30,file=trim(output_dir)//'alltracks_t'//cistep//'.vtk',status='replace')
!calculate part of array to write
write_nstep = (l-input_first)/input_step+1
!write header
write(30,'(a)') '# vtk DataFile Version 3.1'
write(30, '(a)') 'DOUAR tracking cloud output of all tracks'
write(30, '(a)') 'ASCII'
write(30, '(a)') 'DATASET UNSTRUCTURED_GRID'
write(30, '(a)') ''
!write point coordinates
write(30, '(a,i8,a)') 'POINTS ', trcl%np*write_nstep, ' FLOAT'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(3ES12.5E2)') trcl%x(i,k),trcl%y(i,k),trcl%z(i,k)
end do
end do
write(30, '(a)') ''
!write cell data
write(30, '(a,i10,i10)') 'CELLS ',trcl%np, trcl%np*(write_nstep+1)
do i = 0,trcl%np-1
write(30, '(i4)',advance='no') write_nstep
do k = 1,write_nstep
write(30, '(i8)', advance='no') i*write_nstep+k-1
end do
write (30, '(a)') ''
end do
write (30, '(a)') ''
!write cell data type, 4=polyline
write(30, '(a,i10)') 'CELL_TYPES ',trcl%np
do i = 1,trcl%np
write(30, '(i2)') 2
end do
write (30, '(a)') ''
!write point datasets
if (out_track_any > 0) write (30,'(a,i10)') 'POINT_DATA ',trcl%np*write_nstep
if (output_velo==1) then
write (30,'(a)') 'SCALARS Velo FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np do k=1,write_nstep
write(30,'(ES12.5E2)') trcl%velo(i,k)
end do
end do
write(30, '(a)') ''
end if
if (output_veloz==1) then
write (30,'(a)') 'SCALARS VeloZ FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(ES12.5E2)') trcl%veloz(i,k)
end do
end do
write(30, '(a)') ''
end if
if (output_temp==1) then
write (30,'(a)') 'SCALARS Temperature FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(ES12.5E2)') trcl%temp(i,k)
end do
end do
write(30, '(a)') ''
end if
if (output_press==1) then
write (30,'(a)') 'SCALARS Pressure FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(ES12.5E2)') trcl%press(i,k)
end do
end do
write(30, '(a)') ''
end if
if (output_time==1) then
write (30,'(a)') 'SCALARS Time FLOAT 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(ES12.5E2)') trcl%time(k)
end do
end do
write(30, '(a)') ''
end if
if (output_istep==1) then
write (30,'(a)') 'SCALARS Step INTEGER 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(i5)') trcl%step(k)
end do
end do
write(30, '(a)') ''
end if
if (output_ID==1) then
write (30,'(a)') 'SCALARS Particle_ID INTEGER 1'
write (30,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%np
do k=1,write_nstep
write(30,'(i6)') i
end do
end do write(30, '(a)') ''
end if
close (30)
end do
end if
deallocate (trcl%x)
deallocate (trcl%y)
deallocate (trcl%z)
deallocate (trcl%temp)
deallocate (trcl%press)
deallocate (trcl%step)
deallocate (trcl%time)
if (trcl%ntracks>0) then
deallocate (trcl%tracks)
if (output_mode==1 .or. output_mode==2) then
deallocate (ages_AHe)
deallocate (ages_ZHe)
deallocate (ages_AFT)
deallocate (ages_ZFT)
deallocate (ages_Muscovite)
end if
end if
if (output_velo==1) deallocate (trcl%velo)
if (output_veloz==1) deallocate (trcl%veloz)
end if
deallocate (allages)
call mpi_finalize (ierr)
end program post_tracking
subroutine correct_number (x)
!'fix' Nan, Inf or underflow values so that they lie in float (=real(4)) range.
! new NaN Value: 1e38
! new -inf value: -1e37
! new +inf value: 1e37
! underflow set to 0
implicit none
real(4), intent(inout) :: x
if (isnan(x)) then
x=1.0e38
else
if ( x > 1.0e37) x=1.0e37
if ( x < -1.0e37) x=-1.0e37
if (abs(x) < 1.0e-37) x=0.0
end if
end subroutine correct_number