!$Author$ !$Date$ !$Revision$ !$HeadURL$ module asd_mod integer, parameter :: msieve = 26 real, parameter :: mingsd = 2.0 real sdia(msieve) ! array containing sieve size diameters real mnsize ! minimum (imaginary) sieve size to use for computing ! lower sieve cut geometric mean diameter real mxsize ! maximum (imaginary) sieve size to use for computing ! upper sieve cut geometric mean diameter real mdia(msieve+1) ! array containing gmd sieve cut diameters integer nsieve ! number of sieves used integer logcas ! flag to represent which lognormal case to apply contains subroutine asdini() ! + + + PURPOSE + + + ! This subroutine performs the initialization of the asd/sieve ! variables which include the number of sieves and their sizes, ! the geometric mean diameter of each sieve cut and specifies which ! lognormal case will be used to represent aggregate size distributions ! in WERM/WEPS. ! The routine decides which lognormal case to apply based on the ! value of logcas: ! logcas = 0 --> "normal" lognormal case (mnot = 0, minf = infinity) ! logcas = 1 --> "abnormal" lognormal case (mnot != 0, minf = infinity) ! logcas = 2 --> "abnormal" lognormal case (mnot = 0, minf != infinity) ! logcas = 3 --> "abnormal" lognormal case (mnot != 0, minf != infinity) ! + + + KEYWORDS + + + ! aggregate size distribution, asd, sieves, mass fractions ! + + + LOCAL VARIABLES + + + integer :: i ! loop variable for sieve diameters ! + + + END SPECIFICATIONS + + + ! NOTE: using this method generates slightly different sieve sizes ! between debug and optimized compile switches. (and possibly between ! different compilers) To minimize these differences, we should return ! to exactly defined sieve sizes ! specificiations brought here from BLKDAT (see revision 1.4 comment below) ! data logcas / 3 / ! data nsieve / 13 / ! data sdia / 0.018, 0.037, 0.075, 0.15, 0.42, 0.84, 2.0, ! & 6.35, 19.05, 44.45, 76.2, 150.4, 300.8 / ! data mnsize, mxsize / 0.009, 601.2 / logcas = 3 nsieve = msieve - 1 mnsize = 0.005 mxsize = 1000.0 do i = 1, nsieve sdia(i) = exp(log(mnsize) & & + i*(log(mxsize)-log(mnsize))/(nsieve+1)) end do ! compute geometric mean dia. for each sieve cut mdia(1) = sqrt(mnsize*sdia(1)) do 5 i = 2, nsieve mdia(i) = sqrt(sdia(i)*sdia(i-1)) 5 continue mdia(nsieve+1) = sqrt(mxsize*sdia(nsieve)) return end subroutine asdini subroutine asd2m (mnot, minf, gmd, gsd, nlay, mf) ! + + + PURPOSE + + + ! This subroutine performs the inverse of subroutine m2asd. ! asd2m computes the mass fractions for each sieve cut from the ! lognormal representation of the soil aggregate size distribution. ! The routine decides which lognormal case to apply based on the ! value of logcas: ! logcas = 0 --> "normal" lognormal case (mnot = 0, minf = infinity) ! logcas = 1 --> "abnormal" lognormal case (mnot != 0, minf = infinity) ! logcas = 2 --> "abnormal" lognormal case (mnot = 0, minf != infinity) ! logcas = 3 --> "abnormal" lognormal case (mnot != 0, minf != infinity) ! ! + + + KEYWORDS + + + ! aggregate size distribution, asd, sieves, mass fractions ! + + + ARGUMENT DECLARATIONS + + + real mnot(*), minf(*) real gmd(*), gsd(*) integer nlay real, dimension(msieve+1,*) :: mf ! + + + ARGUMENT DEFINITIONS + + + ! mnot - minimum size aggregate (assumed value is known) ! minf - maximum size aggregate (assumed value is known) ! gmd - geometric mean diameter of aggregate size distribution ! (or transformed asd for "modified" lognormal cases) ! gsd - geometric standard deviation of aggregate size distribution ! (or transformed asd for "modified" lognormal cases) ! nlay - number of soil layers used ! mf - mass fractions of aggregates within sieve cuts ! (sum of all mass fractions are expected to = 1.0) ! + + + LOCAL VARIABLES + + + real d(msieve+1) real lngmd, lngsd real prev, this integer i, j ! + + + FUNCTION DEFINITIONS + + + real erf ! + + + LOCAL VARIABLE DEFINITIONS + + + ! ! d - transformed sieve dia. values ! (if "abnormal" lognormal cases) ! lngmd - natural log of gmd ! lngsd - natural log of gsd ! prev - contain previous sieve dia. cumulative prob ! this - contain this sieve dia. cumulative prob ! i - loop variable for sieve sizes ! j - loop variable for soil layers ! ! + + + END SPECIFICATIONS + + + do 20 j = 1, nlay ! compute transformed sieve dia. sizes if (logcas .eq. 3) then do 1 i = 1, nsieve if(sdia(i).lt.minf(j)) then d(i) = (sdia(i)-mnot(j))*(minf(j)-mnot(j))/ & & (minf(j)-sdia(i)) end if 1 continue elseif (logcas .eq. 2) then do 2 i = 1, nsieve if(sdia(i).lt.minf(j)) then d(i) = sdia(i)*minf(j)/(minf(j)-sdia(i)) end if 2 continue elseif (logcas .eq. 1) then do 3 i = 1, nsieve d(i) = sdia(i)-mnot(j) 3 continue elseif (logcas .eq. 0) then do 4 i = 1, nsieve d(i) = sdia(i) 4 continue endif lngmd= log(gmd(j)) lngsd= sqrt(2.0) * log(max(mingsd,gsd(j))) prev= 1.0 ! compute each dia. cumulative probability do 10 i = 1, nsieve if (sdia(i) .le. mnot(j)) then this = 1.0 else if (sdia(i) .lt. minf(j)) then this = 0.5 -0.5*erf((alog(d(i)) - lngmd) / lngsd) else this = 0.0 end if ! compute mass fraction between prev and this dia mf(i,j) = prev - this prev = this ! write(*,*) 'asd2m:',i,sdia(i),this,mf(i,j) ! if roundoff errors or otherwise results in negative ! mass fraction then set to zero mass if (mf(i,j) .lt. 0.0) then mf(i,j) = 0.0 else prev = this endif ! if(j.eq.4) write(*,*) 'asd2m: mf(',i,j,')',mf(i,j) 10 continue ! get mass fraction for upper-most sieve cut mf(nsieve+1,j) = prev ! if( j.eq.1 )write(*,*)'asd2m: mf(',nsieve+1,j,')',mf(nsieve+1,j) ! zero out the rest of the array which is used every where else do i=nsieve+2, msieve+1 mf(i,j) = 0.0 end do 20 continue return end subroutine asd2m subroutine m2asd (mf, nlay, mnot, minf, gmd, gsd) ! + + + PURPOSE + + + ! This subroutine performs the inverse of subroutine asd2m. ! m2asd computes the geometric mean & standard deviation for the ! lognormal representation of the soil aggregate size distribution ! from mf(i,j). ! The routine decides which lognormal case to apply based on the ! value of logcas: ! logcas = 0 --> "normal" lognormal case (mnot = 0, minf = infinity) ! logcas = 1 --> "abnormal" lognormal case (mnot != 0, minf = infinity) ! logcas = 2 --> "abnormal" lognormal case (mnot = 0, minf != infinity) ! logcas = 3 --> "abnormal" lognormal case (mnot != 0, minf != infinity) ! + + + KEYWORDS + + + ! aggregate size distribution, asd, sieves, mass fractions ! + + + ARGUMENT DECLARATIONS + + + real, dimension(msieve+1,*) :: mf integer nlay real mnot(*), minf(*) real gmd(*), gsd(*) ! + + + ARGUMENT DEFINITIONS + + + ! mf - mass fractions of aggregates within sieve cuts ! (sum of all mass fractions are expected to = 1.0) ! nlay - number of soil layers used ! mnot - minimum size aggregate (assumed value is known) ! minf - maximum size aggregate (assumed value is known) ! gmd - geometric mean diameter of aggregate size distribution ! (or transformed asd for "modified" lognormal cases) ! gsd - geometric standard deviation of aggregate size distribution ! (or transformed asd for "modified" lognormal cases) ! + + + LOCAL VARIABLES + + + real tmd(msieve+1) real alpha, beta real mdia_istart, mdia_istop, sdia_temp integer i, j integer istart, istop ! + + + LOCAL VARIABLE DEFINITIONS + + + ! ! tmd - transformed md (later log(tmd)) ! alpha - internal summation variable ! beta - internal summation variable ! i - loop variable for sieve diameters ! istart - loop start variable for sieve diameters ! istop - loop stop variable for sieve diameters ! j - loop variable for soil layers ! + + + END SPECIFICATIONS + + + ! for each soil layer do 200 j=1,nlay ! initialize accumulators alpha = 0.0 beta = 0.0 istart = 1 istop = nsieve + 1 ! check if sieve cut fractions are between mnot and minf ! adjust lower and upper mean diameters if mnot or minf ! fall within the bin range if (logcas .eq. 1 .or. logcas .eq. 3) then do i=nsieve, 1, -1 if (sdia(i) .gt. mnot(j)) then istart = i end if end do ! save value to be restored before exit mdia_istart = mdia(istart) ! set size of lower sieve in bottom bin if( istart.eq.1 ) then sdia_temp = mnsize else sdia_temp = sdia(istart-1) end if ! check if mnot falls within lower sieve bin if( (mnot(j).gt.sdia_temp).or.(mnot(j).lt.mnsize) ) then ! recalculate lower bin mean diameter mdia(istart) = sqrt(sdia(istart)*mnot(j)) if (logcas .eq. 3) then ! check that mdia is greater than mnot, or method fails if( mdia(istart) .lt. mnot(j) * 1.00001 ) then mdia(istart) = mnot(j) * 1.00001 end if end if end if endif if (logcas .ge. 2) then do i=1, nsieve if (sdia(i) .le. minf(j)) then istop = i+1 end if end do ! set size of upper sieve in top bin if( istop.eq.nsieve+1 ) then sdia_temp = mxsize else sdia_temp = sdia(istop) end if ! save value to be restored before exit mdia_istop = mdia(istop) ! check if minf falls within upper sieve bin or outside mxsize if( (minf(j).lt.sdia_temp).or.(minf(j).gt.mxsize) ) then ! recalculate upper bin mean diameter mdia(istop) = sqrt(sdia(istop-1)*minf(j)) if (logcas .ge. 2) then ! check that mdia is less than minf, or method fails if( mdia(istop) .gt. minf(j) * 0.99999 ) then mdia(istop) = minf(j) * 0.99999 end if end if end if else istop = nsieve + 1 end if ! do transformations for "modified" log-normal cases do i= istart, istop if (logcas .eq. 3) then tmd(i) = (mdia(i)-mnot(j))*(minf(j)-mnot(j))/ & & (minf(j)-mdia(i)) else if (logcas .eq. 2) then tmd(i) = mdia(i)*minf(j)/(minf(j)-mdia(i)) else if (logcas .eq. 1) then tmd(i) = mdia(i)-mnot(j) else tmd(i) = mdia(i) end if ! now compute the log of the gmd dia tmd(i) = log(tmd(i)) ! sum diameters & their squares, over all aggregate sizes alpha = alpha + (mf(i,j)*tmd(i)) beta = beta + (mf(i,j)*tmd(i)*tmd(i)) end do ! compute geometric mean and standard deviation gmd(j) = exp(alpha) if( beta-alpha*alpha.le.0.0 ) then gsd(j) = mingsd else gsd(j) = max(mingsd,exp(sqrt(beta-alpha*alpha))) end if ! restore modified geometric mean bin diameters if (logcas .eq. 1 .or. logcas .eq. 3) then mdia(istart) = mdia_istart end if if (logcas .ge. 2) then mdia(istop) = mdia_istop end if 200 continue return end subroutine m2asd end module asd_mod