From 89611a43682aeba770c01c8f6bca956df695505b Mon Sep 17 00:00:00 2001 From: Mickey Rush Date: Wed, 27 Nov 2024 11:36:05 -0700 Subject: [PATCH 1/4] particle exit via wells --- EcoSLIM.f90 | 103 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 95 insertions(+), 8 deletions(-) diff --git a/EcoSLIM.f90 b/EcoSLIM.f90 index 8af69e6..d61f9b1 100644 --- a/EcoSLIM.f90 +++ b/EcoSLIM.f90 @@ -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] @@ -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 @@ -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)) @@ -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 @@ -980,8 +988,8 @@ end subroutine vtk_write_points write(11,*) 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' + WellOut NP_PrecipIn NP_ETOut & + NP_WellOut NP_QOut NP_active_old NP_filtered' flush(11) !! open exited particle file and write header @@ -1355,7 +1363,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) @@ -1394,6 +1402,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 @@ -1427,12 +1437,49 @@ 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 + 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 + endif ! end ET vs. Well if statement goto 999 end if end if ! end-if for evaptrans < 0 @@ -1634,7 +1681,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 @@ -1644,6 +1691,25 @@ end subroutine vtk_write_points end if end if +!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 @@ -1816,6 +1882,27 @@ 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') From d425a1ff1eece1e164be7d629eeaa147c7e80207 Mon Sep 17 00:00:00 2001 From: Mickey Rush Date: Fri, 29 Nov 2024 10:45:58 -0700 Subject: [PATCH 2/4] making well variables shared --- EcoSLIM.f90 | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/EcoSLIM.f90 b/EcoSLIM.f90 index d61f9b1..bd7140d 100644 --- a/EcoSLIM.f90 +++ b/EcoSLIM.f90 @@ -988,7 +988,7 @@ end subroutine vtk_write_points write(11,*) write(11,*) ' **** Transient Simulation Particle Accounting ****' write(11,*) ' Timestep PFTimestep OutStep Time Mean_Age Mean_Comp Mean_Mass Total_Mass PrecipIn ETOut & - WellOut NP_PrecipIn NP_ETOut & + NP_PrecipIn NP_ETOut & NP_WellOut NP_QOut NP_active_old NP_filtered' flush(11) @@ -1263,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 @@ -1438,17 +1439,17 @@ end subroutine vtk_write_points !$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 - else + 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 @@ -1462,10 +1463,9 @@ end subroutine vtk_write_points !$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 @@ -1476,15 +1476,14 @@ end subroutine vtk_write_points !$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 - endif ! end ET vs. Well if statement - goto 999 - end if + 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 @@ -1603,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) @@ -1641,7 +1639,6 @@ end subroutine vtk_write_points call system_clock(T1) - !! format statements for particle output 61 FORMAT(4(e12.5)) 62 FORMAT(4(e12.5)) @@ -1691,6 +1688,11 @@ 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 @@ -1797,7 +1799,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 @@ -1893,7 +1895,6 @@ end subroutine vtk_write_points 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) From c1f70a8b17ccbb1c6358c0a2245fd54b627b0b91 Mon Sep 17 00:00:00 2001 From: Mickey Rush Date: Mon, 2 Dec 2024 17:41:18 -0700 Subject: [PATCH 3/4] transient output format --- EcoSLIM.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/EcoSLIM.f90 b/EcoSLIM.f90 index bd7140d..f81bd3b 100644 --- a/EcoSLIM.f90 +++ b/EcoSLIM.f90 @@ -1640,6 +1640,7 @@ end subroutine vtk_write_points call system_clock(T1) !! format statements for particle output +60 FORMAT(4(e12.10)) 61 FORMAT(4(e12.5)) 62 FORMAT(4(e12.5)) @@ -1654,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,60) P(ii,1), P(ii,2), P(ii,3), P(ii,4), P(ii,11) end do close(14) end if From 7beb8b2292fd65291b5405b42acf33e799c4eab5 Mon Sep 17 00:00:00 2001 From: Mickey Rush Date: Tue, 3 Dec 2024 09:38:18 -0700 Subject: [PATCH 4/4] fixed formatting --- EcoSLIM.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/EcoSLIM.f90 b/EcoSLIM.f90 index f81bd3b..e2d6b0a 100644 --- a/EcoSLIM.f90 +++ b/EcoSLIM.f90 @@ -1640,9 +1640,9 @@ end subroutine vtk_write_points call system_clock(T1) !! format statements for particle output -60 FORMAT(4(e12.10)) 61 FORMAT(4(e12.5)) 62 FORMAT(4(e12.5)) +63 FORMAT(4(e15.8)) write(filenumout,'(i5.5)') outkk @@ -1655,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,60) 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