c file: 'growth.for' subroutine growth (bnslay, bszlyd, bczrtd, bcmrtz, bc0ck, bcrhi, * bcehu0, bczmxc, bc0idc, & bczmrt, bctmin, bctopt,bc0be1, & bc0alf,bc0blf,bc0clf,bc0dlf,bc0arp, * bc0brp,bc0crp,bc0drp, & bc0aht,bc0bht,bc0ssa,bc0ssb,bc0sla,bc0hue, * bwtdav, bwtdmn, bweirr, bhfwsf) c Author : Amare Retta c + + + PURPOSE + + + c This subroutine calculates lai, plant height, biomass, rootmass c distribution, rooting depth, and yield. c + + + KEYWORDS + + + c lai, biomass c + + + ARGUMENT DECLARATIONS + + + integer bnslay, bc0idc real bszlyd(*),bczrtd, bcrhi, bcehu0, bczmxc, bczmrt, bcmrtz(*), c r slai, dmag, r bctmin, bctopt, bc0ck, r bc0be1,bc0alf,bc0blf,bc0clf,bc0dlf, r bc0arp,bc0brp,bc0crp,bc0drp,bc0aht,bc0bht,bc0ssa, * bc0ssb,bc0sla,bc0hue, * bwtdav, bwtdmn, bweirr, bhfwsf,slait logical bm0cif c + + + ARGUMENT DEFINITIONS + + + c bc0alf - leaf partitioning parameter c bc0arp - rprd partitioning parameter c bc0aht - height s-curve parameter c ano3 - amount of applied N (t/ha) c bc0blf - leaf partitioning parameter c bc0brp - rprd partitioning parameter c bc0bht - height s-curve parameter c bc0clf - leaf partitioning parameter c bc0crp - rprd partitioning parameter c bc0ck - extinction coeffficient (fraction) c clfe - leaf area at emergence (cm^2/plant) c bc0dlf - leaf partitioning parameter c bc0drp - rprd partitioning parameter c dmag - stress adjusted cummulated aboveground biomass (t/ha) c bcrhi - yield index of a crop (ratio) c bczmxc - maximum potential plant height (m) c bc0hue - relative heat unit for emergence (fraction) c bc0idc - crop type:annual,perennial,etc c bnslay - number of soil layers c bczrtd - root depth (m) c bcmrtz - root biomass (by depth) c bczmrt - maximum root depth c bc0sla - specific leaf area (cm^2/g) c bc0ssa - parameter of the specific stem area function (cm^2/g) c bc0ssb - parameter of the specific stem area function (cm^2/g) c bctmin - base temperature (deg. C) c bctopt - optimum temperature (deg. C) c bc0be1 - CO2 concentration(ppm) c bszlyd - depth from top of soil to botom of layer, m c bm0cif - flag signalling start of crop growth c + + + COMMON BLOCKS + + + include 'w1clig.inc' include 'p1werm.inc' C *** include 'm1sim.inc' include 'm1subr.inc' include 'm1flag.inc' include 'm1dbug.inc' include 'c1db1.inc' include 'c1glob.inc' c local includes include 'crop/cgrow.inc' include 'crop/cenvr.inc' include 'crop/cparm.inc' include 'crop/csoil.inc' include 'crop/chumus.inc' include 'crop/cfert.inc' include 'crop/p1crop.inc' c + + + LOCAL VARIABLES + + + real slaiy,par,apar,x2, r drw,pddm,pdrw,pddmag,p_lf,p_rp, r pdstwt,pdlfwt,pdrpwt,dryf,pdstarea,pdlfarea,pslaix, r pdht,dht,abr,ws,xxx1,p_st,xsla,dlfwt,dstwt,drpwt, r ddmag,dlfarea,dstarea,hux,ff,pclfwty,ffa,ffw,ax1,ax2, r ffwp,ffap,hui0f, pdrd,xdmag,xclfwt,xcstwt, r xcrpwt,xyld,xrw, ydmag,xw,gif, r pcstwty,yldy integer day, mo, yr, doy integer i ! added declarations after uncommenting some code - LEW real dd ! added declarations after uncommenting some code - LEW real wcg ! added declarations after uncommenting some code - LEW real wft, wfl(mnsz) ! missing variable declarations - 5/15/99 - LEW ! weight fraction total and weight fraction by layer ! used to distribute root mass into the soil layers real xz c + + + LOCAL VARIABLE DEFINITIONS + + + c abr - plant abrasion factor c apar - intercepted photosynthetically active radiation (MJ/m2) c ax1 - leaf weight c ax2 - leaf weight c ceta - cummulative c cht1 - interim plant height variable c chty - plant height on day i-1 c clfwty - leaf dry weight on day i-1 (t/ha) c cta - actual transpiration (mm) c day - day of the month c daye - days after emergence c ddmag - incement in biomass (t/ha) c dht - daily height (m) c dlai - incremnet in lai c dlfarea - increment in leaf area (cm^2/palnt) c dlfwt - increment in leaf dry weight (t/ha) c dmag - stres adjusted above ground biomass (t/ha) c doy - day of year c drd - increment in root depth (m) c drpwt - increment in reproductive mass (t/ha) c drw - increment in root weight (t/ha) c dstarea - increment in stem area (cm^2/plant) c dstwt - increment in dry weight of stem (t/ha) c ea - actual vapor pressure (kPa) c es - saturated vapor pressure (kPa) c ff - senescence factor (ratio) c ffa - leaf senescnce factor (ratio) c ffap - leaf senescence factor under non-stress conditions (ratio) c ffw - leaf weight reduction factor (ratio) c ffwp - lea weight reduction factor for non-stressed conditions (ratio) c frst - frost damage factor c gif - grain index c bcehu0 - relative gdd at start of senescence c hui0f - relative gdd c hux - relative gdd c mo - month c p_lf - leaf partitioning ratio c p_rp - reproductive partitioning ratio c p_st - stem paratitioning ratio c par - photosynthetically active radiation (MJ/m2) c pclfwtx - potential leaf weight (t/ha) c pclfwty - potential leaf weight on day i-1 (t/ha) c pcstwty - potential stem area (day i-1) c pddm - increment in potential dry matter (t/ha) c pddmag - increment in potential aboveground dry matter (t/ha) c pdht - increment in potential height (m) c pdlfarea - increment in potential leaf area (cm^2/plant) c pdlfwt - increment in potential leaf dry weight (t/ha) c pdrd - increment in root length (m) c pdrpwt - increment in reproductive dry matter (t/ha) c pdrw - increment in root dry matter (t/ha) c pdstarea - increment in stem area (cm^2/plant) c pdstwt - increment in stem dry weight (t/ha) c pslaix - maximum potential lai (cm^2/cm^2) c rwl - root weight on day i-1 (t/ha) c slai - green leaf area index c slai0 - slai at start of senescence c slaiy - lai on day i-1 c ws - water stress factor (ratio) c wsf - water stress factor (ratio) c x1 - vapor pressure deficit factor c x2 - CO2 factor c xclfwt - leaf dry weight (g/plant) c xcrpwt - reproductive dry weight (g/plant) c xcstwt - stem dry weight (g/plant) c xdmag - increment in abovegound dry weight (g/plant) c xrw - stem dry weight (g/plant) c xsla - specific leaf area (cm^2/g) c xw - absolute value of minimum temprature c xx - lai adjustment parameter c xxx1 - stress factor (ratio) c xyld - economic yield (g/plant) c ydmag - aboveground biomass (lb/ac) c yr - year c c c + + + FUNCTIONS CALLED + + + c integer dayear c + + + SUBROUTINES CALLED + + + c huc1 c nup c npup c najn c najna c nuts c temps c + + + OUTPUT FORMATS + + + 2130 format(1x,i2,1x,i2,1x,2i4,1x,i4,12(f9.2),1x,f6.3,1x,2i4) c2130 format(1x,i2,'/',i2,'/',2i4,1x,i4,12(f9.2),1x,f6.3,1x,i4) c & 1x,3(f4.2,1x)) 2014 format(i4,1x,i6,1x,5(f6.0,1x),1x,5(f5.2,1x),1x,f5.2,1x,i4) c + + + END OF SPECIFICATIONS + + + call caldatw(day, mo, yr) doy = dayear (day, mo, yr) C c stop growth calculations if a killing harvest occurred or crop has matured c and before emergence or during dormancy C23456789*123456789*123456789*123456789*123456789*123456789*123456789*12 if (am0hrvfl.eq.2.or.hui.gt.1.000.or.hui.lt.bc0hue) then if (am0hrvfl.eq.2) then C *** am0cgf = .false. bczrtd=0. endif c output to the 'junk.out' file for debugging purposes : 5/21/98 c write(61,151)doy, shu, hrlt, am0drmfl,am0hrvfl c151 format (i4,1x,2(f6.1,1x),i2,2i5) write(61,151)doy, shu, phu,hui,bc0hue,am0hrvfl 151 format (i4,1x,4(f8.2,1x),i2) return endif daye=daye+1 c reduce green leaf area and leaf mass during winter dormancy period if (am0drmfl.eq.1) then c initialize heat units c shu=0.25*phu if (bwtdmn.lt.-2.0) then xw=abs(bwtdmn) frst=sqrt((1.-xw/(xw+exp(a_fr-b_fr*xw)))+0.000001) if (frst.gt.1.) frst=1. c reduce green leaf area due to frost damage (10/1/99) clfarea=clfarea*frst slai=clfarea/parea if (slai.lt.0.1) then slai=0.1 clfarea=slai*parea endif c clfwt=clfwt*frst**0.25 c cstwt=cstwt*frst**0.125 c cht=cht*frst**0.125 c rw=rw*frst**0.125 c dmag=clfwt+cstwt c if (dmag.lt.0.03) dmag=0.03 endif goto 887 endif C c if (am0drmfl.eq.1) goto 887 c Start biomass calculations c bweirr is total shortwave radiation and a factor of .5 is assumed C to get to the photosynthetically active radiation par=0.5*bweirr ! C-4 c calculate intercepted PAR, which is the good stuff less what hits the ground apar=par*(1.-exp(-bc0ck*slai)) ! C-4 c calculate biomass conversion efficiency (be) at the given co2 & vpd level c es=0.61078*exp(17.269*bwtdav/(237.3+bwtdav)) c ea=es*rh c vpd=es-ea c x1=amax1(vpd-1.,-.5) c co2=bc0be1 c x2=100.*co2/(co2+exp(a_co-b_co*co2)) C calculate the radiation use efficiency x2=100.*bc0be1/(bc0be1+exp(a_co-b_co*bc0be1)) c be=x2-wavp*x1 c following 1 line added on 4/29/94 c if (be.gt.x2) be=x2 be=x2 ! C-4 c calculate potential biomass conversion (t/ha/day) pddm=0.001*be*apar ! C-4 C partition the increase between the roots and above ground pdrw=pddm*(.4-.2*hui) ! C-5 pddmag=pddm-pdrw pdm=pdm+pddm prw=prw+pdrw pdmag=pdmag+pddmag c new partitioning algorithm added (5/27/94) c calculate leaf, stem, and rprd partitioning parameters. p_lf=bc0alf+bc0blf/(1.+exp(-(hui-bc0clf)/bc0dlf)) if (p_lf.lt.0.) p_lf=0. p_rp=0. c skip reproductive growth calculations for forage crops c if (bc0idc.eq.3) goto 173 c if (hui.gt.0.45) p_rp=bc0arp+bc0brp/(1.+exp(-(hui-bc0crp)/bc0drp)) p_rp=bc0arp+bc0brp/(1.+exp(-(hui-bc0crp)/bc0drp)) if (p_rp.lt.0.) p_rp=0. if (p_rp.gt.1.) p_rp=1. c173 continue p_st=1.-p_lf-p_rp if (p_st.lt.0.) then write(*,*) ' growth: stem partition less than 0' p_lf = p_lf / (p_lf + p_rp) p_rp = p_rp / (p_lf + p_rp) p_st = 0.0 endif c if (bc0idc.eq.3) goto 183 if (hui.ge.bcehu0) then p_lf=0. p_st=0. p_rp=1. endif c183 continue c calculate potential leaf, stem and reproductive masses (t/ha) pdlfwt=p_lf*pddmag pdrpwt=p_rp*pddmag pdstwt=pddmag-pdlfwt-pdrpwt if (pdstwt.lt.0.) then write(*,*) ' growth: stem partition weight less than 0' pdstwt = 0.0 endif pclfwt=pclfwt+pdlfwt pcstwty=pcstwt pcstwt=pcstwt+pdstwt pcrpwt=pcrpwt+pdrpwt c no-linear adjustment of specific leaf area c asla=a_adj*hui**b_adj c if (asla.gt.2.) asla=2. c linear ajustment of specific leaf area during the early growth stages c asla=1. c if (hui.lt..4) then c asla=a_adj+hui*b_adj c endif c xsla=bc0sla*(1.72-1.15*hui) c if (dryf.le.0.) dryf=1. c xsla=bc0sla*dryf c calculate potential leaf and stem areas (cm^2/plant) and area indexes. c 0.01 is a conversion factor from t/ha to g/plant. pdlfarea=pdlfwt*.01*bc0sla*parea c pdstarea=pdstwt*.01*bc0ssa*parea c These two lines replaced by the next three C pdstarea=(pdstwt*.01*parea)*bc0ssa*bc0ssb* C 1 (pcstwt*.01*parea)**(bc0ssb-1.) pdstarea = bc0ssa*(pcstwt*0.01*parea)**bc0ssb & -bc0ssa*(pcstwty*0.01*parea)**bc0ssb if (pdstarea .lt. 0.0) pdstarea = 0.0 pcstarea=pcstarea+pdstarea pcstarea=pcstarea+pdstarea pclfarea=pclfarea+pdlfarea c pcstarea=pcstarea+pdstarea c pdslai=pdlfarea/parea c pdssai=pdstarea/parea pslai=pclfarea/parea if (hui.le.bcehu0) then pslaix=pslai pclfwtx=pclfwt endif pssai=pcstarea/parea c added method (different from EPIC) of calculating plant height c pht=cummulated potential height,pdht=daily potential height c cht=cummulated actual height,adht=daily actual height, bc0aht,bc0bht are c height-scurve parameters (formerly lai parameters) huf=.01+1./(1.+exp((hui-bc0clf)/bc0dlf)) pchty=pcht pcht=bczmxc*huf if (pcht.gt.bczmxc)pcht=bczmxc pdht=pcht-pchty c ------------------------------------------------------------------- c calculate N & P demand and supply c call nuse c calculate N & P uptake with increase in supply if necessary c call najn c call najna c claculate N stress c call nuts (un1,un2,sn) c call nuts (sun,un2,sn) c calculate P stress c call nuts (up1,up2,sp) c calculate temperature stress call temps (bctmin,bctopt, bwtdmn, bwtdav) c write (11,4011) jd,sn,sp,ts c4011 format(i3,1x,3(f6.4,1x)) c calculate water stress c call waters c call abrasion abr=1. sn=1. sp=1. ws=bhfwsf c changed from above statement to the one below : 2/27/96 c ws=(bhfwsf+0.00001)**.75 c strs=min(sn,sp,ts,ws,abr) xxx1=min(ts,ws) strs=xxx1 c the line below was added to test model without stress:8/26/99: A.Retta c strs=1.0 c if (hui.lt.0.25) strs=xxx1**2 c if (hui.gt.huilx) strs=sqrt(xxx1) c strs=sqrt(strs) c calculate stress adjusted height (m) dht=pdht*strs cht=cht+dht c calculate stress adjusted masses (t/ha) dlfwt=pdlfwt*strs dstwt=pdstwt*strs drpwt=pdrpwt*strs clfwt=clfwt+dlfwt cstwt=cstwt+dstwt crpwt=crpwt+drpwt dmag=clfwt+cstwt+crpwt drw=pdrw*strs ddm=dlfwt+dstwt+drpwt+drw ddmag=dlfwt+dstwt+drpwt dm=dm+ddm rw=rw+drw c dmag=dmag+ddmag c c calculate stress adjusted leaf and stem areas, and leaf and stem area c indexes. dlfarea=pdlfarea*strs c if (je.eq.2) dlfarea=pdlfarea*strs*strs c if (je.eq.3) dlfarea=pdlfarea*strs*strs dstarea=pdstarea*strs clfarea=clfarea+dlfarea cstarea=cstarea+dstarea slai=clfarea/parea ssai=cstarea/parea if (hui.le.bcehu0) then slaix=slai clfwtx=clfwt ssaix=ssai endif c skip leaf senescence and leaf mass loss calculations for forage crops c if (bc0idc.eq.3) goto 174 c lag start of senescence hui0f=bcehu0-bcehu0*.1 if (hui.ge.hui0f) then hux=hui-bcehu0 ff=1./(1.+exp(-(hux-bc0clf/2.)/bc0dlf)) ffa=ff**0.25 ffw=ff**0.25 ffap=ff**0.25 ffwp=ff**0.125 c if (je.eq.1) ffw=ffa c slow leaf area senescence for corn and sorghum c if (je.eq.2.or.je.eq.3) ffa=ff**.125 slaiy=slai clfwty=clfwt pclfwty=pclfwt slai=slai*ffa c pslai=pslaix*ffap clfwt=clfwtx*ffw c pclfwt=pclfwtx*ffwp clfarea=slai*parea c pclfarea=pslai*parea ax1=pclfwty-pclfwt ax2=clfwty-clfwt pdmag=pdmag-ax1 dmag=dmag-ax2 pdm=pdm-ax1 dm=dm-ax2 c root senescence : 02/16/2000 (A. Retta) c write(*,*)'root before=',rw rw=rw*0.98 c write(*,*)'root after=',rw endif c174 continue c ***************** end of new algorithm ******************** su=8.0 c calculate rooting depth (eq. 2.203) and check that it is not deeper c than the maximum potential depth, and the depth of the root zone. c drd=2.5*bczmrt*dhuf c bczrtd=bczmrt*huf+0.1 c if (bczrtd.gt.bczmrt)bczrtd=bczmrt c if (bczrtd.gt.bszlyd(bnslay))bczrtd=bszlyd(bnslay) c aczrtd(am0csr)=bczrtd prdy=prd prd=bczmrt*huf+0.1 if (prd.gt.bczmrt)prd=bczmrt if (prd.gt.bszlyd(bnslay))prd=bszlyd(bnslay) pdrd=prd-prdy if (pdrd.lt.0.)pdrd=0. bczrtd=bczrtd+pdrd if (bczrtd.gt.bczmrt) bczrtd=bczmrt C This was commented out, but we turned it back on again (after correcting C for units - depth to bottom of layers is really in "mm" and this routine C thought it was in "m"). Also fixed bug in cinit() due to this units error C (za is now in "m" as expected). C 5/14/99 - LEW c determine botttom layer # where there are roots do 51 i=1,bnslay dd=bczrtd-(bszlyd(i)/1000.0) if (bczrtd.gt.(bszlyd(i)/1000.0)) goto 51 ir=i go to 52 51 continue 52 continue c distribute root weight into each layer (eq. 2.202--epic1910) wft=0.0 do 60 i=1,ir c the root distribution functions were taken from agron. monog. 31, equ. 26 c on page 99. wcg=2.0 for corn and soybeans. wcg=2.0 wfl(i)=(1.0-za(i)/3.0)**wcg wft=wft+wfl(i) 60 continue c distribute root weight into each layer !distribute initial root mass among layers !triggering on top layer having zero live root mass to begin with if (bcmrtz(1) .le. 0.0 ) then do 61 i=1,ir ! corrected units. bcmrtz is in "kg/m^2", rw is in "t/ha" bcmrtz(i)=((rw-drw)*wfl(i)/wft)*0.1 61 continue endif do 62 i=1,ir ! corrected units. bcmrtz is in "kg/m^2", drw is in "t/ha" bcmrtz(i)=bcmrtz(i)+(drw*wfl(i)/wft)*0.1 c write (11,3011) jd,bcmrtz(i),wfl(i),wft,drw c 3011 format(' jd=',i3,' bcmrtz=',f7.5,' wfl=',f7.5,' wft=',f7.5,' drw=', c 1f7.5) 62 continue c calculate yield (t/ha).Note that bcrhi is the harvest index c yld=bcrhi*crpwt c the above 1 statement replaced by the following 9/24/96 gif=-.02+1./(1+exp(-(hui-0.64)/.075)) if (gif.lt.0.) gif=0. if (gif.gt.1.) gif=1. c yld=gif*bcrhi*crpwt c chaff=crpwt-yld c if (chaff.lt.0.)chaff=0. yldy=yld yld=dmag*bcrhi*gif c following 3 lines added on 9/17/99 : Retta if (yld.lt.yldy) yld=yldy if (yld.gt.0.8*crpwt) yld=0.8*crpwt c if (yld.lt.0.9*crpwt) yld=0.9*crpwt chaff=crpwt-yld if (chaff.lt.0.) chaff=0. c convert masses from t/ha to g/plant xdmag=dmag*parea*.01 xclfwt=clfwt*parea*.01 xcstwt=cstwt*parea*.01 xcrpwt=crpwt*parea*.01 xyld=yld*parea*.01 xrw=rw*parea*.01 887 continue c the following write statements are for 'crop.out' c am0cfl is flag to print crop submodel output if (am0cfl .eq. 1) then c write(*,*) day,mo,yr c write(*,131) day,mo,yr c131 format(' growth:day,mo,yr=',3i6) write(17, 2130) day, mo, yr, am0csr, ac0id(am0csr), slai,ssai, & dmag,clfwt,cstwt,crpwt,yld,rw,bczrtd,cht,cta,prcp,hui,doy,daye c write(17, 2130) day, mo, yr, am0csr, ac0id(am0csr), slai,ssai, c & xdmag,xclfwt,xcstwt,xcrpwt,xyld,xrw,bczrtd,cht,cta,ceta,hui,doy end if c write(*, 2130) day, mo, yr, am0csr, ac0id(am0csr), slai, ssai, c & dmag,clfwt,cstwt,crpwt,yld,rw,bczrtd,cht,hui,doy,strs,ws,ts c end if c write (39,2051)jd,cnt,uno3,sunn,un1,un2,cpt,up1,up2,sup,ir c c 2051 format (1x,i3,9f8.4,1x,i3) c write (10,2131)jd c write (10,2133)(ap(i),bcmrtz(i),wno3(i),un(i),up(i),i=1,7) c 2133 format (5(1x,f9.4)) c 2131 format (' jd=',i3) c dmagy = dmag c transfer local variables to globals C *** debugging hack if (slai .lt. 0.04) slai = 0.04 C *** end of hack C comments by LEW - Mon Apr 26 04:59:58 CDT 1999 C I believe that we are setting the daily crop variables here C after a day of growth. The 0.1 constant is probably changing C units from t/ha (EPIC growth model units) back to kg/m^2 (WEPS units) acrlai(am0csr) = slai acmyld(am0csr) = yld * 0.1 c subtract "yield" mass from above ground biomass acmst(am0csr) = (dmag -yld) * 0.1 c acmst(am0csr) = (clfwt+cstwt) * 0.1 acmrt(am0csr) = rw * 0.1 aczht(am0csr) = cht !aczrtd(am0csr) = bczrtd !not necessary - LEW c acmyld(am0csr)=xyld c acmst(am0csr)=xdmag c acmrt(am0csr)=xrw C These appear to be missing - LEW acrsai(am0csr) = ssai acm(am0csr) = acmst(am0csr)+acmrt(am0csr)+acmyld(am0csr) c acmrtz(layer,am0csr) = ? ! now done earlier within the code directly - LEW c acrsaz(layer,am0csr) = ? ! will be included later when needed by erosion - LEW c acrlaz(layer,am0csr) = ? ! will be included later when needed by erosion - LEW c acdstm(am0csr) = ? acffcv(am0csr) = 0.0 c crop canopy cover estimate from light extinction coefficient c if (hui.le.0.8)slait=slaix c if (hui.gt.0.8)slait=slait*0.99 slait=slai if (slait.lt.slaix) slait=slaix c write(*,*) 'growth: ',slait,bc0ck acfscv(am0csr) = 1.0-exp(-bc0ck*slait) if (hui.gt.0.8)acfscv(am0csr)=acfscv(am0csr)*(clfwt/clfwtx)**0.2 c total crop cover acftcv(am0csr) = acfscv(am0csr) + acffcv(am0csr) c convert biomass to lb/acre and output to file part.out ydmag=dmag*893. if (am0cdb.eq.2) then write(57,2014)doy,yr,clfwt*893.,cstwt*893,crpwt*893.,yld*893., & ydmag,p_lf,p_st,p_rp,ws,ts,slai,daye endif C Why do we do this "after" we have updated the WEPS global variable "aczrtd"? C LEW if (am0hrvfl.eq.2) then C *** am0cgf = .false. bczrtd=0. endif c output to the 'junk.out' file for debugging purposes : 5/21/98 c write(61,151)doy, shu, phu, am0drmfl,am0hrvfl,iprint c endif c endif return end