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,dd,mm,yy,bcthud,bcthum)

c     Author : Amare Retta
c     + + + PURPOSE + + +
c     This subroutine initializes parameters

c     + + + KEYWORDS + + +
c     Initialization

c     + + + ARGUMENT DECLARATIONS + + +
      integer bnslay, bc0idc, dd,mm,yy,bcthud
      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,bcthum

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 'c1gen.inc'
      include 'm1subr.inc'
      include 'w1clig.inc'
      include 'c1info.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, y1,dayear,m
      real bsa, dg, dg1, hlmn, wt1, x1, xz, jreal,dxx
      real xx, sdmn, sdmx, hlmx,x(14),y(14),y2(14),phux,sphu,yp1,
     & ypn,bphu,ephu
      real tat, tad, ad, bd, cd,dq,d1(365,3),d2(730,3)

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 can occur
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   - running sum of heat units
c bphu   - sphu at planting date
c ephu   - sphu at harvest date
c tad    - daily heat units
c tat    - dayly average temperature (cubic spline)
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 print crop input data


c     + + + OUTPUT FORMATS + + +
 2110 format (5x,' a_co2=',f6.3,' b_co2=',f6.3,' a_frst=',f6.3,
     &' b_frst=',f6.3)
 2013 format(1x,i4,1x,a6,1x,5(f6.2,1x),1x,5(f6.1,1x),f6.3)
C**print end of season values to file annl.out (moved from crop.for: AR 7/15/00)

c     Initialize
c       initialize variables needed for season heat unit computation: added on
c       3/16/1998 (A. Retta)
        data x /-15,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.
c       For the southern hemisphere, monthly average temperatures should start
c       in July.1?
        do i=1,12
                y(i+1)=awtmav(i)
        end do
        y(1)=y(13)
        y(14)=y(2)
c       print end of season values: moved to this location: 7/15/00
c       y1=yy-1
        y1=yy
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

      if (am0cfl .gt. 0) then
          write(60,*) 'plant population (plants/m2)=', acdpop(am0csr)
      endif
c     convert to m2/plant
      parea=1./acdpop(am0csr)
      clfarea=bc0lfe/10000.
      cstarea=.1*clfarea
      pclfarea=clfarea
      pcstarea=cstarea
      slai=clfarea/parea
      ssai=cstarea/parea
c     slai=clfarea
c     ssai=cstarea
      clfwt=bc0lfe*acdpop(am0csr)/(bc0sla*1000.)
      cstwt=0.2*clfwt
      dmag=clfwt+cstwt
      rw=dmag
      prw=rw
      pdmag=dmag
      pdm=pdmag+prw
      dm=pdm
      pclfwt=clfwt
      pcstwt=cstwt
      pcrpwt=0.
      shu=0.
      hrlt=0.
      hrlty=0.
      hufy=0.
      cht=0.05*bc0lfe/3.0
      pcht=0.05
      pchty=0.0
      crpwt=0.0
      yld=0.0
      prd=0.
      prdy=0.
      cta=0.
      ceta=0.
      prcp=0.
      drmfl=0
      ctp=0.
      chaff=0.
      daye=0
      frs1=bc0fd1
      frs2=bc0fd2
      cropname=ac0nam(am0csr)
      slaix = 0.0
      ssaix = 0.0

c     write values to output file
      if (am0cfl .gt. 0) then
          write(59,2013)y1,cropname,clfwt,
     &    cstwt,crpwt,yld,dmag,cta,ceta,prcp,ctp,slaix,ssaix
      endif

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.
        ephu=0.
        bphu=0.
        tat=0.
        n=14
        yp1=0.
        ypn=0.
        phux=0.
c     day of year
        call caldatw(dd, mm, yy)
        pdate = dayear(dd, mm, yy)
c     skip autonmatic heat unit calculations if heat unit flag=1
      if (bcthud.eq.1) then
         phu=bcthum
         goto 73
      endif
c       phux = 30.0 * 9.5 / (sqrt(amalat + 1.))
c     Replace the above equation(EPIC) with the following to determine
c     the average date of planting. This equation was developed using data
c     fromthe Midwest. Its reliabilty outside of the Midwest is questionable.
c     For this to work properly in the southern hemisphere (eg. Australia)
c     the day of year should be counted starting July 1.
c     pdate=INT(279.9-7126/(amalat+1.0)+0.025*amzele+
c    &(-40.1+2.4*(bctmin)+0.29*(bctmin)**2))
c *** replace the above equation with the following (1/29/01). Under
c *** certain conditions the above equation gave negative values for
c *** pdate! L. Hagen developed the above equation, and the following
c *** relationships that are its replacement.
c     for perennial plants (require more than 1 year)
      if ((bc0idc.eq.3).or.(bc0idc.gt.5)) then
        ad=1.075+0.00347*bctmin**3
        bd=112.+1.37*bctmin**1.5
        cd=38.8-0.037*bctmin**2
        dq=3.94+0.851*bctmin**0.5
        pdate=INT(ad+bd/(1.+exp(-(amalat-cd)/dd))+0.025*amzele)
        if (pdate.le.0)pdate=1
      endif
      hdate=pdate+dtm
c       call cubic spline interpolation routines
        call spline (x, y, n, yp1, ypn, y2)
c ***debugg lines
c     write(*,*)'ad,bd,cd,dd,pdate,bctmin,amalat,amzele,dtm,dxx'
c     write(*,*)ad,bd,cd,dd,pdate,bctmin,amalat,amzele,dtm,dxx
c *** end debug lines
      do i = 1, 365
        jreal = i
c        calculate daily temps. and heat units
        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.
          d1(i,1)=i
          d1(i,2)=tad
          d2(i,1)=i
          d2(i,2)=tad
       end do
c        duplicate the first year into the second year
        do j=1,365
              m=j+365
              d2(m,1)=m
              d2(m,2)=d1(j,2)
        end do
c         running sum of heat units
        do j=1,730
            sphu=sphu+d2(j,2)
            d2(j,3)=sphu
c           if (am0cfl .gt. 0) then
c             print for debugging
c             write(60,*) d2(j,1),d2(j,2),d2(j,3)
c           end if
        end do
        sphu=0.
c            calculate seasonal heat units
        do j=1,730
           if (d2(j,1).eq.pdate)bphu=d2(j,3)
           if (d2(j,1).eq.hdate)ephu=d2(j,3)
        end do
        phu=ephu-bphu
c         adjust heat units for plants whose growth cycle is greater than 365
        dxx=dtm/(365.-FLOAT(pdate))
      if (dtm.gt.365.) phu=phu*dxx
c       increase final phu for winter wheat by 40% to account for fall growth. This was
c       adjustment is based on some test done by A. Retta (6/12/1998)
c       if (bc0idc .eq. 5) then
c         phu = phu + 0.4 * phu
c         If (phu.gt.2100.)phu=2100.0
c       endif
c ...... end of seasonal heat unit calculator algorithm.........
      if (am0cfl .gt. 0) then
          write(60,*) 'computed plant date=',pdate,
     &     ' computed harvest date=', hdate
          write(60,*)'computed total heat units (phu)=',phu
      end if
 73   continue

c     Calculate s-curve parameters
c     write(*,*)bc0be1,cc0be1
c     co2 change eliminated
c      cc0be1 = cc0be1 * 0.01
c      cc0be2 = cc0be2 * 0.01
c      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)
c     if (am0cfl .gt. 0) write (60,2110)a_co,b_co,a_fr,b_fr

c     calculate bc0bn1,bc0bn2,bc0bn3,bn4,bc0bp1,bc0bp2,bc0bp3,bp4
c      call nconc (bc0bn1,bc0bn2,bc0bn3,bn4)
c      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)
      trsd = 0.0
      tfon = 0.0
      tfop = 0.0
      tmp = 0.0
      top = 0.0
      twn = 0.0
      twmn = 0.0
      tp = 0.0
      tap = 0.0
      tno3 = 0.0
      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)
         ids = 1
         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