c$Author: fredfox $ c$Date: 2002-04-29 16:21:48 $ c$Revision: 1.11 $ c$Source: /weru/cvs/weps/weps.src/asd/m2asd.for,v $ c NOTES: c c We will probably want to rewrite what this subroutine does into c several smaller routines for both speed reasons and potentially c modularity. c c The log(md) that we need could be moved to the c the initialization routine and access the log(md) rather than c accessing (md) and computing log(md) multiple times (SPEED SAVINGS). c c We also may want to may want to separate out the looping among c soil layers so that it can be done at a higher level (may make c code more modular - do only one thing extremely well concept). c This should be discussed as to whether this would be beneficial c in the long run. c c Tue Apr 6 14:15:48 CDT 1999 - LEW c ----------------------------------------------------------------- c This routine was simplified and recoded. c It now allows for the sieve cut sizes to lie outside the c range specified by "mnot" and "minf" by checking for this c situation and only using the sieve cuts between "mnot" and "minf" c (this only applies to the pertinent modified log-normal cases). c c Note that: c a) the sieve cut size array, "mdia" must consist of 2 or more sizes c and contain values which increase in size, c b) "mnot" must be greater than or equal to zero, c c) "mnot" must be greater than "minf" (with at least two sieve cut c sizes between them), c d) and the mass fractions, "mf" cannot be less than zero. c These conditions are NOT checked within this code. c c Note also that the return values "gmd" and "gsd" are the c geometric mean and geometric standard deviation of the c "transformed" parameters, based upon the specific "logcas" c used and NOT always the geometric mean and standard deviation c of the aggregate sizes. c ----------------------------------------------------------------- subroutine m2asd i (mf, nlay, i mnot, minf, o gmd, gsd) include 'p1werm.inc' include 'manage/asd.inc' c + + + PURPOSE + + + c This subroutine performs the inverse of subroutine asd2m. c m2asd computes the geometric mean & standard deviation for the c lognormal representation of the soil aggregate size distribution c from mf(i,j). c The routine decides which lognormal case to apply based on the c value of logcas: c c logcas = 0 --> "normal" lognormal case (mnot = 0, minf = infinity) c logcas = 1 --> "abnormal" lognormal case (mnot != 0, minf = infinity) c logcas = 2 --> "abnormal" lognormal case (mnot = 0, minf != infinity) c logcas = 3 --> "abnormal" lognormal case (mnot != 0, minf != infinity) c c + + + KEYWORDS + + + c aggregate size distribution, asd, sieves, mass fractions c c + + + ARGUMENT DECLARATIONS + + + real mf(msieve+1, mnsz) integer nlay real mnot(mnsz), minf(mnsz) real gmd(mnsz), gsd(mnsz) c c c + + + ARGUMENT DEFINITIONS + + + c mf - mass fractions of aggregates within sieve cuts c (sum of all mass fractions are expected to = 1.0) c nlay - number of soil layers used c mnot - minimum size aggregate (assumed value is known) c minf - maximum size aggregate (assumed value is known) c gmd - geometric mean diameter of aggregate size distribution c (or transformed asd for "modified" lognormal cases) c gsd - geometric standard deviation of aggregate size distribution c (or transformed asd for "modified" lognormal cases) c c + + + ACCESSED COMMON BLOCK VARIABLE DEFINITIONS + + + c c nsieve - number of sieves used c mdia - geometric mean dia. for each sieve cut c mnsize - minimum (imaginary) sieve size to use for computing c lower sieve cut geometric mean diameter c mxsize - maximum (imaginary) sieve size to use for computing c upper sieve cut geometric mean diameter c logcas - flag to represent which lognormal case to apply c c + + + PARAMETERS + + + c c + + + LOCAL VARIABLES + + + c real tmd(msieve+1) real alpha, beta real mdia_istart, mdia_istop, sdia_temp integer i, j integer k integer istart, istop c c + + + LOCAL VARIABLE DEFINITIONS + + + c c tmd - transformed md (later log(tmd)) c alpha - internal summation variable c beta - internal summation variable c i - loop variable for sieve diameters c istart - loop start variable for sieve diameters c istop - loop stop variable for sieve diameters c j - loop variable for soil layers c c + + + END SPECIFICATIONS + + + c c for each soil layer do 200 j=1,nlay c initialize accumulators alpha = 0.0 beta = 0.0 istart = 1 istop = nsieve + 1 c check if sieve cut fractions are between mnot and minf c adjust lower and upper mean diameters if mnot or minf c 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 c save value to be restored before exit mdia_istart = mdia(istart) c set size of lower sieve in bottom bin if( istart.eq.1 ) then sdia_temp = mnsize else sdia_temp = sdia(istart-1) end if c check if mnot falls within lower sieve bin if( (mnot(j).gt.sdia_temp).or.(mnot(j).lt.mnsize) ) then c recalculate lower bin mean diameter mdia(istart) = sqrt(sdia(istart)*mnot(j)) 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 c set size of upper sieve in top bin if( istop.eq.nsieve+1 ) then sdia_temp = mxsize else sdia_temp = sdia(istop) end if c save value to be restored before exit mdia_istop = mdia(istop) c check if minf falls within upper sieve bin or outside mxsize if( (minf(j).lt.sdia_temp).or.(minf(j).gt.mxsize) ) then c recalculate upper bin mean diameter mdia(istop) = sqrt(sdia(istop-1)*minf(j)) end if else istop = nsieve + 1 end if c 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 c now compute the log of the gmd dia tmd(i) = log(tmd(i)) c 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 c 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 c 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 c c$Log: not supported by cvs2svn $ cRevision 1.9 2001/09/27 20:36:52 fredfox cmerged in update hydro branch c cRevision 1.8.2.1 2001/08/15 22:12:33 fredfox cCorrected index for zeroing out unused array elements, added headers, edited out debug write statements. c cRevision 1.8 2000/09/08 23:54:11 fredfox cfixed test for valid MNOT for correct case and loop direction c cRevision 1.7 2000/09/08 22:52:04 wagner cAdded check for single 365 day year in cligen file and added code to keep out of an infinite loop when WEPS thinks a simulation year should have a leap day in it c cRevision 1.6 2000/09/08 15:07:20 fredfox ccorrected alpha,beta, i, j declaration statements c cRevision 1.5 2000/09/07 20:25:36 jt cadded file name to print alpha and beta. c cRevision 1.4 1999/04/22 19:02:30 wjr cdebugging write line added c cRevision 1.3 1999/04/06 19:40:29 wagner cSimplified and recoded this routine. cIt now checks for and allows sieve fractions to be outside cthe range from "mnot" to "minf" for the "modified" log-normal ccases. c cRevision 1.1.1.1 1999/03/12 17:05:17 wagner cBaseline version of WEPS with Bill Rust's modifications c c Revision 1.2 1995/09/13 15:49:44 wagner c Necessary changes made to allow FORTRAN src files (*.for) to use the c extended FORTRAN include statement rather than the MICROSOFT $INCLUDE c directive as previously used. This is required to allow use of other c FORTRAN compilers. c c Changes have been made to the prologue.mk, epilogue.mk, and the Unix c master startup.mk files as well as the src files. c c Revision 1.1.1.1 1995/01/18 04:19:56 wagner c Initial checkin c c Revision 1.4 1994/09/01 22:18:54 jt c checking for floating point errors? - LEW c c Revision 1.3 1992/10/10 21:44:14 wagner c Changed names appropriate for submodel name change c from TILLAGE to MANAGEMENT. c c Revision 1.2 1992/04/16 21:41:37 wagner c Uses common memory now. c c Revision 1.1 1992/04/16 13:29:01 wagner c Initial revision c