Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Neptune #11

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Neptune #11

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
89611a4
particle exit via wells
Nov 27, 2024
d425a1f
making well variables shared
Nov 29, 2024
c1f70a8
transient output format
Dec 3, 2024
7beb8b2
fixed formatting
Dec 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 106 additions & 17 deletions EcoSLIM.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ program EcoSLIM
! P(np,7) = Particle source (1=IC, 2=rain, 3=snowmelt...)
! P(np,8) = Particle Status (1=active, 0=inactive)
! P(np,9) = concentration
! P(np,10) = Exit status (1=surface outflow, 2=ET, 3=subsurface outflow)
! P(np,10) = Exit status (1=surface outflow, 2=ET, 3=subsurface outflow, 4=well outflow)
! P(np,11) = Particle Number (This is a unique integer identifier for the particle)
! P(np,12) = Partical Initial X coordinate [L]
! P(np,13) = Partical Initial Y coordinate [L]
Expand Down Expand Up @@ -281,7 +281,8 @@ program EcoSLIM
real*8, allocatable::ET_age(:,:), ET_comp(:,:), water_balance(:,:), ET_mass(:)
real*8, allocatable::Surf_age(:,:), Surf_mass(:,:), Surf_comp(:,:)
real*8, allocatable::Subsurf_age(:,:), Subsurf_mass(:,:), Subsurf_comp(:,:)
integer, allocatable:: ET_np(:), Surf_np(:), Subsurf_np(:)
real*8, allocatable::Well_age(:,:), Well_mass(:), Well_comp(:,:)
integer, allocatable:: ET_np(:), Surf_np(:), Subsurf_np(:), Well_np(:)


real*8 ET_dt, DR_Temp
Expand Down Expand Up @@ -454,7 +455,8 @@ end subroutine vtk_write_points

! number of things written to C array, hard wired at 2 now for Testing
!n_constituents = 5
n_constituents = 9
!n_constituents = 9
n_constituents = 13 ! additional 4 for Well
!allocate arrays
allocate(PInLoc(np,3))
allocate(Sx(nx,ny),Sy(nx,ny), DEM(nx,ny))
Expand Down Expand Up @@ -526,6 +528,12 @@ end subroutine vtk_write_points
Subsurf_comp = 0.0d0
Subsurf_np = 0

allocate(Well_age(pfnt,5), Well_comp(pfnt,3), Well_mass(pfnt), Well_np(pfnt))
Well_age = 0.0d0
Well_mass = 0.0d0
Well_comp = 0.0d0
Well_np = 0

! clear out output particles
npout = 0

Expand Down Expand Up @@ -981,7 +989,7 @@ end subroutine vtk_write_points
write(11,*) ' **** Transient Simulation Particle Accounting ****'
write(11,*) ' Timestep PFTimestep OutStep Time Mean_Age Mean_Comp Mean_Mass Total_Mass PrecipIn ETOut &
NP_PrecipIn NP_ETOut &
NP_QOut NP_active_old NP_filtered'
NP_WellOut NP_QOut NP_active_old NP_filtered'
flush(11)

!! open exited particle file and write header
Expand Down Expand Up @@ -1255,6 +1263,7 @@ end subroutine vtk_write_points
!$OMP& SHARED(kk, pfnt, surf_age, surf_mass, surf_comp, surf_np, dtfrac, et_age, et_mass) &
!$OMP& SHARED(et_comp, et_np, moldiff, efract, C,ipwrite) &
!$OMP& SHARED(subsurf_age, subsurf_mass, subsurf_comp, subsurf_np, clmtrans, velfile) &
!$OMP& SHARED(well_age, well_mass, well_comp, well_np) &
!$OMP& Default(private)

! loop over active particles
Expand Down Expand Up @@ -1355,7 +1364,7 @@ end subroutine vtk_write_points
! check each direction independently
advdt = pfdt
if (Vpx /= 0.0d0) advdt(1) = dabs(dtfrac*(dx/Vpx))
if (Vpy /= 0.0d0) advdt(2) = dabs(dtfrac*(dx/Vpy))
if (Vpy /= 0.0d0) advdt(2) = dabs(dtfrac*(dy/Vpy))
if (Vpz /= 0.0d0) advdt(3) = dtfrac*(dz(Ploc(3)+1)/dabs(Vpz))
!if (Vpx > 0.0d0) advdt(1) = dabs(((1.0d0-Clocx)*dx)/Vpx)
!if (Vpx < 0.0d0) advdt(1) = dabs((Clocx*dx)/Vpx)
Expand Down Expand Up @@ -1394,6 +1403,8 @@ end subroutine vtk_write_points
if (itime_loc >= pfnt) itime_loc = pfnt
Zr = ran1(ir)
if (Zr < ((et_flux*particledt)/water_vol)) then ! check if particle is 'captured' by the roots
! determine whether flux is ET or well based on depth / layers
if (P(ii,3) >= (Zmax - sum(dz(nz-4:nz)))) then ! particle is in top four layers
! this section made atomic since it could inovlve a data race
! that is, each thread can only update the ET arrays one at a time
!$OMP ATOMIC
Expand Down Expand Up @@ -1427,17 +1438,52 @@ end subroutine vtk_write_points
C(8,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) = C(8,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) + P(ii,4)*P(ii,6) ! mass weighted age
!$OMP Atomic
C(9,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) = C(9,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) + P(ii,7)*P(ii,6) ! mass weighted contribution

! now remove particle from domain
P(ii,8) = 0.0d0
!flag as exiting via ET
P(ii,10) = 2.0d0

goto 999
P(ii,8) = 0.0d0
!flag as exiting via ET
P(ii,10) = 2.0d0
goto 999
else
! this section made atomic since it could inovlve a data race
! that is, each thread can only update the ET arrays one at a time
!$OMP ATOMIC
Well_age(itime_loc,1) = Well_age(itime_loc,1) + P(ii,4)*P(ii,6) ! mass weighted age
!$OMP ATOMIC
Well_mass(itime_loc) = Well_mass(itime_loc) + P(ii,6) ! particle mass added to ET
!ET_comp(itime_loc,1) = ET_comp(itime_loc,1) + P(ii,7)*P(ii,6) !mass weighted contribution
if (P(ii,7) == 1.0) then
!$OMP ATOMIC
Well_comp(itime_loc,1) = Well_comp(itime_loc,1) + P(ii,6)
end if
if (P(ii,7) == 2.0) then
!$OMP ATOMIC
Well_comp(itime_loc,2) = Well_comp(itime_loc,2) + P(ii,6)
end if
if (P(ii,7) == 3.0) then
!$OMP ATOMIC
Well_comp(itime_loc,3) = Well_comp(itime_loc,3) + P(ii,6)
end if
!$OMP ATOMIC
Well_np(itime_loc) = Well_np(itime_loc) + 1 ! track number of particles

!outputting spatially distributed Well information
!$OMP Atomic
C(10,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) = C(10,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) + 1 ! Number of Well particles
!$OMP Atomic
C(11,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) = C(11,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) + P(ii,6) ! particle mass added to Well
!$OMP Atomic
C(12,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) = C(12,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) + P(ii,4)*P(ii,6) ! mass weighted age
!$OMP Atomic
C(13,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) = C(13,Ploc(1)+1,Ploc(2)+1,Ploc(3)+1) + P(ii,7)*P(ii,6) ! mass weighted contribution
! now remove particle from domain
P(ii,8) = 0.0d0
!flag as exiting via Well
P(ii,10) = 4.0d0
goto 999
end if ! end-if for ET vs. Well
end if ! end-if for captured
end if ! end-if for evaptrans < 0
end if ! end-if for clmtrans

! Advect particle to new location using Euler advection until next time
P(ii,1) = P(ii,1) + particledt * Vpx
P(ii,2) = P(ii,2) + particledt * Vpy
Expand Down Expand Up @@ -1556,7 +1602,6 @@ end subroutine vtk_write_points

end do ! end of do-while loop for particle time to next time
999 continue ! where we go if the particle is out of bounds

!! concentration routine
! Find the "adjacent" "cell corresponding to the particle's location
Ploc(1) = floor(P(ii,1) / dx)
Expand Down Expand Up @@ -1594,10 +1639,10 @@ end subroutine vtk_write_points

call system_clock(T1)


!! format statements for particle output
61 FORMAT(4(e12.5))
62 FORMAT(4(e12.5))
63 FORMAT(4(e15.8))

write(filenumout,'(i5.5)') outkk

Expand All @@ -1610,7 +1655,7 @@ end subroutine vtk_write_points
open(14,file=trim(runname)//'_transient_particle.'//trim(adjustl(filenumout))//'.3D')
write(14,*) 'X Y Z TIME ID'
do ii = 1, np_active
if (P(ii,8) == 1) write(14,61) P(ii,1), P(ii,2), P(ii,3), P(ii,4), P(ii,11)
if (P(ii,8) == 1) write(14,63) P(ii,1), P(ii,2), P(ii,3), P(ii,4), P(ii,11)
end do
close(14)
end if
Expand All @@ -1634,7 +1679,7 @@ end subroutine vtk_write_points
do i = 1, nx
do j = 1, ny
do k = 1, nz
if (EvapTrans(i,j,k) < 0.0d0) &
if ((EvapTrans(i,j,k) < 0.0d0) .and. (k >= (nz - 4))) &
write(14,'(3(i6), 7(e13.5))') i, j, k, C(6,i,j,k), C(7,i,j,k), C(8,i,j,k), &
C(9,i,j,k), EvapTrans(i,j,k), Saturation(i,j,k), Porosity(i,j,k)
end do
Expand All @@ -1644,6 +1689,30 @@ end subroutine vtk_write_points
end if
end if

! normalize Well ages by mass
where (C(11,:,:,:)>0.0) C(12,:,:,:) = C(12,:,:,:) / C(11,:,:,:)
where (C(11,:,:,:)>0.0) C(13,:,:,:) = C(13,:,:,:) / C(11,:,:,:)


!Write gridded Well outputs to text files
if(etwrite > 0) then
if (mod(kk,etwrite) == 0) then
! open/create/write the 3D output file
open(14,file=trim(runname)//'_Well_summary.'//trim(adjustl(filenumout))//'.txt')
write(14,*) 'X Y Z Well_npart, Well_mass, Well_age, Well_comp, Well_Rate, Saturation, Porosity'
do i = 1, nx
do j = 1, ny
do k = 1, nz
if ((EvapTrans(i,j,k) < 0.0d0) .and. (k < (nz - 4))) &
write(14,'(3(i6), 7(e13.5))') i, j, k, C(10,i,j,k), C(11,i,j,k), C(12,i,j,k), &
C(13,i,j,k), EvapTrans(i,j,k), Saturation(i,j,k), Porosity(i,j,k)
end do
end do
end do
close(14)
end if
end if

! write grid based values ("concentrations")
vtk_file=trim(runname)//'_cgrid'
if(icwrite > 0) then
Expand Down Expand Up @@ -1731,7 +1800,7 @@ end subroutine vtk_write_points
write(11,'(3(i10),3(f12.5),4(1x,e12.5,1x),3(i8),2(i12))') kk, pfkk, outkk, Time_Next(kk), mean_age , mean_comp, mean_mass, &
total_mass, water_balance(kk,1), water_balance(kk,2), &
i_added_particles, &
ET_np(kk), Surf_np(kk), np_active,np_active2
ET_np(kk), Well_np(kk), Surf_np(kk), np_active,np_active2
flush(11)

np_active = np_active2
Expand Down Expand Up @@ -1816,6 +1885,26 @@ end subroutine vtk_write_points
! close ET file
close(13)

!! write Well files
!
open(13,file=trim(runname)//'_Well_output.txt')
write(13,*) 'TIME Well_age Well_comp1 Well_comp2 Well_comp3 Well_mass Well_Np'
do ii = 1, pfnt
if (Well_mass(ii) > 0 ) then
Well_age(ii,1) = Well_age(ii,1)/(Well_mass(ii))
Well_comp(ii,1) = Well_comp(ii,1)/(Well_mass(ii))
Well_comp(ii,2) = Well_comp(ii,2)/(Well_mass(ii))
Well_comp(ii,3) = Well_comp(ii,3)/(Well_mass(ii))
end if
write(13,'(6(e12.5),i12)') float(ii+tout1-1)*ET_dt, Well_age(ii,1), Well_comp(ii,1), &
Well_comp(ii,2), Well_comp(ii,3), Well_mass(ii), Well_np(ii)

end do
flush(13)
! close Well file
close(13)


!! write surface outflow
!
open(13,file=trim(runname)//'_surface_outflow.txt')
Expand Down