diff --git a/src/aed_pathogens.F90 b/src/aed_pathogens.F90 index def536e..4de2a8b 100644 --- a/src/aed_pathogens.F90 +++ b/src/aed_pathogens.F90 @@ -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 # @@ -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 @@ -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') @@ -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) @@ -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 @@ -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) @@ -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_ @@ -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_fraczero_ .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 _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 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 @@ -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 @@ -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