text:
!
!   SimulateBPP.f90
!   Compute insurance coefficients a la Blundell, Pistaferri and Preston (2008)
!
!   Created by Kurt Mitman on 11/08/15.
subroutine SimulateBPP

use params
use globals
use funcs
use procedures

implicit none

INTEGER		:: iseed(2),isize,iAgg,it,Phinds(2),iindex

!Things that don't need to be private
double precision    :: Ph,Pr
double precision                    :: vftol
double precision    :: prandtemp
double precision    :: phptemp(ngpAgg),ephtemp
!Variables that need to be private
INTEGER		:: in,ij,ijret,iW,iR,polind,hjI,index,mcc,hcc,lcc,ljI,ip,ify,izy,iPh,mc
double precision    :: bccmax,ftemp1,ftemp2,bjps,coef1,coef2,coef3
integer     :: inds(2)
double precision    ::  vals(2)
double precision    :: helpWsell,helpWfore
double precision    :: HelpVals(5)
double precision    :: btemp, y, pm, fc
double precision    :: bj,bjp,mj,mjp,hj,lj
double precision    :: Wrentsim,consRpolsim,hconsRpolsim,bRpolsim
double precision    :: Wbuysim
double precision    :: Wrefisim,consNpolsim,hconsNpolsim,bNpolsim,mNpolsim,hNpolsim,lNpolsim
double precision    :: Wpaysim,consPpolsim,hconsPpolsim,bPpolsim,mPpolsim,hPpolsim,lPpolsim
double precision,dimension(nh*nm*nl)    :: helpWB,consB,bB,mPp
integer,dimension(nh*nm*nl) ::hB,mB,hconsB,lB
integer             :: hINpolsim,hIPpolsim,lINpolsim,lIPpolsim
double precision    :: mpcval,bjpmpc,bjmpc,bccmid,bccmin

double precision    :: phh,qmtemp,qltemp,qltemp2,RentVal,BuyVal,Phvals(2)

double precision    :: Ht,htrans

DOUBLE PRECISION, EXTERNAL :: FnVFR, golden, FnVFB, FnVFF, FnVFN,FnVFP,FnXSH
DOUBLE PRECISION, EXTERNAL          :: FnVFRRetJ, FnVFFRetJ, FnVFPRetJ
DOUBLE PRECISION, EXTERNAL          :: FnVFRRet, FnVFBRet, FnVFFRet, FnVFNRet,FnVFPRet

double precision,dimension(ngpPh)   ::Hdemand,Hsupply
double PRECISION, allocatable, dimension(:,:)::hsim,msim,lsim,bsim,hconssim,csim,beqsim
INTEGER, allocatable, dimension(:,:)::hsimI,buysim,ownsim,lsimI,refisim,foresim,sellsim
double precision     :: HInv(ngpAy,ngpZh,ngpPh)
double precision     :: Vmin,Wmin
integer              ::idsimI(nsim)
allocate(hsim(ngpPh,nsim))
allocate(msim(ngpPh,nsim))
allocate(lsim(ngpPh,nsim))
allocate(bsim(ngpPh,nsim))
allocate(hconssim(ngpPh,nsim))
allocate(csim(ngpPh,nsim))
allocate(beqsim(ngpPh,nsim))
allocate(hsimI(ngpPh,nsim))
allocate(lsimI(ngpPh,nsim))
allocate(buysim(ngpPh,nsim))
allocate(ownsim(ngpPh,nsim))
allocate(refisim(ngpPh,nsim))
allocate(foresim(ngpPh,nsim))
allocate(sellsim(ngpPh,nsim))

Vmin=-1.0d10
Wmin=Vmin

Ht=Hsteady
phh=ssph(iAggsteady)
print*,Hsteady,ssph(iAggsteady)
ownsim=0
hsim=0.0d0
hconssim=0.0d0
bsim=0.0d0

vftol=1.0d-5

!$OMP PARALLEL DO
DO in = 1,nsim

    idsimI(in) = mod(in-1,nsimss)+1
    simtrans(tstart)%jsim(in) = (in-1)/nsimss+1
    if(simtrans(tstart+1)%jsim(in) .LT. Jtot) then
       simtrans(tstart+1)%hsimI(in)=1
       simtrans(tstart+1)%lsimI(in)=1
       simtrans(tstart+1)%lsim(in)=0
       simtrans(tstart+1)%hsim(in)=0
       simtrans(tstart+1)%msim(in)=0
       simtrans(tstart+1)%bsim(in)=1.0d0
       simtrans(tstart+1)%ownsim(in)=0
    endif
    simtrans(tstart)%csim(in)=0.35d0
    simtrans(tstart)%hconssim(in)=1.0d0
    simtrans(tstart)%sellsim(in)=0
    simtrans(tstart)%refisim(in)=0
    simtrans(tstart)%foresim(in)=0
    simtrans(tstart)%buysim(in)=0
    fysimtransI(in)=fysimI(in)
    simtrans(tstart)%thetasim(in)=2

    simtrans(tstart)%zysimI(in)=max(zysimI(idsimI(in),simtrans(tstart)%jsim(in)),1)
    simtrans(tstart)%psimI(in) = (simtrans(tstart)%zysimI(in)-1)*ngpfy+fysimtransI(in)

END DO
!$OMP END PARALLEL DO

do it=tstart+1,tend-1

    iAgg=iAggpath(it)
    giAgg=iAgg
    AYpath(it)=GridAY(AtoY(iAgg))
    iAYpath(it)=AtoY(iAgg)
    HDpath(it)=GridDemand(AtoD(iAgg))
    iHDpath(it)=AtoD(iAgg)
    giHD=iHDpath(it)
    giAY=iAYpath(it)
    giZh=AtoZ(iAgg)
    giC=AtoC(iAgg)
    giRf = AtoR(giAgg)
    Ht=Htpath(it)
    Ph=Phpath(it)
    gPh=Ph
    giRf = AtoR(giAgg)
    !Changed 11/20/15 to make grid on Rf on Agg for comovement
    rf=RfGrid(giAgg)
    grf=rf
    rm=Rmgrid(giAgg)
    grm=rm
    rl=RlGrid(giAgg)
    grl=rl
    htrans=gPh*(1.0d0-alpha_h)*(gPh**(alpha_h/(1.0d0-alpha_h)))*((GridZh(giZh))**(1.0d0/(1.0d0-alpha_h)))*((alpha_h)**(alpha_h/(1.0d0-alpha_h)))

    call basefun(Phgrid,ngpPh,Ph,Phvals,Phinds)

    !Simulate each guy

    !$OMP PARALLEL DO PRIVATE(ij,ijret,iW,iR,polind,hjI,index,mcc,mc,hcc,lcc,ljI,ip,ify,izy,bccmax,ftemp1,ftemp2,bjps,inds,vals,helpWsell,helpWfore,HelpVals,btemp, y, pm, fc,bj,bjp,mj,mjp,hj,lj,Wrentsim,consRpolsim,hconsRpolsim,bRpolsim,Wbuysim,Wrefisim,consNpolsim,hconsNpolsim,bNpolsim,mNpolsim,hNpolsim,lNpolsim,Wpaysim,consPpolsim,hconsPpolsim,bPpolsim,mPpolsim,hPpolsim,lPpolsim,helpWB,consB,bB,hB,mB,hconsB,lB,hINpolsim,hIPpolsim,lINpolsim,lIPpolsim,mPp,qmtemp,qltemp,qltemp2,RentVal,BuyVal,coef1,coef2,coef3,bjmpc,mpcval,bjpmpc,bccmin,bccmid) COPYIN(giAy,giHD,giAgg,giC,giZh,giPh,giRf, gPh,grl,grf,grm)

    do in=1,nsim

        IF(simtrans(it-1)%jsim(in) .eq. Jtot) THEN  !die with prob 1
            simtrans(it-1)%diesimI(in) = 1
            simtrans(it)%jsim(in) = 1
        ELSE
            simtrans(it-1)%diesimI(in) = 0
            simtrans(it)%jsim(in) = simtrans(it-1)%jsim(in) + 1
        END IF
        ij = simtrans(it)%jsim(in)
        ijret = ij-Jwork
        rf=grf
        rm=grm
        rl=grl
        !Use $500 to compute the MPC for households
        mpcval = 500.0d0/(DataAvAnnualEarns/Numeraire)

        simtrans(it)%eysimI(in)= 1

        select case(ij)

            case(1)
                CALL DiscreteDist1(simtrans(it)%zysimI(in), ngpzy, zydist(1,:),zyrand(in,it))
                simtrans(it)%ypssim(in) = ypsgrid(ij,fysimtransI(in),simtrans(it)%zysimI(in),giAy)
                simtrans(it)%psimI(in) = 1
                simtrans(it)%psim(in) = 0.0
                simtrans(it)%ysim(in) = ygrid(ij,fysimtransI(in),simtrans(it)%zysimI(in),simtrans(it)%eysimI(in),giAy)
                simtrans(it)%hsimI(in)=1
                simtrans(it)%lsimI(in)=1
                simtrans(it)%ownsim(in)=0
                simtrans(it)%bsim(in)=0
                simtrans(it)%hsim(in)=0
                simtrans(it)%lsim(in)=0
                simtrans(it)%msim(in)=0

            case(2:Jwork)
                CALL DiscreteDist1(simtrans(it)%zysimI(in), ngpzy, zytrans(ij-1,simtrans(it-1)%zysimI(in),:),zyrand(in,it))
                CALL DiscreteDist1(simtrans(it)%thetasim(in), ngptheta, (/thet(ij), 1.0d0-thet(ij) /), thetarand(in,it))
                simtrans(it)%ysim(in) = ygrid(ij,fysimtransI(in),simtrans(it)%zysimI(in),simtrans(it)%eysimI(in),giAy)
                simtrans(it)%ypssim(in) = ypsgrid(ij,fysimtransI(in),simtrans(it)%zysimI(in),giAy)
                simtrans(it)%psimI(in) = 1
                simtrans(it)%psim(in) = 0.0d0

            case(Jwork+1)
               simtrans(it)%psimI(in) = (simtrans(it-1)%zysimI(in)-1)*ngpfy+fysimtransI(in)
               simtrans(it)%psim(in) = pgrid(ij,simtrans(it)%psimI(in))
               simtrans(it)%ysim(in) = 0.0d0
               simtrans(it)%ypssim(in) = 0.0d0

            case(Jwork+2:Jtot)
               simtrans(it)%ysim(in) = 0.0d0
               simtrans(it)%ypssim(in) = 0.0d0
               simtrans(it)%psimI(in) = simtrans(it-1)%psimI(in)
               simtrans(it)%psim(in) = pgrid(ij,simtrans(it)%psimI(in))

        end select

        ip = simtrans(it)%psimI(in)
        ify = fysimtransI(in)
        izy = simtrans(it)%zysimI(in)
        bj=simtrans(it)%bsim(in)

        !Increase cah by $500 to calculate the MPC
        bjmpc=bj+mpcval

        select case(ij)

            case(1:Jwork)

                gij=ij

                iW=Wind(ify,izy)
                giW=iW
                y=simtrans(it)%ysim(in)
                giExo=(giAgg-1)*ngpW+giW

                !Check if comes into period with house
                if(simtrans(it)%ownsim(in) .eq. 0) then

                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wrent(:,giExo,:,gij),bj,gPh,RentVal)
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wbuy(:,giExo,:,gij),bj,gPh,BuyVal)

                    !Check to see if stay renter or buy house
                    if( RentVal .gt. BuyVal  ) then

                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,bRpol(:,giExo,:,gij),bj,gPh,simtrans(it+1)%bsim(in))
                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpol(:,giExo,:,gij),bj,gPh,simtrans(it)%csim(in))
                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,hconsRpol(:,giExo,:,gij),bj,gPh,simtrans(it+1)%hconssim(in))

                        simtrans(it)%buysim(in)=0
                        simtrans(it+1)%ownsim(in)=0
                        simtrans(it+1)%hsim(in)=0.0d0
                        simtrans(it+1)%msim(in)=0.0d0
                        simtrans(it+1)%lsim(in)=0.0d0
                        simtrans(it)%refisim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%sellsim(in)=0
                        simtrans(it+1)%hsimI(in)=1
                        simtrans(it+1)%lsimI(in)=1

                    else
                        simtrans(it)%buysim(in)=1
                        simtrans(it+1)%ownsim(in)=1
                        simtrans(it)%refisim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%sellsim(in)=0

                        EVown(:,:,:,:,gij)=Phvals(1)*EEVown(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVown(:,:,:,:,gij,giExo,Phinds(2))

                        do hcc=1,nh
                            fc=0.0d0
                            do mcc=1,nm
                                do lcc=1,nl
                                    index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                    mB(index)=mcc
                                    hB(index)=hcc
                                    lB(index)=lcc
                                    gih=hcc
                                    gil=lcc
                                    gim=mcc
                                    ghouse=GridH(gih)
                                    gmort=GridM(gim)
                                    gloc=GridL(gil)
                                    gliq=bj+y-FnTax(y)-gPh*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans

                                    bccmax=(gliq+gmort*ghouse)*(1.0d0+rf)
                                    !Changed 10/13/15 to prevent new home buyers from using HELOCs for purchase
                                    if(giC .eq. ngpC) then
                                       bccmin=0.0d0
                                    elseif(HELOCVal .eq. 1) then
                                       bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                    else
                                       bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                    endif

                                    bccmid=(bccmax+bccmin)/2.0d0

                                    if(bccmax .gt. bccmin) then

                                        ftemp1=FnVFB(0.0d0)
                                        ftemp2=golden(bccmin,bccmid,bccmax,FnVFB,vftol,btemp)
                                        if(ftemp2<ftemp1) then
                                            ftemp1=ftemp2
                                        else
                                            btemp=0.0d0
                                        endif
                                    else
                                       helpWB(index)=Wmin
                                       btemp=0.0d0
                                       consB(index)=cmin
                                    endif

                                    bB(index)=btemp
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qm(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,ql(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)
                                    consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                    helpWB(index)=-ftemp1
                                enddo
                                fc=FCGrid(giAgg) !Set fixed cost of mortgage for when mcc>1
                            enddo
                        enddo
                        polind=maxloc(helpWB,dim=1)

                        simtrans(it+1)%hsim(in)=GridH(hB(polind))
                        simtrans(it+1)%hsimI(in)=hB(polind)
                        simtrans(it+1)%msim(in)=GridM(mB(polind))
                        simtrans(it+1)%bsim(in)=bB(polind)
                        simtrans(it+1)%lsim(in)=GridL(lB(polind))
                        simtrans(it+1)%lsimI(in)=lB(polind)
                        simtrans(it)%hconssim(in)=GridH(hB(polind))
                        simtrans(it)%csim(in)=consB(polind)

                    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    !Now do calculations to get MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wrent(:,giExo,:,gij),bjmpc,gPh,RentVal)
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wbuy(:,giExo,:,gij),bjmpc,gPh,BuyVal)

                    !Check to see if stay renter or buy house
                    if( RentVal .gt. BuyVal  ) then

                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpol(:,giExo,:,gij),bjmpc,gPh,simtrans(it)%csim2(in))

                    else

                        EVown(:,:,:,:,gij)=Phvals(1)*EEVown(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVown(:,:,:,:,gij,giExo,Phinds(2))

                        do hcc=1,nh
                            fc=0.0d0
                            do mcc=1,nm
                                do lcc=1,nl
                                    index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                    mB(index)=mcc
                                    hB(index)=hcc
                                    lB(index)=lcc
                                    gih=hcc
                                    gil=lcc
                                    gim=mcc
                                    ghouse=GridH(gih)
                                    gmort=GridM(gim)
                                    gloc=GridL(gil)
                                    gliq=bjmpc+y-FnTax(y)-gPh*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans

                                    bccmax=(gliq+gmort*ghouse)*(1.0d0+rf)
                                    !Changed 10/13/15 to prevent new home buyers from using HELOCs for purchase
                                    if(giC .eq. ngpC) then
                                       bccmin=0.0d0
                                    elseif(HELOCVal .eq. 1) then
                                       bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                    else
                                       bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                    endif

                                    bccmid=(bccmax+bccmin)/2.0d0

                                    if(bccmax .gt. bccmin) then

                                        ftemp1=FnVFB(0.0d0)
                                        ftemp2=golden(bccmin,bccmid,bccmax,FnVFB,vftol,btemp)
                                        if(ftemp2<ftemp1) then
                                            ftemp1=ftemp2
                                        else
                                            btemp=0.0d0
                                        endif
                                    else
                                       helpWB(index)=Wmin
                                       btemp=0.0d0
                                       consB(index)=cmin
                                    endif
                                    bB(index)=btemp
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qm(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,ql(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)
                                    consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                    helpWB(index)=-ftemp1
                                enddo
                                fc=FCGrid(giAgg) !Set fixed cost of mortgage for when mcc>1
                            enddo
                        enddo
                        polind=maxloc(helpWB,dim=1)

                        simtrans(it)%csim2(in)=consB(polind)

                    endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    !END MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                else ! Homeowner

                    if(Modify .eq. 1) then

                        if(it .eq. 166) simtrans(it)%msim(in)=min(simtrans(it)%msim(in),0.95*PhPath(106))
                        if(it .eq. 206) simtrans(it)%msim(in)=min(simtrans(it)%msim(in),0.90*PhPath(106))

                    endif

                    hj=simtrans(it)%hsim(in)
                    lj=simtrans(it)%lsim(in)*(1.0d0+rl)
                    ljI=simtrans(it)%lsimI(in)
                    mj=simtrans(it)%msim(in)*(1.0d0+rm)
                    hjI=simtrans(it)%hsimI(in)

                    bjp=simtrans(it)%bsim(in)+hj*(Ph-mj-lj)-FnAdj(Ph*hj)

                    !Compared to those who started off as renters, additional income from tax deduction
                    bjps=bjp+(Fntax(y)-FnTaxM(y,hj*(simtrans(it)%msim(in)+simtrans(it)%lsim(in)),min(0.0d0,bj)))

                    if(bjps .lt. GridBS(1)) then
                        Wrentsim=Vmin
                        Wbuysim=Vmin
                    elseif(bjps .lt. 0.0d0) then
                        !Renting value

                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WrentS(:,giExo,:,gij),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WbuyS(:,giExo,:,gij),bjps,gPh,Wbuysim)
                    else
                        !Renting value

                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wrent(:,giExo,:,gij),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wbuy(:,giExo,:,gij),bjps,gPh,Wbuysim)
                    endif

                    helpWsell=max(Wbuysim,Wrentsim)

                    ! Figure out foreclosure value
                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wfore(:,ljI,nm,hjI,giExo,:,gij),bj,gPh,helpWfore)
                    !First check if the moving shock hits

                    if(simtrans(it)%thetasim(in) .eq. 1) then !Hit by moving shock

                        !Only options are sell or foreclosure

                        if(helpWsell .ge. helpWfore) then !If selling is better

                            simtrans(it)%sellsim(in)=1
                            simtrans(it)%foresim(in)=0
                            if(Wrentsim .gt. Wbuysim) then !If renting is better

                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,bRpol(:,giExo,:,gij),bjps,gPh,simtrans(it+1)%bsim(in))
                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpol(:,giExo,:,gij),bjps,gPh,simtrans(it)%csim(in))
                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,hconsRpol(:,giExo,:,gij),bjps,gPh,simtrans(it+1)%hconssim(in))

                                simtrans(it)%buysim(in)=0
                                simtrans(it+1)%ownsim(in)=0
                                simtrans(it+1)%hsim(in)=0.0d0
                                simtrans(it+1)%msim(in)=0.0d0
                                simtrans(it+1)%lsim(in)=0.0d0
                                simtrans(it)%refisim(in)=0

                            else

                                EVown(:,:,:,:,gij)=Phvals(1)*EEVown(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVown(:,:,:,:,gij,giExo,Phinds(2))

                                do hcc=1,nh
                                    fc=0.0d0
                                    do mcc=1,nm
                                    do lcc=1,nl
                                        index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                        mB(index)=mcc
                                        hB(index)=hcc
                                        lB(index)=lcc
                                        gih=hcc
                                        gil=lcc
                                        gim=mcc
                                        ghouse=GridH(gih)
                                        gmort=GridM(gim)
                                        gloc=GridL(gil)
                                        gliq=bjps+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans

                                        bccmax=(gliq+gmort*ghouse)*(1.0d0+rf)
                                        !Changed 10/13/15 to prevent new home buyers from using HELOCs for purchase
                                        if(giC .eq. ngpC) then
                                           bccmin=0.0d0
                                        elseif(HELOCVal .eq. 1) then
                                           bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                        else
                                           bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                        endif

                                        bccmid=(bccmax+bccmin)/2.0d0

                                        if(bccmax .gt. bccmin) then

                                            ftemp1=FnVFB(0.0d0)
                                            ftemp2=golden(bccmin,bccmid,bccmax,FnVFB,vftol,btemp)
                                            if(ftemp2<ftemp1) then
                                                ftemp1=ftemp2
                                            else
                                                btemp=0.0d0
                                            endif
                                        else
                                           helpWB(index)=Wmin
                                           btemp=0.0d0
                                           consB(index)=cmin
                                        endif

                                        bB(index)=btemp

                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qm(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,ql(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)
                                        consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                        helpWB(index)=-ftemp1
                                    enddo
                                    fc=mortfc !Set fixed cost of mortgage for when mcc>1
                                    enddo
                                enddo
                                polind=maxloc(helpWB(1:nm*nh*nl),dim=1)

                                simtrans(it+1)%ownsim(in) = 1
                                simtrans(it)%buysim(in) = 1
                                simtrans(it+1)%hsim(in)=GridH(hB(polind))
                                simtrans(it+1)%hsimI(in)=hB(polind)
                                simtrans(it+1)%msim(in)=GridM(mB(polind))
                                simtrans(it+1)%bsim(in)=bB(polind)
                                simtrans(it+1)%lsim(in)=GridL(lB(polind))
                                simtrans(it+1)%lsimI(in)=lB(polind)
                                simtrans(it)%hconssim(in)=simtrans(it+1)%hsim(in)
                                simtrans(it)%csim(in)=consB(polind)
                             endif

                        else  !If foreclosing is better

                            simtrans(it)%buysim(in) = 0
                            simtrans(it+1)%ownsim(in) = 0
                            simtrans(it)%foresim(in) = 1
                            simtrans(it)%sellsim(in)=0
                            simtrans(it)%refisim(in)=0

                            simtrans(it)%hconssim(in) = simtrans(it)%hsim(in)
                            call BiLinInterp1(nb,GridB,ngpPh,PhGrid,bFpol(:,ljI,nm,hjI,giExo,:,gij),bj,gPh,simtrans(it+1)%bsim(in))
                            simtrans(it+1)%bsim(in)=max(0.0d0,simtrans(it+1)%bsim(in))
                            simtrans(it)%csim(in) = max(cmin,bj+y-FnTax(y)-FnQb(simtrans(it+1)%bsim(in),0.0d0)*simtrans(it+1)%bsim(in))
                            simtrans(it+1)%hsim(in) = 0.0d0
                            simtrans(it+1)%hsimI(in) = 1
                            simtrans(it+1)%msim(in) = 0.0d0
                            simtrans(it+1)%lsim(in) = 0.0d0
                            simtrans(it+1)%lsimI(in)=1

                        endif

                    else !Not hit by moving shock

                        if(mj .eq. 0.0d0) then
                            mc=1
                            call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wpay(:,ljI,mc,hjI,giExo,:,gij),bj,gPh,Wpaysim)
                            call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wrefin(:,ljI,mc,hjI,giExo,:,gij),bj,gPh,Wrefisim)

                        else

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,Wpay(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,Wpaysim)
                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,Wrefin(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,Wrefisim)

                        endif

                        helpVals(1) =   Wbuysim
                        helpVals(2) =   Wrentsim
                        helpVals(3) =   helpWfore
                        helpVals(4) =   Wpaysim
                        helpVals(5) =   Wrefisim

                        polind=maxloc(helpVals,dim=1)
                        select case(polind)

                            case(1) !Buying optimal
                                EVown(:,:,:,:,gij)=Phvals(1)*EEVown(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVown(:,:,:,:,gij,giExo,Phinds(2))

                                simtrans(it)%sellsim(in)=1
                                simtrans(it)%foresim(in)=0
                                simtrans(it)%buysim(in)=1
                                simtrans(it+1)%ownsim(in)=1
                                do hcc=1,nh
                                    fc=0.0d0
                                    do mcc=1,nm
                                    do lcc=1,nl
                                        index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                        mB(index)=mcc
                                        hB(index)=hcc
                                        lB(index)=lcc
                                        gih=hcc
                                        gil=lcc
                                        gim=mcc
                                        ghouse=GridH(gih)
                                        gmort=GridM(gim)
                                        gloc=GridL(gil)
                                        gliq=bjps+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))
                                        ftemp1=FnVFB(bmin)
                                        ftemp2=FnVFB(GridB(2))
                                        bccmax=gliq-fc
                                        if(ftemp2<ftemp1) then
                                            ftemp1=golden(bmin,GridB(2),bccmax,FnVFB,vftol,btemp)
                                        else
                                            btemp=bmin
                                        endif
                                        bB(index)=btemp
                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qm(:,gil,gim,gih,giExo,:,gij+1),btemp,Ph,qmtemp)
                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,ql(:,gil,gim,gih,giExo,:,gij+1),btemp,Ph,qltemp)

                                        consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                        helpWB(index)=-ftemp1
                                    enddo
                                    fc=mortfc !Set fixed cost of mortgage for when mcc>1
                                    enddo
                                enddo
                                polind=maxloc(helpWB(1:nh*nl*nm),dim=1)

                                simtrans(it+1)%hsim(in)=GridH(hB(polind))
                                simtrans(it+1)%hsimI(in)=hB(polind)
                                simtrans(it+1)%msim(in)=GridM(mB(polind))
                                simtrans(it+1)%bsim(in)=bB(polind)
                                simtrans(it+1)%lsim(in)=GridL(lB(polind))
                                simtrans(it+1)%lsimI(in)=lB(polind)
                                simtrans(it)%hconssim(in)=simtrans(it+1)%hsim(in)
                                simtrans(it)%csim(in)=consB(polind)

                            case(2) !Renting optimal

                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,bRpol(:,giExo,:,gij),bjps,gPh,simtrans(it+1)%bsim(in))
                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpol(:,giExo,:,gij),bjps,gPh,simtrans(it)%csim(in))
                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,hconsRpol(:,giExo,:,gij),bjps,gPh,simtrans(it+1)%hconssim(in))

                                simtrans(it)%buysim(in)=0
                                simtrans(it+1)%ownsim(in)=0
                                simtrans(it+1)%hsim(in)=0.0d0
                                simtrans(it+1)%msim(in)=0.0d0
                                simtrans(it+1)%lsim(in)=0.0d0
                                simtrans(it)%refisim(in)=0
                                simtrans(it)%foresim(in)=0
                                simtrans(it)%sellsim(in)=0

                            case(3) !Foreclosure optimal
                                simtrans(it)%refisim(in)=0

                                simtrans(it)%buysim(in) = 0
                                simtrans(it+1)%ownsim(in) = 0
                                simtrans(it)%foresim(in) = 1
                                simtrans(it)%hconssim(in) = simtrans(it)%hsim(in)

                                call BiLinInterp1(nb,GridB,ngpPh,PhGrid,bFpol(:,ljI,nm,hjI,giExo,:,gij),bj,gPh,simtrans(it+1)%bsim(in))
                                simtrans(it)%csim(in) = bj+y-FnTax(y)-FnQb(simtrans(it+1)%bsim(in),0.0d0)*simtrans(it+1)%bsim(in)
                                simtrans(it+1)%hsim(in) = 0.0d0
                                simtrans(it+1)%hsimI(in) = 1
                                simtrans(it+1)%msim(in) = 0.0d0
                                simtrans(it+1)%lsim(in) = 0.0d0
                                simtrans(it+1)%lsimI(in)=1

                            case(4) !Payment optimal
                                simtrans(it)%sellsim(in)=0
                                simtrans(it)%foresim(in)=0
                                simtrans(it)%buysim(in)=0
                                simtrans(it+1)%ownsim(in)=1
                                simtrans(it)%refisim(in)=0

                                simtrans(it+1)%hsim(in)=hj
                                simtrans(it+1)%hsimI(in)=hjI
                                simtrans(it)%hconssim(in)=hj

                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bPpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%bsim(in))
                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,mPpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%msim(in))
                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,lPpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%lsim(in))
                                simtrans(it+1)%msim(in)=simtrans(it+1)%msim(in)/hj
                                call basefun(GridL,nl,simtrans(it+1)%lsim(in),vals,inds)
                                simtrans(it+1)%lsimI(in)=inds(1)

                            case(5) !Refinance optimal

                                simtrans(it)%sellsim(in)=0
                                simtrans(it)%foresim(in)=0
                                simtrans(it)%buysim(in)=0
                                simtrans(it+1)%ownsim(in)=1
                                simtrans(it)%refisim(in)=1

                                simtrans(it+1)%hsim(in)=hj
                                simtrans(it+1)%hsimI(in)=hjI
                                simtrans(it)%hconssim(in)=hj

                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bNpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%bsim(in))
                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,mNpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%msim(in))
                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,lNpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%lsim(in))
                                call basefun(GridL,nl,simtrans(it+1)%lsim(in),vals,inds)
                                simtrans(it+1)%lsimI(in)=inds(1)

                         end select

                    endif

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                        !Now check the MPC
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                    !Compared to those who started off as renters, additional income from tax deduction
                    bjps=bjp+(Fntax(y)-FnTaxM(y,hj*(simtrans(it)%msim(in)+simtrans(it)%lsim(in)),min(0.0d0,bj)))+mpcval
                    bjmpc=bj+mpcval

                    if(bjps .lt. GridBS(1)) then
                        Wrentsim=Vmin
                        Wbuysim=Vmin
                    elseif(bjps .lt. 0.0d0) then
                        !Renting value

                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WrentS(:,giExo,:,gij),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WbuyS(:,giExo,:,gij),bjps,gPh,Wbuysim)
                    else
                        !Renting value

                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wrent(:,giExo,:,gij),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wbuy(:,giExo,:,gij),bjps,gPh,Wbuysim)
                    endif

                    helpWsell=max(Wbuysim,Wrentsim)

                    ! Figure out foreclosure value
                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wfore(:,ljI,nm,hjI,giExo,:,gij),bjmpc,gPh,helpWfore)
                    !First check if the moving shock hits

                    if(simtrans(it)%thetasim(in) .eq. 1) then !Hit by moving shock

                        !Only options are sell or foreclosure

                        if(helpWsell .ge. helpWfore) then !If selling is better

                            if(Wrentsim .gt. Wbuysim) then !If renting is better

                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpol(:,giExo,:,gij),bjps,gPh,simtrans(it)%csim2(in))
                            else

                                EVown(:,:,:,:,gij)=Phvals(1)*EEVown(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVown(:,:,:,:,gij,giExo,Phinds(2))

                                do hcc=1,nh
                                    fc=0.0d0
                                    do mcc=1,nm
                                    do lcc=1,nl
                                        index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                        mB(index)=mcc
                                        hB(index)=hcc
                                        lB(index)=lcc
                                        gih=hcc
                                        gil=lcc
                                        gim=mcc
                                        ghouse=GridH(gih)
                                        gmort=GridM(gim)
                                        gloc=GridL(gil)
                                        gliq=bjps+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans
                                        bccmax=(gliq+gmort*ghouse)*(1.0d0+rf)
                                        !Changed 10/13/15 to prevent new home buyers from using HELOCs for purchase
                                        if(giC .eq. ngpC) then
                                           bccmin=0.0d0
                                        elseif(HELOCVal .eq. 1) then
                                           bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                        else
                                           bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                        endif

                                        bccmid=(bccmax+bccmin)/2.0d0

                                        if(bccmax .gt. bccmin) then

                                            ftemp1=FnVFB(0.0d0)
                                            ftemp2=golden(bccmin,bccmid,bccmax,FnVFB,vftol,btemp)
                                            if(ftemp2<ftemp1) then
                                                ftemp1=ftemp2
                                            else
                                                btemp=0.0d0
                                            endif
                                        else
                                           helpWB(index)=Wmin
                                           btemp=0.0d0
                                           consB(index)=cmin
                                        endif
                                        bB(index)=btemp

                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qm(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,ql(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)
                                        consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                        helpWB(index)=-ftemp1
                                    enddo
                                    fc=FCGrid(giAgg) !Set fixed cost of mortgage for when mcc>1
                                    enddo
                                enddo
                                polind=maxloc(helpWB(1:nm*nh*nl),dim=1)

                                simtrans(it)%csim2(in)=consB(polind)
                             endif

                        else  !If foreclosing is better

                            call BiLinInterp1(nb,GridB,ngpPh,PhGrid,bFpol(:,ljI,nm,hjI,giExo,:,gij),bjmpc,gPh,btemp)
                            simtrans(it)%csim2(in) = max(cmin,bjmpc+y-FnTax(y)-FnQb(btemp,0.0d0)*btemp)
                        endif

                    else !Not hit by moving shock

                        if(mj .eq. 0.0d0) then
                            mc=1
                            call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wpay(:,ljI,mc,hjI,giExo,:,gij),bjmpc,gPh,Wpaysim)
                            call BiLinInterp1(nb,GridB,ngpPh,PhGrid,Wrefin(:,ljI,mc,hjI,giExo,:,gij),bjmpc,gPh,Wrefisim)

                        else

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,Wpay(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,Wpaysim)
                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,Wrefin(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,Wrefisim)

                        endif

                        helpVals(1) =   Wbuysim
                        helpVals(2) =   Wrentsim
                        helpVals(3) =   helpWfore
                        helpVals(4) =   Wpaysim
                        helpVals(5) =   Wrefisim

                        polind=maxloc(helpVals,dim=1)
                        !print*,helpVals
                        !pause
                        select case(polind)

                            case(1) !Buying optimal
                                EVown(:,:,:,:,gij)=Phvals(1)*EEVown(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVown(:,:,:,:,gij,giExo,Phinds(2))
                                do hcc=1,nh
                                    fc=0.0d0
                                    do mcc=1,nm
                                    do lcc=1,nl
                                        index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                        mB(index)=mcc
                                        hB(index)=hcc
                                        lB(index)=lcc
                                        gih=hcc
                                        gil=lcc
                                        gim=mcc
                                        ghouse=GridH(gih)
                                        gmort=GridM(gim)
                                        gloc=GridL(gil)
                                        gliq=bjps+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans
                                        bccmax=(gliq+gmort*ghouse)*(1.0d0+rf)
                                        !Changed 10/13/15 to prevent new home buyers from using HELOCs for purchase
                                        if(giC .eq. ngpC) then
                                           bccmin=0.0d0
                                        elseif(HELOCVal .eq. 1) then
                                           bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                        else
                                           bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                        endif

                                        bccmid=(bccmax+bccmin)/2.0d0

                                        if(bccmax .gt. bccmin) then

                                            ftemp1=FnVFB(0.0d0)
                                            ftemp2=golden(bccmin,bccmid,bccmax,FnVFB,vftol,btemp)
                                            if(ftemp2<ftemp1) then
                                                ftemp1=ftemp2
                                            else
                                                btemp=0.0d0
                                            endif
                                        else
                                           helpWB(index)=Wmin
                                           btemp=0.0d0
                                           consB(index)=cmin
                                        endif
                                        bB(index)=btemp
                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qm(:,gil,gim,gih,giExo,:,gij+1),btemp,Ph,qmtemp)
                                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,ql(:,gil,gim,gih,giExo,:,gij+1),btemp,Ph,qltemp)

                                        consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                        helpWB(index)=-ftemp1
                                    enddo
                                    fc=FCGrid(giAgg) !Set fixed cost of mortgage for when mcc>1
                                    enddo
                                enddo
                                polind=maxloc(helpWB(1:nh*nl*nm),dim=1)

                                simtrans(it)%csim2(in)=consB(polind)

                            case(2) !Renting optimal

                                CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpol(:,giExo,:,gij),bjps,gPh,simtrans(it)%csim2(in))

                            case(3) !Foreclosure optimal

                                call BiLinInterp1(nb,GridB,ngpPh,PhGrid,bFpol(:,ljI,nm,hjI,giExo,:,gij),bjmpc,gPh,btemp)
                                simtrans(it)%csim2(in) = bjmpc+y-FnTax(y)-FnQb(btemp,0.0d0)*btemp

                            case(4) !Payment optimal

                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpol(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                            case(5) !Refinance optimal

                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpol(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                        end select
                    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! END MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                endif
                simtrans(it+1)%bsim(in)=max(bmin,simtrans(it+1)%bsim(in))
                simtrans(it+1)%msim(in)=max(0.0d0,simtrans(it+1)%msim(in))
                simtrans(it+1)%lsim(in)=max(0.0d0,simtrans(it+1)%lsim(in))
                simtrans(it)%csim(in)=max(cmin,simtrans(it)%csim(in))
                simtrans(it)%csim2(in)=max(cmin,simtrans(it)%csim2(in))
                simtrans(it)%beqsim(in)=0.0d0
                simtrans(it)%mpcsim(in)=(simtrans(it)%csim2(in)-simtrans(it)%csim(in))/mpcval

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! RETIREMENT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            case(Jwork+1:Jtot-1)

                gij=ijret
                y=simtrans(it)%psim(in)
                giR=ip
                iR=giR

                giExo=(giAgg-1)*ngpR+giR

               !Check if comes into period with house
                if(simtrans(it)%ownsim(in) .eq. 0) then

                    !Check to see if stay renter or buy house
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,WrentRet(:,giExo,:,gij),bj,gPh,RentVal)
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,WbuyRet(:,giExo,:,gij),bj,gPh,BuyVal)

                    if(RentVal .ge. BuyVal ) then

                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,bRpolRet(:,giExo,:,gij),bj,gPh,simtrans(it+1)%bsim(in))
                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpolRet(:,giExo,:,gij),bj,gPh,simtrans(it)%csim(in))
                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,hconsRpolRet(:,giExo,:,gij),bj,gPh,simtrans(it+1)%hconssim(in))

                        simtrans(it)%buysim(in)=0
                        simtrans(it)%refisim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%sellsim(in)=0
                        simtrans(it+1)%ownsim(in)=0
                        simtrans(it+1)%hsim(in)=0.0d0
                        simtrans(it+1)%msim(in)=0.0d0
                        simtrans(it+1)%lsim(in)=0.0d0
                        simtrans(it+1)%lsimI(in)=1

                    else !Buyer
                        EVownRet(:,:,:,:,gij)=Phvals(1)*EEVownRet(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVownRet(:,:,:,:,gij,giExo,Phinds(2))

                        do hcc=1,nh
                            fc=0.0d0
                            do mcc=1,nm
                                do lcc=1,nl

                                    index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                    mB(index)=mcc
                                    hB(index)=hcc
                                    lB(index)=lcc
                                    gih=hcc
                                    gil=lcc
                                    gim=mcc
                                    ghouse=GridH(gih)
                                    gmort=GridM(gim)
                                    gloc=GridL(gil)
                                    gliq=bj+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))
                                    ftemp1=FnVFBRet(bmin)
                                    ftemp2=FnVFBRet(GridB(2))
                                    bccmax=bj+y-FnTax(y)-fc
                                    if(ftemp2<ftemp1) then
                                        ftemp1=golden(bmin,GridB(2),bccmax,FnVFBRet,vftol,btemp)
                                    else
                                        btemp=bmin
                                    endif
                                    bB(index)=btemp
                                    call basefun(GridB,nb,btemp,vals,inds)

                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qmret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qlret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)

                                    consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                    helpWB(index)=-ftemp1

                                enddo
                                fc=FCGrid(giAgg)
                            enddo
                        enddo

                        !maximize over the portfolio choice
                        polind=maxloc(helpWB,dim=1)

                        simtrans(it)%buysim(in)=1
                        simtrans(it+1)%bsim(in)=bB(polind)
                        simtrans(it)%hconssim(in)=GridH(hB(polind))
                        simtrans(it)%csim(in)=consB(polind)
                        simtrans(it)%refisim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%sellsim(in)=0
                        simtrans(it+1)%ownsim(in)=1
                        simtrans(it+1)%hsim(in)=GridH(hB(polind))
                        simtrans(it+1)%msim(in)=GridM(mB(polind))
                        simtrans(it+1)%lsim(in)=GridL(lB(polind))
                        simtrans(it+1)%lsimI(in)=lB(polind)
                        simtrans(it+1)%hsimI(in)=hB(polind)

                    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Now check the MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                    !Check to see if stay renter or buy house
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,WrentRet(:,giExo,:,gij),bjmpc,gPh,RentVal)
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,WbuyRet(:,giExo,:,gij),bjmpc,gPh,BuyVal)
                    if(RentVal .ge. BuyVal ) then

                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpolRet(:,giExo,:,gij),bjmpc,gPh,simtrans(it)%csim2(in))

                    else !Buyer
                        EVownRet(:,:,:,:,gij)=Phvals(1)*EEVownRet(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVownRet(:,:,:,:,gij,giExo,Phinds(2))

                        do hcc=1,nh
                            fc=0.0d0
                            do mcc=1,nm
                                do lcc=1,nl

                                    index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                    mB(index)=mcc
                                    hB(index)=hcc
                                    lB(index)=lcc
                                    gih=hcc
                                    gil=lcc
                                    gim=mcc
                                    ghouse=GridH(gih)
                                    gmort=GridM(gim)
                                    gloc=GridL(gil)
                                    gliq=bjmpc+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans

                                    bccmax=gliq+gmort*ghouse
                                    !Changed 10/13/15 to prevent new home buyers from using HELOCs for downpayment
                                    if(giC .eq. ngpC) then
                                       bccmin=0.0d0
                                    elseif(HELOCVal .eq. 1) then
                                       bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                    else
                                       bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                    endif

                                    bccmid=(bccmax+bccmin)/2.0d0

                                    if(bccmax .gt. bccmin) then
                                        ftemp1=FnVFBRet(0.0d0)
                                        ftemp2=golden(bccmin,bccmid,bccmax,FnVFBRet,vftol,btemp)
                                        if(ftemp2<ftemp1) then
                                            ftemp1=ftemp2
                                        else
                                            btemp=0.0d0
                                        endif
                                        helpWB(index)=-ftemp1
                                    else
                                        helpWB(index)=Wmin
                                        btemp=0.0d0
                                    endif
                                    bB(index)=btemp
                                    call basefun(GridB,nb,btemp,vals,inds)

                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qmret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qlret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)

                                    consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                    helpWB(index)=-ftemp1

                                enddo
                                fc=FCGrid(giAgg)
                            enddo
                        enddo

                        !maximize over the portfolio choice
                        polind=maxloc(helpWB,dim=1)

                        simtrans(it)%csim2(in)=consB(polind)

                    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!END MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                else  ! Homeowner

                    if(Modify .eq. 1) then

                        if(it .eq. 166) simtrans(it)%msim(in)=min(simtrans(it)%msim(in),0.95*PhPath(106))
                        if(it .eq. 206) simtrans(it)%msim(in)=min(simtrans(it)%msim(in),0.90*PhPath(106))

                    endif

                    hj=simtrans(it)%hsim(in)
                    lj=simtrans(it)%lsim(in)*(1.0d0+rl)
                    ljI=simtrans(it)%lsimI(in)
                    mj=simtrans(it)%msim(in)*(1.0d0+rm)
                    hjI=simtrans(it)%hsimI(in)

                    bjp=simtrans(it)%bsim(in)+hj*(Ph-mj-lj)-FnAdj(Ph*hj)

                    call basefun(GridB,nb,bj,vals,inds)

                    if(mj .eq. 0.0d0) then
                        mc=1
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WpayRet(:,ljI,mc,hjI,giExo,:,gij),bj,gPh,Wpaysim)
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WrefinRet(:,ljI,mc,hjI,giExo,:,gij),bj,gPh,Wrefisim)

                    else

                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WpayRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,Wpaysim)
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WrefinRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,Wrefisim)

                    endif

                    !Compared to those who started off as renters, additional income from tax deduction
                    bjps=bjp+(Fntax(y)-FnTaxM(y,hj*(simtrans(it)%msim(in)+simtrans(it)%lsim(in)),min(0.0d0,bj)))
                    if(bjps .lt. GridBS(1)) then
                        Wrentsim=Vmin
                        Wbuysim=Vmin
                    elseif(bjps .lt. 0.0d0) then
                        !Renting value

                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WrentRetS(:,giExo,:,gij+1),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WbuyRetS(:,giExo,:,gij+1),bjps,gPh,Wbuysim)
                    else
                        !Renting value

                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WrentRet(:,giExo,:,gij+1),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WbuyRet(:,giExo,:,gij+1),bjps,gPh,Wbuysim)
                    endif

                    helpWsell=max(Wbuysim,Wrentsim)

                    ! Figure out foreclosure value
                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WforeRet(:,ljI,nm,hjI,giExo,:,gij),bj,gPh,helpWfore)
                    !First check if the moving shock hits

                    helpVals(1) =   Wbuysim
                    helpVals(2) =   Wrentsim
                    helpVals(3) =   helpWfore
                    helpVals(4) =   Wpaysim
                    helpVals(5) =   Wrefisim

                    polind=maxloc(helpVals,dim=1)

                    select case(polind)

                        case(1) !Buying optimal

                            EVownRet(:,:,:,:,gij)=Phvals(1)*EEVownRet(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVownRet(:,:,:,:,gij,giExo,Phinds(2))

                            do hcc=1,nh
                                fc=0.0d0
                                do mcc=1,nm
                                do lcc=1,nl
                                    index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                    mB(index)=mcc
                                    hB(index)=hcc
                                    lB(index)=lcc
                                    gih=hcc
                                    gil=lcc
                                    gim=mcc
                                    ghouse=GridH(gih)
                                    gmort=GridM(gim)
                                    gloc=GridL(gil)
                                    gliq=bjps+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))
                                    ftemp1=FnVFBRet(bmin)
                                    ftemp2=FnVFBRet(GridB(2))
                                    bccmax=gliq-fc
                                    if(ftemp2<ftemp1) then
                                        ftemp1=golden(bmin,GridB(2),bccmax,FnVFBRet,vftol,btemp)
                                    else
                                        btemp=bmin
                                    endif
                                    bB(index)=btemp

                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qmret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qlret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)

                                    consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                    helpWB(index)=-ftemp1

                                enddo
                                fc=mortfc !Set fixed cost of mortgage for when mcc>1
                                enddo
                            enddo
                            polind=maxloc(helpWB(1:nh*nl*nm),dim=1)

                            if(hB(polind) .eq. hjI .and. ( (helpWB(polind) .lt. Wpaysim) .or. (helpWB(polind) .lt. Wrefisim))) then

                                if(Wpaysim .ge. Wrefisim) then
                                    simtrans(it)%sellsim(in)=0
                                    simtrans(it)%foresim(in)=0
                                    simtrans(it)%buysim(in)=0
                                    simtrans(it)%refisim(in)=0
                                    simtrans(it+1)%ownsim(in)=1

                                    simtrans(it+1)%hsim(in)=hj
                                    simtrans(it+1)%hsimI(in)=hjI
                                    simtrans(it)%hconssim(in)=hj

                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%bsim(in))
                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,mPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%msim(in))
                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,lPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%lsim(in))

                                    call basefun(GridL,nl,simtrans(it+1)%lsim(in),vals,inds)
                                    simtrans(it+1)%lsimI(in)=inds(1)
                                    simtrans(it+1)%msim(in)=simtrans(it+1)%msim(in)/hj
                                else
                                    simtrans(it)%sellsim(in)=0
                                    simtrans(it)%foresim(in)=0
                                    simtrans(it)%buysim(in)=0
                                    simtrans(it+1)%ownsim(in)=1

                                    simtrans(it)%refisim(in)=1

                                    simtrans(it+1)%hsim(in)=hj
                                    simtrans(it+1)%hsimI(in)=hjI
                                    simtrans(it)%hconssim(in)=hj

                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%bsim(in))
                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,mNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%msim(in))
                                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,lNpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%lsim(in))

                                    call basefun(GridL,nl,simtrans(it+1)%lsim(in),vals,inds)
                                    simtrans(it+1)%lsimI(in)=inds(1)
                                    simtrans(it+1)%msim(in)=simtrans(it+1)%msim(in)

                                endif

                            else
                                simtrans(it)%sellsim(in)=1
                                simtrans(it)%foresim(in)=0
                                simtrans(it)%buysim(in)=1
                                simtrans(it+1)%ownsim(in)=1
                                simtrans(it+1)%hsim(in)=GridH(hB(polind))
                                simtrans(it+1)%hsimI(in)=hB(polind)
                                simtrans(it+1)%msim(in)=GridM(mB(polind))
                                simtrans(it+1)%bsim(in)=bB(polind)
                                simtrans(it+1)%lsim(in)=GridL(lB(polind))
                                simtrans(it+1)%lsimI(in)=lB(polind)
                                simtrans(it)%hconssim(in)=simtrans(it+1)%hsim(in)
                                simtrans(it)%csim(in)=consB(polind)
                            endif
                    case(2) !Renting optimal

                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,bRpolRet(:,giExo,:,gij),bjps,gPh,simtrans(it+1)%bsim(in))
                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpolRet(:,giExo,:,gij),bjps,gPh,simtrans(it)%csim(in))
                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,hconsRpolRet(:,giExo,:,gij),bjps,gPh,simtrans(it+1)%hconssim(in))

                        simtrans(it)%buysim(in)=0
                        simtrans(it)%refisim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%sellsim(in)=0
                        simtrans(it+1)%ownsim(in)=0
                        simtrans(it+1)%hsim(in)=0.0d0
                        simtrans(it+1)%msim(in)=0.0d0
                        simtrans(it+1)%lsim(in)=0.0d0
                        simtrans(it+1)%lsimI(in)=1

                    case(3) !Foreclosure optimal

                        simtrans(it)%buysim(in) = 0
                        simtrans(it+1)%ownsim(in) = 0
                        simtrans(it)%foresim(in) = 1
                        simtrans(it)%hconssim(in) = simtrans(it)%hsim(in)
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,bFpolRet(:,ljI,nm,hjI,giExo,:,ijret),bj,gPh,simtrans(it+1)%bsim(in))
                        simtrans(it)%csim(in) = bj+y-FnTax(y)-FnQb(simtrans(it+1)%bsim(in),0.0d0)*simtrans(it+1)%bsim(in)
                        simtrans(it+1)%hsim(in) = 0.0d0
                        simtrans(it+1)%hsimI(in) = 1
                        simtrans(it+1)%msim(in) = 0.0d0
                        simtrans(it+1)%lsim(in) = 0.0d0
                        simtrans(it+1)%lsimI(in)=1

                    case(4) !Payment optimal
                        simtrans(it)%sellsim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%buysim(in)=0
                        simtrans(it)%refisim(in)=0
                        simtrans(it+1)%ownsim(in)=1

                        simtrans(it+1)%hsim(in)=hj
                        simtrans(it+1)%hsimI(in)=hjI
                        simtrans(it)%hconssim(in)=hj

                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%bsim(in))
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,mPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%msim(in))
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,lPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%lsim(in))

                        call basefun(GridL,nl,simtrans(it+1)%lsim(in),vals,inds)
                        simtrans(it+1)%lsimI(in)=inds(1)
                        simtrans(it+1)%msim(in)=simtrans(it+1)%msim(in)/hj

                    case(5) !Refinance optimal

                        simtrans(it)%sellsim(in)=0
                        simtrans(it)%foresim(in)=0
                        simtrans(it)%buysim(in)=0
                        simtrans(it+1)%ownsim(in)=1

                        simtrans(it)%refisim(in)=1

                        simtrans(it+1)%hsim(in)=hj
                        simtrans(it+1)%hsimI(in)=hjI
                        simtrans(it)%hconssim(in)=hj

                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%bsim(in))
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,mNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%msim(in))
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,lNpol(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it+1)%lsim(in))

                        call basefun(GridL,nl,simtrans(it+1)%lsim(in),vals,inds)
                        simtrans(it+1)%lsimI(in)=inds(1)
                        simtrans(it+1)%msim(in)=simtrans(it+1)%msim(in)

                     end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Now check the MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                    bjpmpc=simtrans(it)%bsim(in)+hj*(Ph-mj-lj)-FnAdj(Ph*hj)+mpcval

                    call basefun(GridB,nb,bjmpc,vals,inds)

                    if(mj .eq. 0.0d0) then
                        mc=1
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WpayRet(:,ljI,mc,hjI,giExo,:,gij),bjmpc,gPh,Wpaysim)
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WrefinRet(:,ljI,mc,hjI,giExo,:,gij),bjmpc,gPh,Wrefisim)

                    else

                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WpayRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,Wpaysim)
                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WrefinRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,Wrefisim)

                    endif

                    !Compared to those who started off as renters, additional income from tax deduction
                    bjps=bjpmpc+(Fntax(y)-FnTaxM(y,hj*(simtrans(it)%msim(in)+simtrans(it)%lsim(in)),min(0.0d0,bj)))
                    if(bjps .lt. GridBS(1)) then
                        Wrentsim=Vmin
                        Wbuysim=Vmin
                    elseif(bjps .lt. 0.0d0) then
                        !Renting value

                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WrentRetS(:,giExo,:,gij+1),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nbs,GridBS,ngpPh,PhGrid,WbuyRetS(:,giExo,:,gij+1),bjps,gPh,Wbuysim)
                    else
                        !Renting value

                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WrentRet(:,giExo,:,gij+1),bjps,gPh,Wrentsim)

                        !Buying value
                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WbuyRet(:,giExo,:,gij+1),bjps,gPh,Wbuysim)
                    endif

                    helpWsell=max(Wbuysim,Wrentsim)

                    ! Figure out foreclosure value
                    call BiLinInterp1(nb,GridB,ngpPh,PhGrid,WforeRet(:,ljI,nm,hjI,giExo,:,gij),bjmpc,gPh,helpWfore)
                    !First check if the moving shock hits

                    helpVals(1) =   Wbuysim
                    helpVals(2) =   Wrentsim
                    helpVals(3) =   helpWfore
                    helpVals(4) =   Wpaysim
                    helpVals(5) =   Wrefisim

                    polind=maxloc(helpVals,dim=1)

                    select case(polind)

                    case(1) !Buying optimal

                        EVownRet(:,:,:,:,gij)=Phvals(1)*EEVownRet(:,:,:,:,gij,giExo,Phinds(1))+Phvals(2)*EEVownRet(:,:,:,:,gij,giExo,Phinds(2))

                        do hcc=1,nh
                            fc=0.0d0
                            do mcc=1,nm
                            do lcc=1,nl
                                index=(hcc-1)*nm*nl+(mcc-1)*nl+lcc
                                mB(index)=mcc
                                hB(index)=hcc
                                lB(index)=lcc
                                gih=hcc
                                gil=lcc
                                gim=mcc
                                ghouse=GridH(gih)
                                gmort=GridM(gim)
                                gloc=GridL(gil)
                                gliq=bjps+y-FnTax(y)-Ph*GridH(hcc)-fc-FnMaint(GridH(hcc))+htrans
                                bccmax=gliq+gmort*ghouse
                                !Changed 10/13/15 to prevent new home buyers from using HELOCs for downpayment
                                !bccmin=0.0d0
    !                            bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                if(giC .eq. ngpC) then
                                   bccmin=0.0d0
                                elseif(HELOCVal .eq. 1) then
                                   bccmin=-HELOCGrid(giAgg)*gPh*ghouse
                                else
                                   bccmin=min(-HELOCGrid(giAgg)*((gPh-gmort)*ghouse),0.0d0)
                                endif

                                bccmid=(bccmax+bccmin)/2.0d0

                                if(bccmax .gt. bccmin) then
                                    ftemp1=FnVFBRet(0.0d0)
                                    ftemp2=golden(bccmin,bccmid,bccmax,FnVFBRet,vftol,btemp)
                                    if(ftemp2<ftemp1) then
                                        ftemp1=ftemp2
                                    else
                                        btemp=0.0d0
                                    endif
                                    helpWB(index)=-ftemp1
                                else
                                    helpWB(index)=Wmin
                                    btemp=0.0d0
                                endif
                                bB(index)=btemp

                                call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qmret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qmtemp)
                                call BiLinInterp1(nb,GridB,ngpPh,PhGrid,qlret(:,gil,gim,gih,giExo,:,gij+1),btemp,gPh,qltemp)

                                consB(index)=gliq+(qltemp*gloc+qmtemp*gmort)*ghouse-FnQb(btemp,0.0d0)*btemp
                                helpWB(index)=-ftemp1

                            enddo
                            fc=FCGrid(giAgg) !Set fixed cost of mortgage for when mcc>1
                            enddo
                        enddo
                        polind=maxloc(helpWB(1:nh*nl*nm),dim=1)

                        if(hB(polind) .eq. hjI .and. ( (helpWB(polind) .lt. Wpaysim) .or. (helpWB(polind) .lt. Wrefisim))) then

                            if(Wpaysim .ge. Wrefisim) then

                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                            else

                                call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                            endif

                        else

                           simtrans(it)%csim2(in)=consB(polind)

                        endif
                    case(2) !Renting optimal

                        CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpolRet(:,giExo,:,gij),bjps,gPh,simtrans(it)%csim2(in))

                    case(3) !Foreclosure optimal

                        call BiLinInterp1(nb,GridB,ngpPh,PhGrid,bFpolRet(:,ljI,nm,hjI,giExo,:,ijret),bjmpc,gPh,btemp)
                        simtrans(it)%csim2(in) = bjmpc+y-FnTax(y)-FnQb(btemp,0.0d0)*btemp

                    case(4) !Payment optimal

                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                    case(5) !Refinance optimal

                        call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                     end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!END MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                  endif

                simtrans(it)%beqsim(in)=0.0d0
                simtrans(it+1)%bsim(in)=max(bmin,simtrans(it+1)%bsim(in))
                simtrans(it+1)%msim(in)=max(0.0d0,simtrans(it+1)%msim(in))
                simtrans(it+1)%lsim(in)=max(0.0d0,simtrans(it+1)%lsim(in))
                simtrans(it)%csim(in)=max(cmin,simtrans(it)%csim(in))
                simtrans(it)%csim2(in)=max(cmin,simtrans(it)%csim2(in))
                simtrans(it)%mpcsim(in)=(simtrans(it)%csim2(in)-simtrans(it)%csim(in))/mpcval
                simtrans(it)%beqsim(in)=0.0d0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! LAST PERIOD OF LIFE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            case(Jtot)

                iR=ip
                giR=iR
                y=simtrans(it)%psim(in)

                gij=ijret

                giExo=(giAgg-1)*ngpR+giR

               !Check if comes into period with house
                if(simtrans(it)%ownsim(in) .eq. 0) then

                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,bRpolRet(:,giExo,:,gij),bj,gPh,simtrans(it)%beqsim(in))
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpolRet(:,giExo,:,gij),bj,gPh,simtrans(it)%csim(in))
                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,hconsRpolRet(:,giExo,:,gij),bj,gPh,simtrans(it+1)%hconssim(in))

                    simtrans(it)%buysim(in)=0
                    simtrans(it)%refisim(in)=0
                    simtrans(it)%foresim(in)=0
                    simtrans(it)%sellsim(in)=0

                    CALL BiLinInterp1(nb,GridB,ngpPh,PhGrid,consRpolRet(:,giExo,:,gij),bjmpc,gPh,simtrans(it)%csim2(in))

                else

                    if(Modify .eq. 1) then

                        if(it .eq. 166) simtrans(it)%msim(in)=min(simtrans(it)%msim(in),0.95*PhPath(106))
                        if(it .eq. 206) simtrans(it)%msim(in)=min(simtrans(it)%msim(in),0.90*PhPath(106))

                    endif

                    hj=simtrans(it)%hsim(in)
                    lj=simtrans(it)%lsim(in)*(1.0d0+rl)
                    ljI=simtrans(it)%lsimI(in)
                    mj=simtrans(it)%msim(in)*(1.0d0+rm)
                    hjI=simtrans(it)%hsimI(in)

                    bjp=bj+hj*(Ph-mj-lj)-FnAdj(Ph*hj)

                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    !Start off with the sellers problem
                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    !Compute value of liquid assets after the sale

                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WforeRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,helpWfore)
                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WforeRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,WpaySim)
                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WforeRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,helpWsell)

                    HelpVals(1)=helpWsell
                    HelpVals(2)=Wpaysim
                    HelpVals(3)=helpWfore

                    polind=maxloc(HelpVals(1:3),dim=1)

                    simtrans(it)%buysim(in)=0
                    simtrans(it)%refisim(in)=0

                    select case(polind)

                        case(1)   ! Sell

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,hconsNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%hconssim(in))
                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bNpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%beqsim(in))

                            simtrans(it)%sellsim(in)=1
                            simtrans(it)%foresim(in)=0
                            simtrans(it)%buysim(in)=0
                            simtrans(it)%refisim(in)=0

                        case(2)   ! Pay
                            simtrans(it)%buysim(in)=0
                            simtrans(it)%refisim(in)=0
                            simtrans(it)%sellsim(in)=0

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bPpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%beqsim(in))

                            simtrans(it)%beqsim(in)=simtrans(it)%beqsim(in)+gPh*hj
                            simtrans(it)%hconssim(in)=hj
                            simtrans(it)%foresim(in)=0

                        case(3)   ! Foreclose
                            simtrans(it)%sellsim(in)=0

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consFpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%csim(in))
                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,bFpolRet(:,ljI,:,hjI,giExo,:,gij),bj,simtrans(it)%msim(in),gPh,simtrans(it)%beqsim(in))

                            simtrans(it)%hconssim(in)=hj*0.5d0/(dble(BiAnnual)+1.0d0)+(1.0d0-0.5d0/(dble(BiAnnual)+1.0d0))*GridHR(1)
                            simtrans(it)%foresim(in)=1
                            simtrans(it)%buysim(in)=0
                            simtrans(it)%refisim(in)=0

                    end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Do MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WforeRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,helpWfore)
                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WforeRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,WpaySim)
                    call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,WforeRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,helpWsell)

                    HelpVals(1)=helpWsell
                    HelpVals(2)=Wpaysim
                    HelpVals(3)=helpWfore

                    polind=maxloc(HelpVals(1:3),dim=1)

                    select case(polind)

                        case(1)   ! Sell

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consNpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                        case(2)   ! Pay

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consPpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                        case(3)   ! Foreclose

                            call TriLinInterp1(nb,GridB,nm,GridM,ngpPh,PhGrid,consFpolRet(:,ljI,:,hjI,giExo,:,gij),bjmpc,simtrans(it)%msim(in),gPh,simtrans(it)%csim2(in))

                    end select

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!END MPC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            endif

            simtrans(it)%csim(in)=max(cmin,simtrans(it)%csim(in))
            simtrans(it)%csim2(in)=max(cmin,simtrans(it)%csim2(in))
            simtrans(it)%mpcsim(in)=(simtrans(it)%csim2(in)-simtrans(it)%csim(in))/mpcval

            !Then you die

         end select
      end do
      !$OMP END PARALLEL DO

   end do

end subroutine SimulateBPP

max_stars_repo_name: lnsongxf/housing-boom-bust

max_issues_repo_name: lnsongxf/housing-boom-bust

max_forks_repo_name: kurtmitman/housing-boom-bust


text:
      SUBROUTINE RCOVLS (LASTSS)
C
C     THIS ROUTINE CREATES THE SOLN ITEM FOR A LOWER LEVEL SUBSTRUCTURE,
C     LASTSS, BY EDITING THAT OF THE SOLUTION SUBSTRUCTURE FSS.
C
      INTEGER         DRY        ,STEP       ,FSS        ,RFNO       ,
     1                UINMS      ,UA         ,LASTSS(2)  ,EQSS       ,
     2                RC         ,SOLN       ,SRD        ,SWRT       ,
     3                IZ(3)      ,EOI        ,EOG        ,NAME(2)
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /BLANK / DRY        ,LOOP       ,STEP       ,FSS(2)     ,
     1                RFNO       ,NEIGV      ,LUI        ,UINMS(2,5) ,
     2                NOSORT     ,UTHRES     ,PTHRES     ,QTHRES
      COMMON /RCOVCR/ ICORE      ,LCORE      ,BUF1       ,BUF2       ,
     1                BUF3       ,BUF4       ,SOF1       ,SOF2       ,
     2                SOF3
      COMMON /RCOVCM/ MRECVR     ,UA         ,PA         ,QA         ,
     1                IOPT       ,RSS(2)     ,ENERGY     ,UIMPRO     ,
     2                RANGE(2)   ,IREQ       ,LREQ       ,LBASIC
      COMMON /SYSTEM/ SYSBUF     ,NOUT
      COMMON /ZZZZZZ/ Z(1)
      EQUIVALENCE     (Z(1),IZ(1))
      DATA    NAME  / 4HRCOV, 4HLS  /
      DATA    EQSS  , SOLN  / 4HEQSS,4HSOLN /
      DATA    SRD   , SWRT  / 1,2   /
      DATA    EOG   , EOI   / 2,3   /
C
C     CREATE SOLN ITEM FOR THE RECOVERED SUBSTRUCTURE
C
      IF (RFNO .EQ. 3) GO TO 490
C
C     OBTAIN LIST OF CONTRIBUTING BASIC SUBSTRUCTURES FROM EQSS.
C     STORE IN OPEN CORE AT ICORE.
C
      CALL SFETCH (LASTSS,EQSS,SRD,RC)
      CALL SUREAD (Z(ICORE),4,NWDS,RC)
      NSS = IZ(ICORE+2)
      IF (LCORE .LE. ICORE+2*NSS-1) GO TO 9008
      CALL SUREAD (Z(ICORE),2*NSS,NWDS,RC)
C
C     CONSTRUCT SOLN GROUP 0 IN OPEN CORE AT IG0.  TWO SLOTS FOR THE
C     NUMBER OF LOADS ON EACH SUBSTRUCTURE.  FIRST IS FOR OLD FSS
C     SOLN, SECOND FOR NEW ONE.
C
      IG0 = ICORE + 2*NSS
      CALL SFETCH (FSS,SOLN,SRD,RC)
      IF (RC .EQ. 1) GO TO 462
      CALL SMSG (RC-2,SOLN,FSS)
      GO TO 498
  462 CALL SUREAD (Z(IG0),5,NWDS,RC)
      ISOL = IZ(IG0+2)
      IF (ISOL .NE. RFNO) GO TO 6369
      NS = IZ(IG0+3)
      NC = IZ(IG0+4)
      IF (IG0+4+4*NS .GT. LCORE) GO TO 9008
      DO 465 I = 1,NS
      CALL SUREAD (Z(IG0+1+4*I),3,NWDS,RC)
      IZ(IG0+4*I+4) = -65535
      DO 463 J = 1,NSS
      IF (IZ(IG0+4*I+1) .NE. IZ(ICORE+2*J-2)) GO TO 463
      IF (IZ(IG0+4*I+2) .NE. IZ(ICORE+2*J-1)) GO TO 463
      IZ(IG0+4*I+4) = IZ(IG0+4*I+3)
      GO TO 465
  463 CONTINUE
  465 CONTINUE
      IF (RFNO.EQ.8 .OR. RFNO.EQ.9) GO TO 600
      I = 1
      CALL SJUMP (I)
C
C     STATICS SOLUTION ITEM
C
C     READ ALL GROUPS OF THE OLD FSS SOLN INTO OPEN CORE AT IGS.
C     AS EACH ONE IS READ, ELIMINATE LOAD VECTORS WHICH DO NOT
C     APPLY TO THE NEW SOLN BY SETTING THEIR LOAD VECTOR
C     NUMBERS TO -65535.
C     UPDATE THE NUMBER OF LOAD VECTORS WHICH DO APPLY.
C
      IGS = IG0 + 4*NS + 5
      JGS = IGS
      DO 478 I = 1,NC
      CALL SUREAD (Z(JGS),1,NWDS,RC)
      N = IABS(IZ(JGS))
      IF (JGS+N*2 .GT. LCORE) GO TO 9008
      CALL SUREAD (Z(JGS+1),-1,NWDS,RC)
      NL = 0
      IF (N .EQ. 0) GO TO 477
      DO 475 J = 1,N
      LVN = IZ(JGS+2*J-1)
C
C     FIND SUBSTRUCTURE WHERE LVN IS APPLIED FOR FSS SOLN ITEM.
C
      L1 = 0
      L2 = 0
      DO 470 K = 1,NS
      IF (LVN .GT. L2+IZ(IG0+4*K+3)) GO TO 468
      IF (IZ(IG0+4*K+4) .LT. 0) GO TO 471
      LVN = LVN - L1
      NL  = NL + 1
      GO TO 472
  468 IF (IZ(IG0+4*K+4) .LT. 0) L1 = L1 + IZ(IG0+4*K+3)
      L2 = L2 + IZ(IG0+4*K+3)
  470 CONTINUE
  471 LVN = -65535
  472 IZ(JGS+2*J-1) = LVN
  475 CONTINUE
      IF (IZ(JGS) .LT. 0) NL = -NL
      IZ(JGS) = NL
  477 JGS = JGS + 2*N + 1
  478 CONTINUE
C
C     WRITE THE NEW SOLN FOR THE RECOVERED SUBSTRUCTURE ON THE SOF.
C     IN CASE USER FORGOT TO EDIT OUT THIS SOLN FROM A PREVIOUS
C     RUN, DELETE IT TO AVOID LOSING OR SCREWING UP THE RECOVERED
C     DISPLACEMENTS.
C
      CALL DELETE (LASTSS,SOLN,RC)
      IZ(IG0+3) = NSS
      RC = 3
      CALL SFETCH (LASTSS,SOLN,SWRT,RC)
      CALL SUWRT (Z(IG0),5,1)
      DO 480 I = 1,NS
      IF (IZ(IG0+4*I+4) .LT. 0) GO TO 480
      CALL SUWRT (Z(IG0+4*I+1),2,1)
      CALL SUWRT (Z(IG0+4*I+4),1,1)
  480 CONTINUE
      CALL SUWRT (0,0,EOG)
      JGS = IGS
      DO 488 I = 1,NC
      K  = 0
      NL = IZ(JGS)
      JGS= JGS + 1
      CALL SUWRT (NL,1,1)
      IF (NL .EQ. 0) GO TO 485
      NL = IABS(NL)
  482 IF (IZ(JGS) .EQ. -65535) GO TO 484
      CALL SUWRT (IZ(JGS),2,1)
      K = K + 1
  484 JGS = JGS + 2
      IF (K .LT. NL) GO TO 482
  485 IF (IZ(JGS) .NE. -65535) GO TO 486
      JGS = JGS + 2
      GO TO 485
  486 CALL SUWRT (0,0,EOG)
  488 CONTINUE
      CALL SUWRT (0,0,EOI)
      GO TO 498
C
C     MODAL SOLUTION ITEM
C
C     FOR MODAL COPY THE SOLN UNCHANGED.  IN CASE THE USER FORGOT
C     TO EDIT OUT THIS SOLN FROM A PREVIOUS RUN, DELETE IT TO AVOID
C     LOSING OR SCREWING UP THE RECOVERED DISPLACEMENTS.
C
  490 CALL DELETE (LASTSS,SOLN,RC)
      CALL SFETCH (FSS,SOLN,SRD,RC)
      CALL SUREAD (IZ(ICORE),-1,NWDS,RC)
      ISOL = IZ(ICORE+2)
      IF (ISOL .NE. RFNO) GO TO 6369
      IF (IZ(ICORE+3) .GT. 0) GO TO 492
      RC = 3
      CALL SFETCH (LASTSS,SOLN,SWRT,RC)
      CALL SUWRT (Z(ICORE),4,EOG)
      CALL SUWRT (0,0,EOI)
      GO TO 498
  492 IF (LCORE .LT. ICORE+7*IZ(ICORE+3)+3) GO TO 9008
      CALL SUREAD (IZ(ICORE+4),-1,NWDS,RC)
      RC = 3
      CALL SFETCH (LASTSS,SOLN,SWRT,RC)
      CALL SUWRT (Z(ICORE),4,EOG)
      CALL SUWRT (Z(ICORE+4),7*IZ(ICORE+3),EOG)
      CALL SUWRT (0,0,EOI)
      GO TO 498
C
C     DYNAMIC SOLUTION ITEM
C
C     READ IN STATIC LOAD SETS
C
  600 INCR = 1
      IF (RFNO .EQ. 8) INCR = 2
      IGS = IG0 + 4*NS + 5
      CALL SUREAD (Z(IGS),1,NWDS,RC)
      NSL = IZ(IGS)
      LSL = NSL*INCR
      NSLL= 0
      IF (NSL .EQ. 0) GO TO 660
      IF (IGS+NSL .GT. LCORE) GO TO 9008
      CALL SUREAD (Z(IGS+1),NSL,NWDS,RC)
C
C     FLAG THOSE STATIC LOAD IDS THAT ARE NOT IN THE LOWER LEVEL
C     SUBSTRUCTURE AND RENUMBER THOSE THAT ARE LEFT
C
      DO 650 J = 1,NSL
      LVN = IZ(IGS+J)
      L1  = 0
      L2  = 0
      DO 620 K = 1,NS
      IF (LVN .GT. L2+IZ(IG0+4*K+3)) GO TO 610
      IF (IZ(IG0+4*K+4) .LT. 0) GO TO 630
      LVN  = LVN - L1
      NSLL = NSLL + 1
      GO TO 640
  610 IF (IZ(IG0+4*K+4) .LT. 0) L1 = L1 + IZ(IG0+4*K+3)
      L2 = L2 + IZ(IG0+4*K+3)
  620 CONTINUE
  630 LVN = -65535
  640 IZ(IGS+J) = LVN
  650 CONTINUE
C
C     COPY THE FREQUENCY OR TIME STEP RECORD INTO CORE
C
  660 I = 1
      CALL SJUMP (I)
      ISTEP = IGS + NSL + 1
      IF (ISTEP+NC .GT. LCORE) GO TO 9008
      CALL SUREAD (IZ(ISTEP),-1,NWDS,RC)
C
C     COPY IN ALL LOAD FACTOR DATA
C
      IF (NSLL .EQ. 0) GO TO 675
      JGS = ISTEP + NC
      DO 670 I = 1,NC
      IF (JGS+LSL .GT. LCORE) GO TO 9008
      CALL SUREAD (Z(JGS),-1,NWDS,RC)
  670 JGS = JGS + LSL
C
C     WRITE THE NEW SOLN ITEM FOR THE RECOVERED SUBSTRUCTURE.  IN CASE
C     THE USER FORGOT TO EDIT OUT THIS SOLN FROM A PREVIOUS RUN,
C     DELETE IT TO AVOID LOSING OR SCREWING UP THE RECOVERED
C     DISPLACEMENTS
C
  675 CALL DELETE (LASTSS,SOLN,RC)
      RC = 3
      CALL SFETCH (LASTSS,SOLN,SWRT,RC)
      IZ(IG0+3) = NSS
      CALL SUWRT (Z(IG0),5,1)
      DO 680 I = 1,NS
      IF (IZ(IG0+4*I+4) .LT. 0) GO TO 680
      CALL SUWRT (Z(IG0+4*I+1),2,1)
      CALL SUWRT (Z(IG0+4*I+4),1,1)
  680 CONTINUE
      CALL SUWRT (NSLL,1,1)
      IF (NSLL .EQ. 0) GO TO 700
      DO 690 I = 1,NSL
      IF (Z(IGS+I) .LT. 0) GO TO 690
      CALL SUWRT (Z(IGS+I),1,1)
  690 CONTINUE
  700 CALL SUWRT (0,0,EOG)
C
C     COPY THE TIME OR FREQUENCY STEP INFO TO SOF.
C
      CALL SUWRT (Z(ISTEP),NC,EOG)
C
C     COPY LOAD FACTORS FOR EACH STEP TO SOF EDITING OUT THOSE
C     THAT NO LONGER PARTICIAPTE
C
      IF (NSLL .EQ. 0) GO TO 730
      KGS = ISTEP + NC
      DO 720 I = 1,NC
      K = 1
      DO 710 J = 1,NSL
      IF (Z(IGS+J) .LT. 0) GO TO 710
      CALL SUWRT (Z(KGS+K-1),INCR,1)
      K = K + INCR
  710 CONTINUE
      CALL SUWRT (0,0,EOG)
      KGS = KGS + LSL
  720 CONTINUE
C
  730 CALL SUWRT (0,0,EOI)
      GO TO 498
C
C     NORMAL RETURN
C
  498 RETURN
C
C     ERROR PROCESSING
C
 6369 WRITE (NOUT,63690) UFM,ISOL,RFNO
      GO TO 9100
 9008 N = 8
      CALL MESAGE (N,0,NAME)
 9100 IOPT = -1
      CALL SOFCLS
      RETURN
C
63690 FORMAT (A23,' 6369.  SOLN ITEM HAS INCORRECT RIGID FORMAT NUMBER',
     1       /31X,'SOLUTION RIGID FORMAT WAS',I5,
     2       ' AND CURRENT NASTRAN EXECUTION RIGID FORMAT IS',I5)
      END

ext: f


text:
!
! AtmProfile_Inspect
!
! Program to inspect the contents of an AtmProfile file.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 24-May-2010
!                       paul.vandelst@noaa.gov
!

PROGRAM AtmProfile_Inspect

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module usage
  USE File_Utility     , ONLY: File_Exists
  USE Message_Handler  , ONLY: SUCCESS, FAILURE, Program_Message, Display_Message
  USE AtmProfile_Define, ONLY: AtmProfile_type, &
                               AtmProfile_Destroy, &
                               Inspect => AtmProfile_Inspect
  USE AtmProfile_IO    , ONLY: AtmProfile_InquireFile, &
                               AtmProfile_ReadFile
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'AtmProfile_Inspect'
  CHARACTER(*), PARAMETER :: PROGRAM_RCS_ID = &

  ! ---------
  ! Variables
  ! ---------
  INTEGER :: err_stat
  CHARACTER(256) :: msg
  CHARACTER(256) :: Filename
  INTEGER :: n_Profiles, m
  TYPE(AtmProfile_type) :: A(1)

  ! Output program header
  CALL Program_Message( PROGRAM_NAME, &
                        'Program to display the contents of an '//&
                        'AtmProfile file to stdout.', &
                        '$Revision$' )

  ! Get the filename
  WRITE( *,FMT='(/5x,"Enter the AtmProfile filename: ")',ADVANCE='NO' )
  READ( *,'(a)' ) Filename
  Filename = ADJUSTL(Filename)
  IF ( .NOT. File_Exists( TRIM(Filename) ) ) THEN
    CALL Display_Message( PROGRAM_NAME,'File '//TRIM(Filename)//' not found.',FAILURE )
    STOP
  END IF

  ! Inquire the data file for the number of profiles
  err_stat = AtmProfile_InquireFile( Filename, n_Profiles = n_Profiles )
  IF ( err_stat /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME,'Error reading '//TRIM(Filename),FAILURE )
    STOP
  END IF

  ! Begin open loop for inspection
  Inspect_Loop: DO

    ! Enter profile to inspect
    WRITE( *,'(//4x,"Enter profile number [1-",i0,"] : ")',ADVANCE='NO') n_Profiles
    READ( *,* ) m
    ! ...Exit if requested profile out of range
    IF ( m < 1 .OR. m > n_Profiles ) EXIT Inspect_Loop

    ! Read the requested profile
    err_stat = AtmProfile_ReadFile( Filename, A, Profile_List=(/m/) )
    IF ( err_stat /= SUCCESS ) THEN
      WRITE( msg,'("Error reading profile #",i0," from ",a)' ) m, TRIM(Filename)
      CALL Display_Message( PROGRAM_NAME,msg,FAILURE )
      STOP
    END IF

    ! Display the contents
    CALL Inspect( A(1) )

  END DO Inspect_Loop

  ! Clean up
  CALL AtmProfile_Destroy( A )
  WRITE(*,'(/2x,"Goodbye...")')

END PROGRAM AtmProfile_Inspect

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm