c file: cinit.for C NOTE! depth to bottom of soil layers coming in are in "mm" and not "m" as C this routine thinks. I believe I have "quick fixed" the problems here. C 5/14/99 - LEW subroutine cinit(bnslay, bszlyt, bszlyd, bsdblk, bsfcce, bsfcec, * bsfsmb, bsfom, bsfcla, bs0ph, & bc0bn1, bc0bn2, bc0bn3, bc0bp1, bc0bp2, * bc0bp3, bsmno3, & bc0fd1, bc0fd2,bctmin, & cc0fd1, cc0fd2, bc0be1, bc0be2, cc0be1, * cc0be2, & bc0sla,bc0lfe,bc0idc) c Author : Amare Retta c + + + PURPOSE + + + c This subroutine initializes parameters c + + + KEYWORDS + + + c Initialization c + + + ARGUMENT DECLARATIONS + + + integer bnslay, bc0idc real bszlyt(*), ! added so a local variable would be set correctly - LEW * bszlyd(*), bsdblk(*), bsfcce(*), bsfcec(*), bsfsmb(*), * bsfom(*), bsfcla(*), bs0ph(*), r bc0fd1, bc0fd2, bctmin, r cc0fd1, cc0fd2, bc0be1, bc0be2, cc0be1, cc0be2, r bc0bn1, bc0bn2, bc0bn3, bn4, bc0bp1, bc0bp2, bc0bp3, * bp4, bsmno3, r bc0sla, bc0lfe c + + + ARGUMENT DEFINITIONS + + + c bsmno3 - fertilizer applied as NO3 c bsdblk - bulk density of a layer (g/cm^3=t/m^3) c bsfcce - calcium carbonate (%) c bsfcla - % clay c bsfom - percent organic matter c bsfcec - cation exchange capacity (cmol/kg) c dmag - above ground biomass (t/ha) c bnslay - number of soil layers c bs0ph - soil pH c bsfsmb - sum of bases (cmol/kg) c bszlyd - depth from top of soil to botom of layer, m c + + + GLOBAL COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'm1flag.inc' include 'm1dbug.inc' include 'm1sim.inc' include 'c1info.inc' include 'c1gen.inc' include 'm1subr.inc' include 'w1clig.inc' c + + + LOCAL COMMON BLOCKS + + + include 'crop/csoil.inc' include 'crop/chumus.inc' include 'crop/cfert.inc' include 'crop/cgrow.inc' include 'crop/cenvr.inc' include 'crop/cparm.inc' include 'crop/p1crop.inc' *$reference c + + + LOCAL VARIABLES + + c character*20 cpnm integer i,n, pdate,hdate,j real bsa, dg, dg1, hlmn, wt1, x1, xz, jreal real xx, sdmn, sdmx, hlmx,x(13),y(13),y2(13),phux,sphu,yp1, & ypn real tat, tad c + + + LOCAL VARIABLE DEFINITIONS + + + c alt - idex of crop toleramce to aluminum saturation (1-5;1=sensitive, c 5=tolerant) c bc0be1 - co2 concentration(ppm) (1st x-coordinate in the biomass c conversion efficiency s-curve) c bc0be2 - co2 concentration(ppm) (2nd x-coordinate in the biomass c conversion efficiency s-curve) c bc0fd1 - minimum temperature below zero (c) (1st x-coordinate in the c frost damage s-curve) c bc0fd2 - minimum temperature below zero (c) (2nd x-coordinate in the c frost damage s-curve) c bc0idc - crop type:annual,perennial,etc c bctmin - base temperature (deg. c) c bn - N uptake parametres at different stages of growth(ratio) c bp - P uptake parametres at different stages of growth(ratio) c bsa - base saturation (%?) c cc0be1 - biomass conversion efficiency (kg/ha/mj) (1st y-coordinate c in the biomass efficiency s-curve) c cc0be2 - biomass conversion efficiency (kg/ha/mj) (2nd y-coordinate c in the biomass efficiency s-curve) c cc0fd1 - fraction of biomass lost each day due to frost c (1st y-coordinate in the frost damge s-curve) c cc0fd2 - fraction of biomass lost each day due to frost c (2nd y-coordinate in the frost damge s-curve) c ch - interim variable in calculations of day length c ck - extinction coeffficient (fraction) c co2 - co2 concentration in the atmosphere (ppm) c cpnm - crop name c dg - soil layer thickness (mm) c dg1 - previous value of dg c dlax1 - fraction of grwing season (1st x-coordinate in lai s-curve) c dlay1 - fraction of maximum lai (1st y-coordinate in the lai s-curve) c dlax2 - fraction of grwing season (2nd x-coordinate in lai s-curve) c dlay2 - fraction of maximum lai (2nd y-coordinate in the lai s-curve) c dtm - days to maturity c frs1 - same as bc0fd1 (value needed in cfrost.for:11/20/96) c frs2 - same as bc0fd2 (value needed in cfrost.for:11/20/96) c h - interim variable in calculations of day length c hdate - day of year harvest occurs c hi - harvest index of a crop c hix1 - ratio of actual to potential et (1st x-coordinate in the c water stress - harvest index s-curve) c hiy1 - fraction of reduction in harvest index (1st x-coordinate in c the water stress - harvest index s-curve) c hix2 - ratio of actual to potential et (2nd x-coordinate in the c water stress - harvest index s-curve) c hiy2 - fraction of reduction in harvest index (2nd x-coordinate c in the water stress - harvest index s-curve) c hlmn - minimum daylength for a site (hr) c hlmx - maximum daylength for a site (hr) c hlmn0 - minimum daylength for a site plus one (hr) c hmx - maximum potential plant height (m) c hui0 - heat unit index when leaf senescence starts. c irint - flag for printing end-of-season values : 10/6/99 c nc - crop number c pdate - day of year planting oddurs c phu - potential heat units for crop maturity (deg. c) c phux - heat unit accumulated before before planting c rbmd - biomass-energy ratio decline factor c rdmx - maximum root depth c rlad - crop parameter that governs leaf area index decline rate c s11x1 - soil labile p concentration (ppm) (1st x-coordinate in the c p uptake reduction s-curve) c s11y1 - p uptake restriction factor (1st y-coordinate in the p c uptake reduction s-curve) c s11x2 - soil labile p concentration (ppm) (2nd x-coordinate in the c p uptake reduction s-curve) c s11y2 - p uptake restriction factor (2nd y-coordinate in the p c uptake reduction s-curve) c s8x1 - scaled ratio of actaul to potential n or p (1st x-coordinate c in the n stress factor s-curve) c s8y1 - n or p stress factor (1st ycoordinate in the n or p stress c s-curve) c s8x2 - scaled ratio of actaul to potential n or p (2nd x-coordinate c in the n stress factor s-curve) c s8y2 - n or p stress factor (2nd ycoordinate in the n or p stress c s-curve) c slaix - maximum potential lai c sphu - c tad - c tat - c topt - optimum temperature (deg. c) c wavp - parameter relating vapor pressure deficit to be c wsyf - minimum crop harvest index under drought c wt1 - convert N, P, etc. conc. from g/t to kg/t c x1 - temporary bulk density value c xx - depth of above layer c xz - c ytn - tangent of the sun's declination angle c + + + SUBROUTINES CALLED + + + c scrv1 c sdst c nconc c + + + INPUT FORMATS + + + c + + + OUTPUT FORMATS + + + c2020 format ('day/mo/year subr je slai pdm dm rw c & rd hia dmg yld shu hui fleaf leafwt stem c &wt rprwt cht cht2') c2016 format(/,' hum ap wno3 wp rsd fop c : rtn pnm bk psp fon wn wmn c : op') c2017 format (1x,14f9.3) c2018 format (/,1x,'tap=',f8.2,' tno3=',f8.2,' tp=',f8.2,' twmn=',f8.2, c 1 ' top=',f8.2,' tmp=',f8.2,' trsd=',f8.2) c 2110 format (5x,' a_co2=',f6.3,' b_co2=',f6.3,' a_frst=',f6.3, &' b_frst=',f6.3) c 1' b_s11=',f6.3,' a_s8=',f6.3,' b_s8=',f6.3,/) c2115 format(1x,'bc0bn1=',f7.6,' bc0bn2=',f9.6,' bc0bn3=',f7.6,' bn4=',f9.6, c 1 ' bc0bp1=',f7.6,' bc0bp2=',f7.6,' bc0bp3=',f7.6,' bp4=',f9.6,/) c c2990 format (' bszlyd bsdblk bsfcla wp bs0ph bsfsmb bsfom c 1 bsfcce bsfcec wno3 ap rsd psp wn') c c2992 format (13f8.3,1x,f8.1) c + + + END OF SPECIFICATIONS + + + c not sure if this is where this should go but we need to initialize c this and other variables when planting a new crop - jt 6/3/94 c parea=2000. c reg=1. c write (21,4002) c 4002 format(' jd hui slai ssai clfwt cstwt crpwt c 1 dmag rw dm cht clfarea cstarea strs c 2 xsla') c write (34,4003) c 4003 format(' jd hui pslai pssai pclfwt pcstwt pcrpwt c 1 pdmag prw pdm pcht pclfarea pcstarea strs c 2 xsla') c write (35,4004) c 4004 format(' jd hui pdlfwt dlfwt pdstwt dstwt pdrpwt c 1 drpwt pdrw drw pddmag ddmag p_lf p_st c 2 p_rp') c write (36,4005) c 4005 format(' jd tmn tmx rad rain humd hlmn0 c c hrlt') c write(21,1311) rdmx, bctmin, topt, wavp c write(21,1312) bc0fd1, bc0fd2,cc0fd1, cc0fd2 c write(21,1313) bc0be1, bc0be2, cc0be1, cc0be2 c write(21,1314) ssaa,ssab,bc0sla,bc0lfe c c1311 format(' rdmx=',f9.4,' bctmin=',f9.4,' topt=',f9.4,' wavp=',f9.4) c1312 format('bc0fd1=',f9.4,' bc0fd2=',f9.4,' cc0fd1=',f9.4,' cc0fd2='f9.4) c1313 format('bc0be1=',f8.3,' bc0be2=',f8.3,' cc0be1=',f8.3,' cc0be2=',f8.3) c1314 format('ssaa=',f8.3,' ssab=',f8.3,' bc0sla=',f8.3,' bc0lfe=',f8.3) c open debug file c open (unit=56, file='initdb.out',status='unknown') c Initialize c initialize variables needed for season heat unit computation: added on c 3/16/1998 (A. Retta) data x / 15,45,74,105,135,166,196,227,258,288,319,349,380 / c transfer average monthly temperatures from the global array to a local array do i=1,12 y(i)=awtmav(i) end do y(13)=y(1) c added algorithm to compute yield adjustment factor; 10/18/99 c tdmag=15.0 c if (dmag.le.0.)yaf=1.0 c if (dmag.gt.0.)yaf=log(tdmag/0.00001)/log(dmag/0.00001) c write (*,*)'yaf, dmag, tdmag ' ,yaf,dmag,tdmag c initialize variables phu=0. dmag=0. clfwt=0. cstwt=0. pdate=0 write(60,*) 'plant population (plants/m2)=', acdpop(am0csr) c convert to m2/plant parea=10000./acdpop(am0csr) clfarea=bc0lfe cstarea=.1*clfarea pclfarea=clfarea pcstarea=cstarea slai=clfarea/parea ssai=cstarea/parea clfwt=(clfarea/bc0sla)*100./parea cstwt=0.2*clfwt dmag=clfwt+cstwt rw=dmag prw=rw pdmag=dmag pdm=pdmag+prw dm=pdm pclfwt=clfwt pcstwt=cstwt shu=0. hrlt=0. hrlty=0. hufy=0. cht=0.05 pcht=0.05 pchty=0.0 crpwt=0.0 prd=0. prdy=0. cta=0. ceta=0. prcp=0. c am0hrvfl=0 !9/28/99 ctp=0. chaff=0. daye=0 frs1=bc0fd1 frs2=bc0fd2 iprint = 0 !10/6/99 c Calculate minimum and maximum daylength for a location xlat = amalat c sd=0.41021*sin(2.*3.1416*(354.-80.25)/365.) c hlmn=7.64*acos(-1.*tan(2.*3.1416*amalat/365.)*tan(sd)) c if (amalat.lt.0.) then c sd=0.41021*sin(2.*3.1416*(171.-80.25)/365.) c hlmn=7.64*acos(-1.*tan(2.*3.1416*amalat/365.)*tan(sd)) c end if c hlmn0=hlmn+1.0 c above 7 lines replaced by folowing 8 lines sdmn=0.41021*sin(2.*3.1416*(354.-80.25)/365.) sdmx=0.41021*sin(2.*3.1416*(173.-80.25)/365.) if (amalat.lt.0.) then sdmn=0.41021*sin(2.*3.1416*(173.-80.25)/365.) sdmx=0.41021*sin(2.*3.1416*(354.-80.25)/365.) end if hlmn=7.64*acos(-1.*tan(2.*3.1416*amalat/360.)*tan(sdmn)) hlmx=7.64*acos(-1.*tan(2.*3.1416*amalat/360.)*tan(sdmx)) hlmn0=hlmn+(hlmx-hlmn)*0.18 c start calculation of seasonal heat unit requirement c first calculate heat units needed to accumulate before planting can be done c The equation for computing phux was taken from the PHU prgram of EPIC. sphu=0. tat=0. n=13 yp1=0. ypn=0. phux=0. sphu=0. phux = 30.0 * 9.5 / (sqrt(amalat + 1.)) c call cubic spline interpolation routines call spline (x, y, n, yp1, ypn, y2) c ....... c do i=1,13 c write(56,*)'y2(i)=',y2(i) c end do c write (56,*) 'bctmin=',bctmin,' phux=',phux c ...... do j = 1, 365 jreal = j call splint(x, y, y2, n, jreal, tat) c......write(56,*)'n=',n,' j=',j,' tat=',tat tad = tat - bctmin if (tad .lt. 0.) tad = 0. if (sphu .lt. phux) then sphu = sphu + tad if (sphu .ge. phux) then pdate = j hdate = pdate + dtm endif endif c calculate heat units after planting if (pdate.gt.0) then if (j .gt. hdate) goto 27 phu = phu + tad endif end do 27 continue c increase final phu for winter wheat by 20% to account for fall growth. This was c adjustment is based on some test done by A. Retta (6/12/1998) if (bc0idc .eq. 5) then phu = phu + 0.5 * phu If (phu.gt.2000.)phu=2000.0 endif c ...... end of seasonal heat unit calculator algorithm......... write(60,*) 'computed plant date=',pdate, & ' computed harvest date=', hdate write(60,*)'computed total heat units (phu)=',phu c Calculate s-curve parameters c write(*,*)bc0be1,cc0be1 cc0be1 = cc0be1 * 0.01 cc0be2 = cc0be2 * 0.01 call scrv1(bc0be1,cc0be1,bc0be2,cc0be2,a_co,b_co) call scrv1(bc0fd1,cc0fd1,bc0fd2,cc0fd2,a_fr,b_fr) call scrv1(s11x1,s11y1,s11x2,s11y2,a_s11,b_s11) call scrv1(s8x1,s8y1,s8x2,s8y2,a_s8,b_s8) if (am0cdb.eq.2) write (60,2110)a_co,b_co,a_fr,b_fr c calculate bc0bn1,bc0bn2,bc0bn3,bn4,bc0bp1,bc0bp2,bc0bp3,bp4 call nconc (bc0bn1,bc0bn2,bc0bn3,bn4) call nconc (bc0bp1,bc0bp2,bc0bp3,bp4) c write (33,2115) bc0bn1,bc0bn2,bc0bn3,bn4,bc0bp1,bc0bp2,bc0bp3,bp4 c This section estimates soil test data and calculates initial c values of no3, labile p, fresh organic matter, humus, etc in each layer. c write (33,2990) c write (33,2992) (bszlyd(i),bsdblk(i),bsfcla(i),wp(i),bs0ph(i),bsfsmb(i), c 1bsfom(i),bsfcce(i),bsfcec(i),wno3(i),ap(i),rsd(i),psp(i),wn(i),i=1,10) do 117 i=1,bnslay x1=bsdblk(i) c dg=soil layer thickness in mm; xx=bszlyd(i-1); wt(i)=wt of soil(t/ha-mm)(conve c rt to t/ha ??) wt1=convert N, P, etc. conc. from g/t to kg/t c bsfcla(i)=100.-san(i)-sil(i) ! NOTE: bszlyd is in "mm" and not "m" - 5/14/99 - LEW ! Commenting out unnecessary code here C*** c this if statment was added so that subscript will not = 0 when i = 1 C*** if (i .eq. 1) then C*** xx = 0. C*** else C*** xx = bszlyd(i-1) C*** end if ! Since "dg" is just the layer thickness in "mm", we have set it ! to the now included WEPS global layer thickness variable here. !dg=1000.*(bszlyd(i)-xx) dg = bszlyt(i) wt(i)=bsdblk(i)*dg*10. wt1=wt(i)/1000. c calculate depth to the middle of a layer ! NOTE: local crop global variable za is expected to be in "m" ! so we have "fixed" it here by dividing by 1000.0 in the ! appropriate places. - 5/14/99 - LEW !za(i)=za(i)/2. ! I think he wanted to set za(1) = bszlyd(1)/2.0 !if (i.gt.1) za(i)=bszlyd(i-1)+(bszlyd(i)-bszlyd(i-1))/2. if (i.eq.1) then za(i)=(bszlyd(i)/2.0)/1000.0 else za(i)= (bszlyd(i-1)+(bszlyd(i)-bszlyd(i-1))/2.0)/1000.0 endif c estimate initial values of: rsd(residue:t/ha);ap(labile P conc.:g/t);wno3 c (no3 conc.:g/t). dg1=previous value of dg. c call sdst (rsd,dg,dg1,i) call sdst (ap,dg,dg1,i) call sdst (wno3,dg,dg1,i) trsd=trsd+rsd(i) dg1=dg c calculate ratio (rtn) of active(wmn) to total(wn) N pools associated with c humus. yc=years of cultivation. c yc was set to 150. as suggested by A. Retta - jt 1/10/94 yc = 150. if (i.eq.1) rtn(1)=0.4*exp(-0.0277*yc)+0.1 if (i.gt.1) rtn(i)=rtn(1)*.4 c estimate bsfcec, bsa c there are some very questionable things done here - jt 1/10/94 c xx=bszlyd(i) if (bsfcec(i).gt.0.) then if (bsfcce(i).gt.0.) bsfsmb(i)=bsfcec(i) if (bsfsmb(i).gt.bsfcec(i)) bsfsmb(i)=bsfcec(i) bsa=bsfsmb(i)*100./(bsfcec(i)+1.e-20) bsfcec(i)=bs0ph(i) bsfsmb(i)=bsfcec(i) endif c calculate aluminum saturation of the soil c ALS=aluminum saturation of the soil layer (%) c RWT=root weight (t/ha) c als is not used any where else in the crop submodel, c commented out 1/10/94 - jt c als(i)=154.2-1.017*bsa-3.173*bsfom(i)-14.23*bs0ph(i) c if (als(i).lt.0.) als(i)=0. c if (als(i).gt.95.) als(i)=95. c if (bs0ph(i).gt.5.6) als(i)=0. c c root weithr is initialized in main/input.for c rwt(i)=0. c calculate amounts of N & P (kg/ha) from fresh organic matter(rsd), c assuming that N content of residue is 0.8%; fon & fop are in kg/ha. fon(i)=rsd(i)*8. fop(i)=rsd(i)*1.1 tfon=tfon+fon(i) tfop=tfop+fop(i) c initial organic(humus) N & P concentrations (g/t) in the soil if (wn(i).eq.0.) wn(i)=1000.*bsfom(i) if (wp(i).eq.0.) wp(i)=0.125*wn(i) c c Estimate psp(P sorption ratio), which is the fraction of fertilizer P that c remains in the labile form after incubation for different soil conditions. c The weathering status code ids is inputted. c ids=1:esitmate psp for calcareous soils without weathering information c ids=2:estimate psp for noncalcareous slightly weathered soils c ids=3:estimate psp for noncalcareous moderately weathered soils c ids=4:estimate psp for noncalcareous highly weathered soils c ids=5: input value of psp c finally estimate the flow coefficient bk between the active and stable c P pools (1/d) if (ids.eq.1) then psp(i)=0.5 if (bsfcce(i).gt.0.) psp(i)=0.58-0.0061*bsfcce(i) endif if (ids.eq.2) psp(i)=0.02+0.0104*ap(i) if (ids.eq.3) psp(i)=0.0054*bsa+0.116*bs0ph(i)-0.73 if (ids.eq.4) psp(i)=0.46-0.0916*alog(bsfcla(i)) if (psp(i).lt.0.05) psp(i)=0.05 if (psp(i).gt.0.75) psp(i)=0.75 bk(i)=0.0076 if (bsfcce(i).gt.0.) bk(i)=exp(-1.77*psp(i)-7.05) c calculate initial amount of active(pmn) and stable(op) mineral P pools c ap=initial amount of labile P(g/t);wt1=conversion factor to kg/ha pmn(i)=ap(i)*(1.-psp(i))/psp(i)*wt1 tmp=tmp+pmn(i) op(i)=4.*pmn(i) top=top+op(i) c calculate amount of active(readily mineralizable) humus N pool(wmn) and c total(active+stable) humus N pool(wn) wmn(i)=rtn(i)*wn(i) wn(i)=wn(i)-wmn(i) wn(i)=wn(i)*wt1 twn=twn+wn(i) wmn(i)=wmn(i)*wt1 twmn=twmn+wmn(i) c convert total(active+stable) humus P pool to kg/ha wp(i)=wp(i)*wt1 tp=tp+wp(i) c convert initial no3 & labile p in the soil to kg/ha ap(i) = ap(i) * wt1 tap = tap + ap(i) wno3(i) = wno3(i) * wt1 tno3 = tno3 + wno3(i) c calculate amount of humus in a layer (t/ha) c moved from original location xz = bsfom(i) * .0172 hum(i) = xz * wt(i) 117 continue c add applied fertilizer to the top layer wno3(1) = wno3(1) + bsmno3 tno3 = tno3 + bsmno3 c write (*,*) 'past cinit slai=',slai return end