Skip to content

Consolidates remaining three pathogen pull requests #66

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

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Consolidates remaining three pathogen pull requests #66

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
84 changes: 63 additions & 21 deletions src/aed_pathogens.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
!# #
!# http://aquatic.science.uwa.edu.au/ #
!# #
!# Copyright 2017 - 2022 - The University of Western Australia #
!# Copyright 2017 - 2023 - The University of Western Australia #
!# #
!# AED is free software: you can redistribute it and/or modify #
!# it under the terms of the GNU General Public License as published by #
Expand Down Expand Up @@ -101,9 +101,10 @@ MODULE aed_pathogens

! Diagnostic IDs for processes
INTEGER,ALLOCATABLE :: id_growth(:), id_mortality(:), id_sunlight(:), id_grazing(:), id_total(:)
INTEGER,ALLOCATABLE :: id_pth_f_sed(:), id_pth_d_sed(:), id_pth_a_sed(:), id_attachment(:)

INTEGER :: id_oxy, id_pH, id_doc, id_tss ! Dependency ID
INTEGER :: id_tem, id_sal ! Environmental IDs (3D)
INTEGER :: id_tem, id_sal, id_dz ! Environmental IDs (3D)
INTEGER :: id_par, id_nir, id_uva, id_uvb ! Environmental IDs (3D)
INTEGER :: id_I_0 ! Environmental ID (2D)
INTEGER :: resuspension
Expand Down Expand Up @@ -253,6 +254,7 @@ SUBROUTINE aed_define_pathogens(data, namlst)
data%id_uva = aed_locate_global('uva')
data%id_uvb = aed_locate_global('uvb')
data%id_I_0 = aed_locate_sheet_global('par_sf')
data%id_dz = aed_locate_global('layer_ht')
IF ( resuspension > 0 ) &
data%id_taub = aed_locate_sheet_global('taub')

Expand Down Expand Up @@ -379,8 +381,7 @@ SUBROUTINE aed_pathogens_load_params(data, dbase, count, list)
status = load_csv(dbase, pd)
CASE (NML_TYPE)
print*,"nml format parameter file is deprecated. Please update to CSV format"
tfil = find_free_lun()
open(tfil,file=dbase, status='OLD',iostat=status)
open(NEWUNIT=tfil,file=dbase, status='OLD',iostat=status)
IF (status /= 0) STOP 'Error opening namelist pathogen_data'
read(tfil,nml=pathogen_data,iostat=status)
close(tfil)
Expand All @@ -403,6 +404,10 @@ SUBROUTINE aed_pathogens_load_params(data, dbase, count, list)
ALLOCATE(data%id_growth(count))
ALLOCATE(data%id_sunlight(count))
ALLOCATE(data%id_mortality(count))
ALLOCATE(data%id_pth_f_sed(count))
ALLOCATE(data%id_pth_d_sed(count))
ALLOCATE(data%id_pth_a_sed(count))
ALLOCATE(data%id_attachment(count))
ENDIF

DO i=1,count
Expand Down Expand Up @@ -456,6 +461,16 @@ SUBROUTINE aed_pathogens_load_params(data, dbase, count, list)
data%id_growth(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_g', 'orgs/m3/day', 'growth')
data%id_sunlight(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_l', 'orgs/m3/day', 'sunlight')
data%id_mortality(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_m', 'orgs/m3/day', 'mortality')
data%id_pth_f_sed(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_f_set', &
'orgs/m3/d', 'alive sedimentation')
data%id_pth_d_sed(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_d_set', &
'orgs/m3/d', 'dead sedimentation')
IF (data%pathogens(i)%coef_sett_fa > zero_) THEN
data%id_pth_a_sed(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_a_set', &
'orgs/m3/d', 'attached sedimentation')
data%id_attachment(i) = aed_define_diag_variable( TRIM(data%pathogens(i)%p_name)//'_att', &
'orgs/m3/d', 'attachment rate')
END IF
ENDIF
ENDDO
DEALLOCATE(pd)
Expand Down Expand Up @@ -525,7 +540,10 @@ SUBROUTINE aed_calculate_pathogens(data,column,layer_idx)
pth_f = _STATE_VAR_(data%id_pf(pth_i))
IF (data%pathogens(pth_i)%coef_sett_fa > zero_) THEN
pth_a = _STATE_VAR_(data%id_pa(pth_i))
ELSE
pth_a = 0.0
END IF
pth_d = _STATE_VAR_(data%id_pd(pth_i))

growth = zero_
predation = zero_
Expand Down Expand Up @@ -579,21 +597,33 @@ SUBROUTINE aed_calculate_pathogens(data,column,layer_idx)


!-- Attachment of free orgs to particles (as impacted by SS and desired attachment ratio)
! 24/02/2023
! The previous formulation used _FLUX_VAR_as a test, but this is zero if resus_flux is zero
! (which is common and will apply to cells away from the bed anyway)
! This meant that attachment was also being set to zero in these cells
! The formulation below uses existing conditions and stores to work out
! attachment and detachment kinetics, assuming unlimited sediment.
! For a proper balance on the availability to meet attachment, the test should be
! attachment > pth_f/dt and
! attachment > pth_a/dt
! which is a comparison of two orgs/m3/s rates
! but dt cannot be seen in this modulle so as a stop gap I used a typical
! dt of 15 minutes. This should be fixed by a coder who knows what they are doing!

attachment = zero_
IF (data%pathogens(pth_i)%coef_sett_fa > zero_ .AND. (pth_f+pth_a) > 1e-2) THEN
! First check if ratio at last time step is less than desired (ie coef_sett_fa)
! Compute ratio at last time step for compariosn with coef_sett_fa
att_frac = pth_a/(pth_a+pth_f)
IF (att_frac<data%pathogens(pth_i)%coef_sett_fa) THEN
! Assume rate of attachment is slow based on att_ts (orgs/m3/s)
! Small fraction, so the "attachment deficit" is redistributed each time step
attachment = data%att_ts * (data%pathogens(pth_i)%coef_sett_fa*pth_a - pth_f)
IF(attachment>zero_ .AND. attachment > _FLUX_VAR_(data%id_pf(pth_i)))THEN
! proposed attachment is more than is available.
attachment = 0.5 * _FLUX_VAR_(data%id_pf(pth_i))
ELSEIF(attachment<zero_ .AND. attachment > _FLUX_VAR_(data%id_pa(pth_i)))THEN
! proposed de-attachment is more than is available.
attachment = 0.5 * _FLUX_VAR_(data%id_pa(pth_i))
ENDIF
! Compute attachment defecit/credit
! This was: attachment = data%att_ts * (data%pathogens(pth_i)%coef_sett_fa*pth_a - pth_f)
attachment = data%att_ts * (data%pathogens(pth_i)%coef_sett_fa*(pth_a + pth_f) - pth_a)
! Compare current frac with available pools
! Not enough pth_f to meet attachment need
IF ((att_frac<data%pathogens(pth_i)%coef_sett_fa).AND.(attachment > pth_f/900.0)) THEN
attachment = 0.5 * pth_f/900.0
! Not enough pth_a to meet detachment need
ELSEIF ((att_frac>data%pathogens(pth_i)%coef_sett_fa).AND.(abs(attachment) > pth_a/900.0)) THEN
attachment = 0.5 * pth_a/900.0
ENDIF
ENDIF

Expand All @@ -617,11 +647,15 @@ SUBROUTINE aed_calculate_pathogens(data,column,layer_idx)

!-----------------------------------------------------------------
! SET DIAGNOSTICS
_DIAG_VAR_(data%id_total(pth_i)) = pth_f + pth_a + pth_d ! orgs/m3/s
_DIAG_VAR_(data%id_total(pth_i)) = pth_f + pth_a + pth_d ! orgs/m3
IF ( diag_level >= 10 ) THEN
_DIAG_VAR_(data%id_growth(pth_i)) = growth*(pth_f + pth_a) ! orgs/m3/s
_DIAG_VAR_(data%id_sunlight(pth_i)) = light*pth_f + (light/2.)*pth_a ! orgs/m3/s
_DIAG_VAR_(data%id_mortality(pth_i)) = mortality*(pth_f + pth_a) ! orgs/m3/s
_DIAG_VAR_(data%id_growth(pth_i)) = growth*(pth_f + pth_a) * secs_per_day ! orgs/m3/d
_DIAG_VAR_(data%id_sunlight(pth_i)) = (light*pth_f + (light/2.)*pth_a) * secs_per_day ! orgs/m3/d
_DIAG_VAR_(data%id_mortality(pth_i)) = mortality*(pth_f + pth_a) * secs_per_day ! orgs/m3/d
! Attached
IF (data%pathogens(pth_i)%coef_sett_fa > zero_) THEN
_DIAG_VAR_(data%id_attachment(pth_i)) = attachment * secs_per_day ! orgs/m3/d
END IF
ENDIF
ENDDO
END SUBROUTINE aed_calculate_pathogens
Expand Down Expand Up @@ -761,16 +795,24 @@ SUBROUTINE aed_mobility_pathogens(data,column,layer_idx,mobility)
AED_REAL,INTENT(inout) :: mobility(:)
!
!LOCALS
AED_REAL :: temp, vel
AED_REAL :: temp, vel, dz
INTEGER :: ss_i,pth_i
!
!-------------------------------------------------------------------------------
!BEGIN
temp = _STATE_VAR_(data%id_tem)
dz = _STATE_VAR_(data%id_dz)

_DIAG_VAR_(data%id_pth_f_sed) = zero_
_DIAG_VAR_(data%id_pth_d_sed) = zero_

! First set velocity for free pathogen groups
DO pth_i=1,data%num_pathogens
mobility(data%id_pf(pth_i)) = data%pathogens(pth_i)%coef_sett_w_path
IF ( diag_level >= 10 ) THEN
_DIAG_VAR_(data%id_pth_f_sed(pth_i)) = (mobility(data%id_pf(pth_i))/dz) * _STATE_VAR_(data%id_pf(pth_i)) * secs_per_day ! orgs/m3/d !PRequest
_DIAG_VAR_(data%id_pth_d_sed(pth_i)) = (mobility(data%id_pd(pth_i))/dz) * _STATE_VAR_(data%id_pd(pth_i)) * secs_per_day ! orgs/m3/d !PRequest
END IF
ENDDO

! Compute settling rate of particles
Expand Down