Skip to content
Snippets Groups Projects
Commit 207fe367 authored by Matthias Schmiddunser's avatar Matthias Schmiddunser
Browse files

Self-scheduling age calculation 

Thermochron Ages are calculated in a self-scheduling manner. Write
output moved to subroutine. Output only in *.vtk
parent 06f2e5c4
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
src/Makefile
+ 3
− 0
View file @ 207fe367
......@@ -128,6 +128,7 @@ subdirs:
cd NN; make all
cd OCTREE; make all
cd RESAMPLE; make all
cd TRPOST; make all
vtk:
cd OCTREE; make all
......@@ -139,6 +140,7 @@ objclean:
cd NN; make $@
cd OCTREE; make $@
cd RESAMPLE; make $@
cd TRPOST; make $@
rm -f *.o *.mod $(COMPILE_OUT)
.PHONY: distclean
......@@ -147,6 +149,7 @@ distclean: objclean
cd NN; make $@
cd OCTREE; make $@
cd RESAMPLE; make $@
cd TRPOST; make $@
rm -f $(BIN) $(BIN).link.stdout $(BIN).link.stderr module_gitversion.f90
.PHONY: more-output
......
src/Makefile.inc.esd.argand
+ 2
− 2
View file @ 207fe367
......@@ -22,10 +22,10 @@ INCLUDE = $(MPI_INC)
# Define compiler flags:
#F90FLAGS = -O3 -Wall -Werror -Warray-bounds -Wunderflow -fbacktrace -fdump-core -ffree-line-length-none -fcheck=all
F90FLAGS = -O2 -Wall -Warray-bounds -Wunderflow -fbacktrace -fdump-core -ffree-line-length-none -fcheck=all
F90FLAGS = -O3 -Wall -Warray-bounds -Wunderflow -fbacktrace -fdump-core -ffree-line-length-none -fcheck=all
F77FLAGS = $(F90FLAGS)
# CFLAGS = -O3 -Wall -Werror
CFLAGS = -O2 -Wall
CFLAGS = -O3 -Wall
WSMP = -L ../../wsmp/wsmp-Linux64-GNU/lib/openmpi/ -lpwsmp64
......
src/TRPOST/Makefile
+ 13
− 5
View file @ 207fe367
......@@ -2,24 +2,32 @@
NAME=post_track
include ../Makefile.inc
INCLUDES=
INCLUDE=-I.. $(MPI_INC)
BIN=$(NAME)$(BITS)
# object (make VTK output for visualization)
OBJECTS = \
../module_definitions.o \
../toolbox.o \
post_tracking.o
tridag.o \
ageCalculation.o \
correct_number.o \
module_trdefinitions.o \
post_tracking.o \
write_post_track.o
OCTREE = ../OCTREE/libOctree$(BITS).a
LIBS = \
-L../OCTREE -lOctree$(BITS) \
$(MPI_LIB)
# make the program
#
all: $(BIN)
$(BIN): $(OBJECTS)
@echo "$(F90) $(FFLAGS) $(OPTIONS) $(OBJECTS) $(MPI_LIB) -o $(BIN)" \
@echo "$(F90) $(FFLAGS) $(INCLUDE) $(OPTIONS) $(OBJECTS) $(LIBS) -o $(BIN)" \
1>>$(BIN).link.stdout 2>>$(BIN).link.stderr
$(F90) $(FFLAGS) $(OPTIONS) $(OBJECTS) $(OCTREE) $(MPI_LIB) -o $(BIN) \
$(F90) $(FFLAGS) $(INCLUDE) $(OPTIONS) $(OBJECTS) $(LIBS) -o $(BIN) \
1>>$(BIN).link.stdout 2>>$(BIN).link.stderr
.PHONY: objclean
......
src/TRPOST/correct_number.f90 0 → 100644
+ 41
− 0
View file @ 207fe367
subroutine correct_numberdp (x)
!'fix' Nan, Inf or underflow values so that they lie in real range.
! new NaN Value: 1e38
! new -inf value: -1e37
! new +inf value: 1e37
! underflow set to 0
implicit none
real(8), 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_numberdp
subroutine correct_numbersp (x)
!'fix' Nan, Inf or underflow values so that they lie in real range.
! new NaN Value: 1e38
! new -inf value: -1e37
! new +inf value: 1e37
! underflow set to 0
implicit none
real, 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_numbersp
src/TRPOST/module_trdefinitions.f90 0 → 100644
+ 57
− 0
View file @ 207fe367
module trdefinitions
type tracking_cloud
integer np,ntracks
real, dimension(:), pointer :: time
real,dimension(:),pointer :: AHe,ZHe,AFT,ZFT,Muscovite
real, dimension(:,:),pointer :: x,y,z,press,temp,velo,veloZ
integer,dimension(:),pointer :: step
integer,dimension(:,:),pointer :: tracks
end type tracking_cloud
type input_parameters
character(150) :: dir
integer :: first
integer :: last
integer :: step
integer :: nstep
end type input_parameters
type output_parameters
character(150) :: dir
integer :: mode
integer :: scaled
integer :: first
integer :: last
integer :: step
integer :: nstep
integer :: keep_nstep
integer :: min_track
real :: pad_age
integer :: AHe
integer :: ZHe
integer :: AFT
integer :: ZFT
integer :: Muscovite
integer :: ageID
integer :: nages
integer :: velo
integer :: veloz
integer :: temp
integer :: press
integer :: time
integer :: istep
integer :: trackID
integer :: ntrackvars
end type output_parameters
interface correct_number
subroutine correct_numberdp (x)
real(8), intent(inout) :: x
end subroutine correct_numberdp
subroutine correct_numbersp (x)
real, intent(inout) :: x
end subroutine correct_numbersp
end interface correct_number
end module trdefinitions
src/TRPOST/post_tracking.f90
+ 405
− 759
View file @ 207fe367
This diff is collapsed.
src/TRPOST/write_post_track.f90 0 → 100644
+ 449
− 0
View file @ 207fe367
subroutine write_post_track (trcl,inparam,outparam)
use trdefinitions
implicit none
include 'mpif.h'
type (tracking_cloud) :: trcl
type (input_parameters) :: inparam
type (output_parameters) :: outparam
integer :: write_nstep,write_points,write_tracks
integer :: i,k,l,ipoint
character(4) :: cistep
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
if ((outparam%mode==1.or.outparam%mode==2) .and. trcl%ntracks>0) then
writeages: do l=outparam%first,outparam%last,outparam%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(outparam%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',outparam%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 (outparam%nages > 0) write (40,'(a,i10)') 'POINT_DATA ',write_tracks
if (outparam%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)') trcl%AHe(i)
end if
end do
write(40, '(a)') ''
end if
if (outparam%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)') trcl%ZHe(i)
end if
end do
write(40, '(a)') ''
end if
if (outparam%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)') trcl%AFT(i)
end if
end do
write(40, '(a)') ''
end if
if (outparam%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)') trcl%ZFT(i)
end if
end do
write(40, '(a)') ''
end if
if (outparam%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)') trcl%Muscovite(i)
end if
end do
write(40, '(a)') ''
end if
if (outparam%ageID==1) then
write (40,'(a)') 'SCALARS trackID FLOAT 1'
write (40,'(a)') 'LOOKUP_TABLE default'
do i=1,trcl%ntracks
if (trcl%tracks(i,5)==l) then
write(40,'(i5)') trcl%tracks(i,7)
end if
end do
write(40, '(a)') ''
end if
close (40)
end if
end do writeages
end if
outparam%ntrackvars = outparam%velo+outparam%veloz+outparam%temp+outparam%press+ &
outparam%time+outparam%istep+outparam%trackID
if (outparam%mode==2 .and. trcl%ntracks>0) then
!output ages and surfaced tracks for each output timestep
do l=outparam%first,outparam%last,outparam%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(outparam%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 (outparam%ntrackvars > 0) write (30,'(a,i10)') 'POINT_DATA ',write_points
if (outparam%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 (outparam%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 (outparam%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 (outparam%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 (outparam%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 (outparam%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 (outparam%trackID==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,7)
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 (outparam%mode==3) then
!output all tracks
do l=outparam%first,outparam%last,outparam%step
call int_to_char (cistep,4,l)
open (unit=30,file=trim(outparam%dir)//'alltracks_t'//cistep//'.vtk',status='replace')
!calculate part of array to write
write_nstep = (l-inparam%first)/inparam%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 (outparam%ntrackvars > 0) write (30,'(a,i10)') 'POINT_DATA ',trcl%np*write_nstep
if (outparam%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 (outparam%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 (outparam%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 (outparam%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 (outparam%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 (outparam%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 (outparam%trackID==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
end subroutine write_post_track
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment