program bgtcfi
!
! SIMPLE OLG MODEL WITH GROWTH AND UNCERTAIN LIFETIMES
!
!    Dynamic program of individuals
!
!       V(j,a) = max  u(c,1-l) + beta ps(j,j+1)* V(j+1,a')
!               c,l,a'
!
!                s.t. (1+gz) a' = (1+i) a + w eps(j) l -c+ psi-taxes
!                                          
!       where growing variables are detrended relative to technology 
!       trend (1+gz)^t
!
!    Functional forms: 
!       u(c,1-l)   = log c + gam log(1-l)
!       F(KT,KI,L) = KT^thetT *KI^thetI * (ZL)^{1-thetT-thetI}, 2 sectors
!
!    Input file requires values for
!
!       alb        = lower bound on asset holdings 
!       aub        = upper bound on asset holdings 
!       alpha      = sector 1 share
!       beta       = discount factor
!       delt1      = depreciation rate on tangible capital, sector 1
!       deli1      = depreciation rate on intangible capital, sector 1
!       delt2      = depreciation rate on tangible capital, sector 2
!       deli2      = depreciation rate on intangible capital, sector 2
!       gam        = parameter in utility
!       gn         = growth rate of population
!       gz         = growth rate of technology
!       lam        = degree of annuitization
!       tran       = common transfers to all ages 
!       taud       = tax rate on distributions
!       taul       = tax rate on labor
!       taupc      = tax rate on corporate profits
!       taupu      = tax rate on unincorporated profits (distributed in full)
!       tfp        = TFP parameter = Z^{1-theta}
!       thett1     = tangible capital share, sector 1
!       theti1     = intangible capital share, sector 1
!       thett2     = tangible capital share, sector 2
!       theti2     = intangible capital share, sector 2
!       debt       = government borrowing 
!       gspend     = government spending
!       idebt      = 0 if debt in levels, 1 if share of GDP
!       ispend     = 0 if gspend in levels, 1 if share of GDP
!       irate      = initial guess for interest rate
!       tauc       = initial guess for the tax rate on consumption
!       J x 3 
!         eps      = labor productivity for ages j, j=1:J
!         ps       = probability of survival from j to j+1, j=1:J (Jth ignored)
!         xi       = age-dependent transfers for ages j=1:J
!       beq,jbeq   = nonaccidental bequest, age of bequester
!
!       NOTE:        edit file to change maximum age (jmax) or number of 
!                    points grid points in asset grid (n).
!
!    Output file is loaded into Matlab -- run plotbg to view results.
!
!    The algorithm used for solving the fixed point is Gauss-Seidel.
!

! Ellen McGrattan, 6-26-08
! Revised, 7-25-12

  implicit none
  integer, parameter           :: jmax=101,n=3000,npar=24,tmax=240,nx=2
  integer                      :: i,j,t,idone,it,maxit,jbeq
  integer, dimension(jmax)     :: jpi
  integer, dimension(jmax,n)   :: opt
  integer, dimension(2)        :: ipar
  real, dimension(npar)        :: par
  real, dimension(jmax)        :: cj,lj,aj,apj
  real, dimension(jmax)        :: eps,ps,xi
  real, dimension(nx)          :: x,x1,res
  real                         :: kapt1,kapt2,kapi1,kapi2,lab1,lab2,c,y,     &
                                  crit,beq,mu0,sum1
  external agg

  open(unit=5, file='bgtc.inp')
  open(unit=7, file='bgtc.nxt')
  open(unit=8, file='bg.dat')
  open(unit=9, file='tran.inp')

  !
  ! Read in parameters
  !
  do i=1,npar
    read(5,*) par(i)
  enddo
  read(5,*) ipar(1)
  read(5,*) ipar(2)

  !
  ! Read in guess for the interest rate and the tax rate on consumption
  !
  read(5,*) x(1)
  read(5,*) x(2)
  x1       = x

  !
  ! Read in labor productivities, survival probabilities, and transfers
  !
  do j=1,jmax
    read(5,*) eps(j),ps(j),xi(j)
  enddo
  ps(jmax) = 0.
  read(5,*) beq,jbeq

  ! 
  ! Compute steady state 
  !
  maxit   = 10
  crit    = 1.e-8
  idone   = 0
  it      = 0 
  opt     = -1

  !
  ! 1. Hold tauc fixed and compute irate 
  !

  do while ((idone==0).and.(it<maxit))
    call agg(x,nx,jmax,n,eps,ps,xi,par,npar,ipar,beq,jbeq,opt, res,jpi,cj,   &
             lj,aj,apj,kapt1,kapt2,kapi1,kapi2,lab1,lab2,c,y)
    x1(1) = .95*x(1)+.05*(x(1)-res(1))
    write(*,'(10(1x,e15.9))') x1,maxval(abs(x1-x)/(1+abs(x)))
    if (abs(x1(1)-x(1))/(1+abs(x(1)))<crit) idone=1
    it    = it+1
    x(1)  = x1(1)
  enddo
  !
  ! 2. Take solution for irate and compute tauc
  !
  maxit   = 10
  idone   = 0
  it      = 0 
  opt     = -1
  do while ((idone==0).and.(it<maxit))

    call agg(x,nx,jmax,n,eps,ps,xi,par,npar,ipar,beq,jbeq,opt, res,jpi,cj,   &
             lj,aj,apj,kapt1,kapt2,kapi1,kapi2,lab1,lab2,c,y)
    x1(2) = x(2)-res(2)
    write(*,'(10(1x,e15.9))') x1,maxval(abs(x1-x)/(1+abs(x)))
    if (maxval(abs(x1-x)/(1+abs(x)))<crit) idone=1
    it    = it+1
    x(2)  = x1(2)
  enddo

  if (maxit==0) then
    write(*,*) 'Results based on user-supplied inputs...'
    call agg(x1,nx,jmax,n,eps,ps,xi,par,npar,ipar,beq,jbeq,opt, res,jpi,cj,  &
             lj,aj,apj,kapt1,kapt2,kapi1,kapi2,lab1,lab2,c,y)
  endif

  !
  ! Write out new input file for steady state calculations
  !
  write(7,'(1x,f14.8,5x,a3)') par(1),'alb'
  write(7,'(1x,f14.8,5x,a3)') par(2),'aub'
  write(7,'(1x,f14.8,5x,a4)') par(3),'alpha'
  write(7,'(1x,f14.8,5x,a4)') par(4),'beta'
  write(7,'(1x,f14.8,5x,a5)') par(5),'delt1'
  write(7,'(1x,f14.8,5x,a5)') par(6),'deli1'
  write(7,'(1x,f14.8,5x,a5)') par(7),'delt2'
  write(7,'(1x,f14.8,5x,a5)') par(8),'deli2'
  write(7,'(1x,f14.8,5x,a3)') par(9),'gam'
  write(7,'(1x,f14.8,5x,a2)') par(10),'gn'
  write(7,'(1x,f14.8,5x,a2)') par(11),'gz'
  write(7,'(1x,f14.8,5x,a3)') par(12),'lam'
  write(7,'(1x,f14.8,5x,a4)') par(13),'tran'
  write(7,'(1x,f14.8,5x,a4)') par(14),'taud'
  write(7,'(1x,f14.8,5x,a4)') par(15),'taul'
  write(7,'(1x,f14.8,5x,a5)') par(16),'taupc'
  write(7,'(1x,f14.8,5x,a5)') par(17),'taupu'
  write(7,'(1x,f14.8,5x,a3)') par(18),'tfp'
  write(7,'(1x,f14.8,5x,a6)') par(19),'thett1'
  write(7,'(1x,f14.8,5x,a6)') par(20),'theti1'
  write(7,'(1x,f14.8,5x,a6)') par(21),'thett2'
  write(7,'(1x,f14.8,5x,a6)') par(22),'theti2'
  write(7,'(1x,f14.8,5x,a4)') par(23),'debt'
  write(7,'(1x,f14.8,5x,a6)') par(24),'gspend'
  write(7,'(1x,i14,5x,a5)')   ipar(1),'idebt'
  write(7,'(1x,i14,5x,a6)')   ipar(2),'ispend'
  write(7,'(1x,e14.8,5x,a5)') x1(1),'irate'
  write(7,'(1x,e14.8,5x,a4)') x1(2),'tauc'
  do j=1,jmax
    write(7,'(3(1x,e14.6))') eps(j),ps(j),xi(j)
  enddo
  write(7,'(1x,e14.8,1x,i3)') beq,jbeq

  !
  ! Write out new input file for transition calculations
  !
  write(9,'(1x,f14.8,5x,a3)') par(1),'alb'
  write(9,'(1x,f14.8,5x,a3)') par(2),'aub'
  write(9,'(1x,f14.8,5x,a5)') par(3),'alpha'
  write(9,'(1x,f14.8,5x,a4)') par(4),'beta'
  write(9,'(1x,f14.8,5x,a5)') par(5),'delt1'
  write(9,'(1x,f14.8,5x,a5)') par(6),'deli1'
  write(9,'(1x,f14.8,5x,a5)') par(7),'delt2'
  write(9,'(1x,f14.8,5x,a5)') par(8),'deli2'
  write(9,'(1x,f14.8,5x,a3)') par(9),'gam'
  write(9,'(1x,f14.8,5x,a2)') par(11),'gz'
  write(9,'(1x,f14.8,5x,a3)') par(12),'lam'
  write(9,'(1x,f14.8,5x,a6)') par(19),'thett1'
  write(9,'(1x,f14.8,5x,a6)') par(20),'theti1'
  write(9,'(1x,f14.8,5x,a6)') par(21),'thett2'
  write(9,'(1x,f14.8,5x,a6)') par(22),'theti2'
  if (ipar(1)==0) then
    write(9,'(1x,f14.8,5x,a5)') par(23),'debt0'
  else
    write(9,'(1x,f14.8,5x,a5)') par(23)*(y-((1.+par(10))*(1.+par(11))-1.+   &
         par(6))*kapi1-((1.+par(10))*(1.+par(11))-1.+par(8))*kapi2),'debt0'
  endif  
  write(9,'(1x,i14,5x,a5)')   ipar(1),'idebt'
  write(9,'(1x,i14,5x,a6)')   ipar(2),'ispend'
  mu0            = 1.
  sum1           = 0.
  do j=2,jmax
    mu0          = mu0*ps(j-1)/(1.+par(10))
    sum1         = sum1+mu0
  enddo
  mu0            = 1./(1.+sum1)
  write(9,'(2x,i13,1x,e14.6)') jpi(1),mu0
  do j=2,jmax
    mu0          = mu0*ps(j-1)/(1.+par(10))
    write(9,'(2x,i13,1x,e14.6)') jpi(j),mu0
  enddo
  do t=1,tmax
    write(9,'(9(1x,e14.6))') par(23),par(24),par(13:18),par(10)
  enddo
  do j=1,jmax
    write(9,'(1x,e14.6)') eps(j)
  enddo
  do j=1,jmax
    do t=1,tmax
      write(9,'(2(1x,e14.6))') ps(j),xi(j)
    enddo
  enddo
  do t=1,tmax
    write(9,'(3(1x,e14.6))') x1(1),(1.-par(19)-par(20))*par(3)*y/lab1,x1(2)
  enddo
  !
  ! Write out results for plotbg.m
  !
  do i=1,12
    write(8,'(f13.8)') par(i)
  enddo
  write(8,'(f13.8)') x1(2)
  do i=14,npar-2
    write(8,'(f13.8)') par(i)
  enddo
  if (ipar(1)==0) then
    write(8,'(f13.8)') par(npar-1)
  else
    write(8,'(f13.8)') par(npar-1)*(y-((1.+par(10))*(1.+par(11))-1.+par(6))* &
                       kapi1-((1.+par(10))*(1.+par(11))-1.+par(8))*kapi2)
  endif
  if (ipar(2)==0) then
    write(8,'(f13.8)') par(npar)
  else
    write(8,'(f13.8)') par(npar)*(y-((1.+par(10))*(1.+par(11))-1.+par(6))*   &
                       kapi1-((1.+par(10))*(1.+par(11))-1.+par(8))*kapi2)
  endif
  write(8,'(f13.8)') kapt1
  write(8,'(f13.8)') kapt2
  write(8,'(f13.8)') kapi1
  write(8,'(f13.8)') kapi2
  write(8,'(f13.8)') lab1
  write(8,'(f13.8)') lab2
  write(8,'(f13.8)') c
  write(8,'(f13.8)') ((1.+par(10))*(1.+par(11))-1.+par(5))*kapt1
  write(8,'(f13.8)') ((1.+par(10))*(1.+par(11))-1.+par(7))*kapt2
  write(8,'(f13.8)') ((1.+par(10))*(1.+par(11))-1.+par(6))*kapi1
  write(8,'(f13.8)') ((1.+par(10))*(1.+par(11))-1.+par(8))*kapi2
  write(8,'(f13.8)') y
  write(8,'(f13.8)') x1(1)
  write(8,'(f13.8)') par(13)
  write(8,'(i13)')   jmax
  write(8,'(i13)')   n
  do j=1,jmax
    write(8,'(f13.8)') eps(j)
  enddo 
  do j=1,jmax
    write(8,'(f13.8)') ps(j)
  enddo 
  do j=1,jmax
    write(8,'(f13.8)') xi(j)
  enddo 
  do j=1,jmax
    write(8,'(i13)') jpi(j)
  enddo
  do j=1,jmax
    write(8,'(f13.8)') cj(j)
  enddo
  do j=1,jmax
    write(8,'(f13.8)') lj(j)
  enddo
  do j=1,jmax
    write(8,'(f13.8)') aj(j)
  enddo
  do j=1,jmax
    write(8,'(f13.8)') apj(j)
  enddo
  do i=1,n
    do j=1,jmax
      write(8,'(i13)') opt(j,i)
    enddo
  enddo
  write(8,'(f13.8)') beq
  write(8,'(i13)') jbeq

end program bgtcfi

subroutine agg(x,nx,jmax,n,eps,ps,xi0,par,npar,ipar,beq,jbeq,opt, resx,jpi,  &
               cj,lj,aj,apj,kapt1,kapt2,kapi1,kapi2,lab1,lab2,c,y)
  !
  ! Dynamic programs for the OLG model
  !
  implicit none
  integer, intent(in)                      :: jmax,n,npar,jbeq,nx
  integer                                  :: i,j,k,k1,k2,kzero
  integer, dimension(1)                    :: kk
  integer, dimension(2),intent(in)         :: ipar
  integer, dimension(jmax),intent(out)     :: jpi
  integer, dimension(jmax,n),intent(inout) :: opt
  real, dimension(nx),intent(in)           :: x
  real, dimension(nx),intent(out)          :: resx
  real, dimension(npar),intent(in)         :: par
  real, intent(in)                         :: beq
  real, dimension(jmax),intent(in)         :: eps,ps,xi0
  real, dimension(jmax),intent(out)        :: cj,lj,aj,apj
  real, intent(out)                        :: kapt1,kapt2,kapi1,kapi2,lab1,  &
                                              lab2,c,y
  real, dimension(jmax)                    :: mu,xi
  real, dimension(n)                       :: vnext,v,rhs
  real, dimension(4)                       :: input
  real, dimension(3)                       :: output
  real                                     :: alb,aub,alpha,alph1,beta,      &
                                              delt1,deli1,delt2,deli2,gn,    &
                                              gz,gam,tfp,tauc,taud,taul,     &
                                              taupc,taupu,thett1,theti1,     &
                                              thetl1,thett2,theti2,thetl2,   &
                                              w,sum1,mu0,tran,tkap,a,at,atw, &
                                              irate,cdiv,cprof,uprof,asset,  &
                                              dshare,gshare,gdp,rt1,rt2,ri1, &
                                              ri2,tranj,invt1,invt2,invi1,   &
                                              invi2,p1,p2,y1,y2,krat1,krat2, &
                                              kaprat,kdivlt1,kdivlt2,        &
                                              kdivli1,kdivli2,ydivl,ldivl1,  &
                                              ldivl2,lab,wtmp,coef1,lam
  
  !
  ! Initialize parameters
  !
  alb            = par(1)
  aub            = par(2)
  alpha          = par(3)
  beta           = par(4)
  delt1          = par(5) 
  deli1          = par(6) 
  delt2          = par(7) 
  deli2          = par(8) 
  gam            = par(9)
  gn             = par(10)
  gz             = par(11)
  lam            = par(12)
  tran           = par(13)
  taud           = par(14)
  taul           = par(15)
  taupc          = par(16)
  taupu          = par(17)
  tfp            = par(18)
  thett1         = par(19)
  theti1         = par(20)
  thett2         = par(21)
  theti2         = par(22)
  alph1          = 1.-alpha
  thetl1         = 1.-thett1-theti1
  thetl2         = 1.-thett2-theti2
  irate          = x(1)
  tauc           = x(2)
  rt1            = (irate-1.)/(1.-taupc)+delt1
  rt2            = (irate-1.)/(1.-taupu)+delt2
  ri1            = irate-1.+deli1
  ri2            = irate-1.+deli2
  ldivl1         = thetl1*alpha/(thetl1*alpha+thetl2*alph1)
  ldivl2         = thetl2*alph1/(thetl1*alpha+thetl2*alph1)
  kdivlt1        = thett1/(thetl1*rt1)
  kdivlt2        = thett2/(thetl2*rt2)
  kdivli1        = theti1/(thetl1*ri1)
  kdivli2        = theti2/(thetl2*ri2)
  ydivl          = 2.*(tfp*kdivlt1**thett1*kdivli1**theti1*ldivl1)**alpha*   &
                      (tfp*kdivlt2**thett2*kdivli2**theti2*ldivl2)**alph1
  w              = ((thetl1*alpha+thetl2*alph1)*ydivl)**                     &
                   (1./(1.-((1.-thetl1)*alpha+(1.-thetl2)*alph1)))
  kdivlt1        = kdivlt1*w
  kdivlt2        = kdivlt2*w
  kdivli1        = kdivli1*w
  kdivli2        = kdivli2*w
  ydivl          = ydivl*w**((1.-thetl1)*alpha+(1.-thetl2)*alph1)
  atw            = (1.-taul)*w
  mu0            = 1.
  sum1           = 0.
  do j=2,jmax
    mu0          = mu0*ps(j-1)/(1.+gn)
    sum1         = sum1+mu0
  enddo
  mu(1)          = 1./(1.+sum1)
  do j=2,jmax
    mu(j)        = mu(j-1)*ps(j-1)/(1.+gn)
  enddo
  xi             = xi0
  xi(jbeq)       = xi(jbeq) - beq
  xi(1)          = xi(1) +mu(jbeq)*beq/mu(1)
  kzero          = (aub-alb*n)/(aub-alb)
  input(3)       = gam
  input(4)       = 1.+tauc


  !
  ! Iterate backwards for ages j=jmax to 1, solving
  ! dynamic programs recursively 
  !
  input(2)       = atw*eps(jmax)
  do i=1,n
    a            = (aub*float(i-1)+alb*float(n-i))/float(n-1)
    input(1)     = tran+xi(jmax)+irate*a
    call static(input, output)
    vnext(i)     = output(1)
  enddo
  opt(jmax,:)    = kzero

  do j=jmax-1,1,-1
    input(2)     = atw*eps(j)
    do i=1,n
      a          = (aub*float(i-1)+alb*float(n-i))/float(n-1)
      if (opt(j,i)>0) then
        k1       = max(1,floor(float(opt(j,i))-max(2.,.1*n)))
        k2       = min(floor(float(opt(j,i))+max(2.,.1*n)),n)
      else
        k1       = 1
        k2       = n
      endif
      rhs        = -1.e+9
      do k=k1,k2
        at       = (aub*float(k-1)+alb*float(n-k))/float(n-1)
        input(1) = tran+xi(j)+irate*a-(1.+gz)*at *(1.-lam*(1.-ps(j)))
        call static(input, output)
        rhs(k)   = output(1)+beta*ps(j)*vnext(k)
      enddo
      kk         = maxloc(rhs)
      k          = kk(1)
      opt(j,i)   = k
      v(i)       = rhs(k)
    enddo
    vnext        = v
  enddo
  ! 
  ! Update distributions of asset holdings
  !
  jpi(1)         = kzero
  do j=1,jmax-1
    jpi(j+1)     = opt(j,jpi(j))
  enddo
  !
  ! Add everything up 
  !
  asset          = 0.
  tkap           = 0.
  lab            = 0.
  c              = 0.
  y              = 0.
  cj             = 0.
  lj             = 0.
  aj             = 0.
  apj            = 0.

  do j=1,jmax
    input(2)     = atw*eps(j)
    i            = jpi(j)
    a            = (aub*float(i-1)+alb*float(n-i))/float(n-1)
    if (j<jmax) then
      k          = jpi(j+1)
      at         = (aub*float(k-1)+alb*float(n-k))/float(n-1)
    else
      at         = 0.
    endif
    input(1)     = tran+xi(j)+irate*a-(1.+gz)*at *(1.-lam*(1.-ps(j)))
    call static(input,output)
    asset        = asset+ mu(j)*a
    tkap         = tkap + mu(j)*(1.-ps(j))*(1.-lam)*at
    c            = c    + mu(j)*output(2)
    lab          = lab  + mu(j)*eps(j)*output(3)
    cj(j)        = cj(j)+ output(2)
    lj(j)        = lj(j)+ output(3)
    aj(j)        = aj(j)+ a
    apj(j)       = apj(j)+at
  enddo
  tranj          = sum(mu*xi)
  y              = ydivl*lab
  krat1          = kdivli1/kdivlt1
  krat2          = kdivli2/kdivlt2
  kaprat         = (kdivlt2/kdivlt1)*(ldivl2/ldivl1)
  coef1          = (1.-taud)*(1.+(1.-taupc)*krat1)+(1.+(1.-taupu)*krat2)*kaprat
  if (ipar(1)==0) then
    kapt1        = (asset-par(23))/coef1
    kapt2        = kaprat*kapt1
    kapi1        = krat1*kapt1
    kapi2        = krat2*kapt2
    invi1        = ((1.+gz)*(1.+gn)-1.+deli1)*kapi1
    invi2        = ((1.+gz)*(1.+gn)-1.+deli2)*kapi2
    gdp          = y-invi1-invi2
    dshare       = par(23)/gdp
  else
    dshare       = par(23)
    kapt1        = (asset-dshare*y)/(coef1-dshare*((1.+gz)*(1.+gn)-1.+       &
                   deli1)*krat1-dshare*((1.+gz)*(1.+gn)-1.+deli2)*krat2*     &
                   kaprat) 
    kapt2        = kaprat*kapt1
    kapi1        = krat1*kapt1
    kapi2        = krat2*kapt2
    invi1        = ((1.+gz)*(1.+gn)-1.+deli1)*kapi1
    invi2        = ((1.+gz)*(1.+gn)-1.+deli2)*kapi2
    gdp          = y-invi1-invi2
  endif
  invt1          = ((1.+gz)*(1.+gn)-1.+delt1)*kapt1
  invt2          = ((1.+gz)*(1.+gn)-1.+delt2)*kapt2
  lab1           = ldivl1*lab
  lab2           = ldivl2*lab
  y1             = tfp*(kapt1)**thett1*(kapi1)**theti1*(lab1)**thetl1
  y2             = tfp*(kapt2)**thett2*(kapi2)**theti2*(lab2)**thetl2
  p1             = alpha*y/y1
  p2             = alph1*y/y2
  cprof          = p1*y1-w*lab1-delt1*kapt1-invi1
  uprof          = p2*y2-w*lab2-delt2*kapt2-invi2
  cdiv           = p1*y1-w*lab1-invt1-invi1-taupc*cprof

  !
  ! Update guess for capital/labor ratio and transfers
  !
  if (ipar(2)==0) then
    gshare       = par(24)/gdp
  else
    gshare       = par(24)
  endif
  resx(1)        = rt1-thett1*alpha*y/kapt1
  resx(2)        = tauc-(tran+tranj+gshare*gdp-(1.+gz)*tkap-taul*w*lab-      &
                   taud*cdiv-taupc*cprof-taupu*uprof-((1.+gz)*(1.+gn)-       &
                   irate)*dshare*gdp)/c
end subroutine agg

subroutine static(input, output)
  !
  ! Two equations in two unknowns (c,l)
  !    (du/dl)/(du/dc) = w
  !    c = (nlinc+w*l)/pc 
  !
  ! where
  !    u(c,1-l) = log(c)+gam*log(1-l)
  !    nlinc = nonlabor income
  !    l in [0,1] is checked

  implicit none
  real,dimension(4),intent(in)   :: input
  real,dimension(3),intent(out)  :: output
  real                           :: nlinc,w,gam,pc,u,c,l

  nlinc      = input(1)
  w          = input(2)
  gam        = input(3)
  pc         = input(4)
  if (w<1.0e-10) then
    l        = 0.
    c        = nlinc/pc
  else
    l        = (w-gam*nlinc)/(w+w*gam)
    c        = (nlinc+w*l)/pc
  endif
  if (l<=0.) then
    l        = 0.
    c        = nlinc/pc
  endif
  if (l>=1.) then
    l        = 1.
    c        = (nlinc+w)/pc
  endif
  if (c<=0.) then
    u        = -1.e+9
    c        = 0.
  else
    u        = log(c)+gam*log(1.-l)
  endif
  output     = (/ u,c,l /)

end subroutine static