CFPP$ UNROLL (32) F c c TITLE: implementation of finite-element method for model c in ``The Macroeconomics of War and Peace,'' mimeo, c Federal Reserve Bank of Minneapolis. c c PROBLEM: find c^l=c^l(k,kg,z), h^l=h^l(k,kg,z) that solve c the functional equation: c u_c(c^l,1-h^l,a^l) - beta sum_{j=1:n} pi(l,j) c * int u_c(c^j(x~),1-h^j(x~)-a^j) c *(1-delta+z~*F_k(k~,kg~,1-l(x~)-a^j) c *f(eps) d(eps), where x~=(k~,kg~,z~) c for l=1,..,n c k~ = z*F(k,kg,h^l(k,kg,z))-c^l(k,kg,z) c +(1-delta)k-ig^l-bg^l c z~ = z^psi*eps c c(0,0,z) = 0 c c ALGORITHM: the Galerkin method with piecewise trilinear shape c functions N_a(x,y,z), e.g., c c(x,y,z) = sum N_a(x,y,z)*c_a, c_a constant c c MODEL: stochastic growth model with government capital; c the government chooses c[t],i[t],h[t] to maximize c c E[ sum beta^t u(c[t],l[t]) | k[0],kg[0],z[0]] c t c subject to c c c[t]+i[t] = z[t]*f(k[t],kg[t],h[t])-ig[t]-bg[t] c k[t+1] = (1-delta)*k[t]+i[t] c kg[t+1] = (1-delta)*kg[t]+ig[t] c z[t+1] = z[t]^psi*eps[t] c c FUNCTIONS: u(c,l) = log(c)+alpha*log(l) c F(k,kg,h) = lam*z*(phi*k^rho+(1-phi)*kg^rho)^ ... c (theta/rho)*h^(1-theta) c c TEST CASE: delta=1, phi=1, alpha=0, ig[t]=0, bg[t]=0 c ==> c(k,kg,z)=(1-beta*theta)*lam*z*k^theta c c INPUTS: (a) war3da.inp: alpha,beta,delta,lam,phi,psi,rho,sig,theta c a0(i,j) i=1:na, j=1:n c xa(i),nxp(i) i=1:nx c xa(nx+1) c ya(i),nyp(i) i=1:ny c ya(ny+1) c za(i),nzp(i) i=1:nz c za(nz+1) c where a0's are initial guesses for policy rules c xa's,ya's,za's are grid pts. c nwp(i) are # of points per interval for c integrating over `w' interval, w=x,y,z c c OUTPUTS: (b) war3da.out: general results c (c) war3da.plt: x,y,z,c(x,y,z),n(x,y,z),k~(x,y,z),res c (d) war3da.nxt: new input file c (e) war3d.dat: input file for war3d.m (Matlab file) c c SUBROUTINES: (LINPACK) blas.f dgefa.f, dgesl.f, dgeco.f c (KMN) dgefs.f, machcon.f, xerror.f c (McGrattan) qgausl.f c c REFERENCES: Eischen, Jeff, (notes on the Finite Element Method) c Kahaner, Moler, and Nash (1989), NUMERICAL METHODS c AND SOFTWARE (KMN) c Saad, Youcef (1993), ``Lecture notes: Scientific c Computation 8002,'' Univ. Minnesota c c SEE ALSO: war3db.f (which solves Ax=b differently) c c Ellen McGrattan, 3-1-93 c Revised, 9-29-93, ERM c parameter (na=96,nx=15,ny=1,nz=2,ns=5,nas=480,nxm=3,nym=3, c nzm=3,ne1=10,nt=10000,nw=10000) implicit double precision (a-h,o-z) double precision f(nas), df(nas,nas), work(nw), c xa(nx+1), ya(ny+1), za(nz+1), xx(nxm,nx), wx(nxm,nx), c yy(nym,ny), wy(nym,ny), zz(nzm,nz), wz(nzm,nz), c x(nxm), u(nxm), y(nym), v(nym), z(nzm), w(nzm), c res(nt,ns), a(nt,ns), b(nt,ns), c(nt,ns), d(nt,ns), c a0(na,ns), a1(na,ns), zzt(ne1), e1(ne1), w1(ne1), c wgt(ne1), sum(9), wrk(ne1*ns,8), mcp(ns,ns), ginv(ns), c bg(ns), army(ns), lam integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), e(nx*ny*nz,8), c kk(ne1), ind(ne1*ns) c tiny = 0.000000001d0 pi = 3.14159265d0 crit = 0.0000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 open(unit=5, file='war3da.inp') open(unit=7, file='war3da.out') open(unit=8, file='war3da.plt') open(unit=9, file='war3da.nxt') open(unit=10,file='war3d.dat') read(5,*) alpha,beta,delta,lam,phi,psi,rho,sig,theta do 5 i=1,ns read(5,*) ginv(i),bg(i),army(i) 5 continue do 10 i=1,ns read(5,*) (mcp(i,j),j=1,ns) 10 continue do 20 i=1,na read(5,*) (a0(i,j),j=1,ns) 20 continue do 30 i=1,nx read(5,*) xa(i),nxp(i) 30 continue read(5,*) xa(nx1) do 40 i=1,ny read(5,*) ya(i),nyp(i) 40 continue read(5,*) ya(ny1) do 50 i=1,nz read(5,*) za(i),nzp(i) 50 continue read(5,*) za(nz1) c thet1= 1.d0-theta delt1= 1.d0-delta sum1 = 0.d0 sum2 = 0.d0 call qgausl(ne1,-2.88d0*sig,2.88d0*sig,e1,w1) do 60 i=1,ne1 wgt(i) = dexp(-0.5d0*e1(i)*e1(i)/(sig*sig))/ c dsqrt(2.d0*pi*sig*sig)*w1(i) sum1 = sum1 + wgt(i) 60 continue do 70 i=1,ne1 wgt(i) = wgt(i)/sum1 sum2 = sum2 + wgt(i) 70 continue n = 0 do 80 k=1,nz do 80 j=1,ny do 80 i=1,nx ll = (k-1)*nx*ny+(j-1)*nx+i e(ll,1) = (k-1)*nx1*ny1+(j-1)*nx1+i e(ll,2) = i e(ll,3) = j e(ll,4) = k e(ll,5) = nxp(i) e(ll,6) = nyp(j) e(ll,7) = nzp(k) e(ll,8) = n n = n+nxp(i)*nyp(j)*nzp(k) 80 continue do 90 i=1,nx call qgausl(nxp(i),xa(i),xa(i+1),x,u) do 90 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 90 continue do 100 i=1,ny call qgausl(nyp(i),ya(i),ya(i+1),y,v) do 100 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 100 continue do 110 i=1,nz call qgausl(nzp(i),za(i),za(i+1),z,w) do 110 j=1,nzp(i) zz(j,i) = z(j) wz(j,i) = w(j) 110 continue ne = nx*ny*nz dd = 1.d0 do 120 it=1,maxit do 130 i=1,nas f(i) = 0.d0 do 130 j=1,nas df(i,j) = 0.d0 130 continue do 140 n=1,ne n1 = e(n,5) n2 = e(n,6) n3 = e(n,7) i1 = e(n,8) x1 = xa(e(n,2)) x2 = xa(e(n,2)+1) dx = x2-x1 y1 = ya(e(n,3)) y2 = ya(e(n,3)+1) dy = y2-y1 z1 = za(e(n,4)) z2 = za(e(n,4)+1) dz = z2-z1 nd1= e(n,1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 150 j=1,n2 y(j) = yy(j,e(n,3)) 150 continue do 160 k=1,n3 z(k) = zz(k,e(n,4)) 160 continue c create vectors a,b,c for consumption, leisure, next capital do 170 i=1,n1 x(i) = xx(i,e(n,2)) xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx do 170 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 170 k=1,n3 zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz ll = i1+(k-1)*n1*n2+(j-1)*n1+i do 170 l=1,ns a(ll,l) = a0(nd1,l)*xp*yp*zp+a0(nd2,l)*xm*yp*zp+ c a0(nd3,l)*xm*ym*zp+a0(nd4,l)*xp*ym*zp+ c a0(nd5,l)*xp*yp*zm+a0(nd6,l)*xm*yp*zm+ c a0(nd7,l)*xm*ym*zm+a0(nd8,l)*xp*ym*zm a(ll,l) = dmax1(a(ll,l),tiny) c(ll,l) = z(k)*lam*(phi*x(i)**rho+(1.d0-phi)* c y(j)**rho)**(theta/rho) b(ll,l) = thet1*c(ll,l)/a(ll,l) d(ll,l) = -a(ll,l)+delt1*x(i)-ginv(l)-bg(l) 170 continue nn = n1*n2*n3 do 180 ll=i1+1,i1+nn do 180 l=1,ns tem1 = b(ll,l) b0 = 0.5d0 do 190 m=1,5 b0 = (tem1*(1.d0-army(l))-alpha*thet1*b0**theta)/ c (tem1+alpha*theta*b0**(theta-1.d0)) b0 = dmax1(b0,tiny) 190 continue b(ll,l) = dmin1(b0,1.d0) c(ll,l) = dmax1(c(ll,l)*b(ll,l)**thet1+d(ll,l),tiny) 180 continue 140 continue do 200 l=1,ns do 210 n=1,ne n1 = e(n,5) n2 = e(n,6) n3 = e(n,7) i1 = e(n,8) x1 = xa(e(n,2)) x2 = xa(e(n,2)+1) dx = x2-x1 y1 = ya(e(n,3)) y2 = ya(e(n,3)+1) dy = y2-y1 z1 = za(e(n,4)) z2 = za(e(n,4)+1) dz = z2-z1 nd1= e(n,1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 220 i=1,n1 x(i) = xx(i,e(n,2)) u(i) = wx(i,e(n,2)) 220 continue do 230 j=1,n2 y(j) = yy(j,e(n,3)) v(j) = wy(j,e(n,3)) 230 continue do 240 k=1,n3 z(k) = zz(k,e(n,4)) w(k) = wz(k,e(n,4)) zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz do 240 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 240 i=1,n1 ll = (k-1)*n1*n2+(j-1)*n1+i+i1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx xt = c(ll,l) yt = (1-delta)*y(j)+ginv(l) ii = 1 do 250 m=1,nx j1 = int( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 250 continue ii = max(1,ii-1) jj = 1 do 260 m=1,ny j1 = int( dmax1(yt-ya(m),0.d0)/(yt-ya(m))+.2) jj = jj + j1 260 continue jj = max(1,jj-1) CFPP$ UNROLL do 270 m=1,ne1 zzt(m) = dexp(psi*dlog(z(k))+e1(m)) kk(m) = 1 do 280 mm=1,nz j1 = int( dmax1(zzt(m)-za(mm),0.d0)/ c (zzt(m)-za(mm))+.2) kk(m) = kk(m) + j1 280 continue kk(m) = max(1,kk(m)-1) 270 continue do 290 m=1,9 sum(m) = 0.d0 290 continue do 300 m=1,nas work(m) = 0.d0 300 continue do 310 mm=1,ne1 do 310 m=1,ns nct = (m-1)*ne1+mm zt = zzt(mm) nn = (kk(mm)-1)*nx*ny+(jj-1)*nx+ii x1t = xa(e(nn,2)) x2t = xa(e(nn,2)+1) dxt = x2t-x1t y1t = ya(e(nn,3)) y2t = ya(e(nn,3)+1) dyt = y2t-y1t z1t = za(e(nn,4)) z2t = za(e(nn,4)+1) dzt = z2t-z1t xpt = (x2t-xt)/dxt xmt = (xt-x1t)/dxt ypt = (y2t-yt)/dyt ymt = (yt-y1t)/dyt zpt = (z2t-zt)/dzt zmt = (zt-z1t)/dzt nd1t = e(nn,1) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 nd5t = nd1t+nx1*ny1 nd6t = nd5t+1 nd8t = nd5t+nx1 nd7t = nd8t+1 ind(nct) = nd1t at = a0(nd1t,m)*xpt*ypt*zpt+ c a0(nd2t,m)*xmt*ypt*zpt+ c a0(nd3t,m)*xmt*ymt*zpt+ c a0(nd4t,m)*xpt*ymt*zpt+ c a0(nd5t,m)*xpt*ypt*zmt+ c a0(nd6t,m)*xmt*ypt*zmt+ c a0(nd7t,m)*xmt*ymt*zmt+ c a0(nd8t,m)*xpt*ymt*zmt at = dmax1(at,tiny) tem1 = zt*lam*(phi*xt**rho+(1.d0-phi)*yt**rho)** c (theta/rho)*thet1/at b0 = 0.5d0 do 320 mt=1,5 b0 = (tem1*(1.d0-army(m))-alpha*thet1* c b0**theta)/(tem1+alpha*theta* c b0**(theta-1.d0)) 320 continue bt = dmax1(b0,tiny) bt = dmin1(bt,1.d0) tem1 = phi*xt**rho+(1.d0-phi)*yt**rho tem2 = beta*zt*lam*tem1**(theta/rho)*bt**thet1 tem3 = xt**(rho-1.d0) tem4 = (-a0(nd1t,m)*ypt*zpt+a0(nd2t,m)*ypt*zpt c +a0(nd3t,m)*ymt*zpt-a0(nd4t,m)*ymt*zpt c -a0(nd5t,m)*ypt*zmt+a0(nd6t,m)*ypt*zmt c +a0(nd7t,m)*ymt*zmt-a0(nd8t,m)*ymt*zmt) c /dxt tem5 = 1.d0/(theta*(1.d0-bt-army(m))+bt) tem6 = -alpha*b(ll,l)/(theta*(1.d0- c b(ll,l)-army(l))+b(ll,l))-1.d0 tem7 = phi*theta*tem3*tem2/tem1*((rho-1.d0)/ c xt+rho*phi*(theta/rho-1.d0)*tem3/tem1+ c thet1*(1.d0-bt-army(m))*tem5*(phi*theta* c tem3/tem1-tem4/at))*tem6 tem9 = -phi*theta*tem3*tem2/tem1*thet1* c (1.d0-bt-army(m))*tem5/at tem1 = beta*delt1+tem2/tem1*phi*theta*tem3 tem2 = 1.d0/at tem8 = -tem2*tem2*tem4*tem6 tem0 = -tem2*tem2 sum(1) = sum(1) + mcp(l,m)*wgt(mm)*tem1*tem2 tem4 = mcp(l,m)*wgt(mm)*(tem7*tem2+tem1*tem8) sum(2) = sum(2) + tem4*xp*yp*zp sum(3) = sum(3) + tem4*xm*yp*zp sum(4) = sum(4) + tem4*xm*ym*zp sum(5) = sum(5) + tem4*xp*ym*zp sum(6) = sum(6) + tem4*xp*yp*zm sum(7) = sum(7) + tem4*xm*yp*zm sum(8) = sum(8) + tem4*xm*ym*zm sum(9) = sum(9) + tem4*xp*ym*zm c tem4 = mcp(l,m)*wgt(mm)*(tem9*tem2+tem1*tem0) wrk(nct,1)= tem4*xpt*ypt*zpt wrk(nct,2)= tem4*xmt*ypt*zpt wrk(nct,3)= tem4*xmt*ymt*zpt wrk(nct,4)= tem4*xpt*ymt*zpt wrk(nct,5)= tem4*xpt*ypt*zmt wrk(nct,6)= tem4*xmt*ypt*zmt wrk(nct,7)= tem4*xmt*ymt*zmt wrk(nct,8)= tem4*xpt*ymt*zmt 310 continue tem1= 1.d0/(a(ll,l)*a(ll,l)) m = (l-1)*na work(m+nd1) = sum(2)+tem1*xp*yp*zp work(m+nd2) = sum(3)+tem1*xm*yp*zp work(m+nd3) = sum(4)+tem1*xm*ym*zp work(m+nd4) = sum(5)+tem1*xp*ym*zp work(m+nd5) = sum(6)+tem1*xp*yp*zm work(m+nd6) = sum(7)+tem1*xm*yp*zm work(m+nd7) = sum(8)+tem1*xm*ym*zm work(m+nd8) = sum(9)+tem1*xp*ym*zm do 330 mm=1,ne1 do 330 m=1,ns nct = (m-1)*ne1+mm nd1t = ind(nct) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 nd5t = nd1t+nx1*ny1 nd6t = nd5t+1 nd8t = nd5t+nx1 nd7t = nd8t+1 mt = (m-1)*na work(mt+nd1t) = work(mt+nd1t) + wrk(nct,1) work(mt+nd2t) = work(mt+nd2t) + wrk(nct,2) work(mt+nd3t) = work(mt+nd3t) + wrk(nct,3) work(mt+nd4t) = work(mt+nd4t) + wrk(nct,4) work(mt+nd5t) = work(mt+nd5t) + wrk(nct,5) work(mt+nd6t) = work(mt+nd6t) + wrk(nct,6) work(mt+nd7t) = work(mt+nd7t) + wrk(nct,7) work(mt+nd8t) = work(mt+nd8t) + wrk(nct,8) 330 continue res(ll,l) = -1.d0/a(ll,l)+sum(1) tem1 = u(i)*v(j)*w(k)*res(ll,l) m = (l-1)*na f(m+nd1) = f(m+nd1) + tem1*xp*yp*zp f(m+nd2) = f(m+nd2) + tem1*xm*yp*zp f(m+nd3) = f(m+nd3) + tem1*xm*ym*zp f(m+nd4) = f(m+nd4) + tem1*xp*ym*zp f(m+nd5) = f(m+nd5) + tem1*xp*yp*zm f(m+nd6) = f(m+nd6) + tem1*xm*yp*zm f(m+nd7) = f(m+nd7) + tem1*xm*ym*zm f(m+nd8) = f(m+nd8) + tem1*xp*ym*zm c tem1 = u(i)*v(j)*w(k) do 340 mm=1,nas tem2 = tem1*work(mm) df(m+nd1,mm) = df(m+nd1,mm) + tem2*xp*yp*zp df(m+nd2,mm) = df(m+nd2,mm) + tem2*xm*yp*zp df(m+nd3,mm) = df(m+nd3,mm) + tem2*xm*ym*zp df(m+nd4,mm) = df(m+nd4,mm) + tem2*xp*ym*zp df(m+nd5,mm) = df(m+nd5,mm) + tem2*xp*yp*zm df(m+nd6,mm) = df(m+nd6,mm) + tem2*xm*yp*zm df(m+nd7,mm) = df(m+nd7,mm) + tem2*xm*ym*zm df(m+nd8,mm) = df(m+nd8,mm) + tem2*xp*ym*zm 340 continue 240 continue 210 continue 200 continue c do 344 j=1,nas c write(75,801) f(j) c do 344 i=1,nas c write(76,801) df(i,j) c344 continue c stop itask = 1 call dgefs (df, nas, nas, f, itask, ind1, work, iwork, rcond) sum3 = 0.d0 do 350 j=1,ns do 350 i=1,na tem1 = f((j-1)*na+i) a1(i,j) = a0(i,j) - tem1 sum3 = sum3 + tem1*tem1 350 continue dd = dsqrt(sum3)/float(nas) write(*,*) ' at iteration ', it, ' the residual is ', dd if (dd.lt.crit) go to 999 do 360 j=1,ns do 360 i=1,na a0(i,j) = a1(i,j) 360 continue 120 continue 999 continue c c print results c write(7,*) 'WAR-3D-A RESULTS' write(7,*) write(7,*) ' Parameters ' write(7,'('' alpha = '',F10.5,'' '')') alpha write(7,'('' beta = '',F10.5,'' '')') beta write(7,'('' delta = '',F10.5,'' '')') delta write(7,'('' lam = '',F10.5,'' '')') lam write(7,'('' phi = '',F10.5,'' '')') phi write(7,'('' psi = '',F10.5,'' '')') psi write(7,'('' rho = '',F10.5,'' '')') rho write(7,'('' sig = '',F10.5,'' '')') sig write(7,'('' theta = '',F10.5,'' '')') theta write(7,*) write(7,*) ' Markov Chain ' write(7,*) ' transition matrix: ' if (ns.eq.1) then write(7,801) mcp(1,1) else do 370 i=1,ns if (ns.eq.2) write(7,802) (mcp(i,j),j=1,ns) if (ns.eq.3) write(7,803) (mcp(i,j),j=1,ns) if (ns.eq.4) write(7,804) (mcp(i,j),j=1,ns) if (ns.eq.5) write(7,805) (mcp(i,j),j=1,ns) if (ns.eq.6) write(7,806) (mcp(i,j),j=1,ns) if (ns.eq.7) write(7,807) (mcp(i,j),j=1,ns) if (ns.eq.8) write(7,808) (mcp(i,j),j=1,ns) 370 continue endif write(7,*) write(7,*) ' values of ginv, bg, army in each state: ' do 380 i=1,ns write(7,803) ginv(i),bg(i),army(i) 380 continue write(7,*) write(7,*) ' Integrals: ' write(7,'('' sum1 = '',E9.3,'' '')') sum1 write(7,'('' sum2 = '',E9.3,'' '')') sum2 write(7,*) write(7,*) write(7,*) ' Results from dgefs.f ' if (ind1 .eq. -10) then write (7,*) ' error code =', ind1 else if (ind1 .lt. 0) then write (7,*) ' error code =', ind1 stop else write (7,*) ' number of accurate digits =', ind1 end if c write(7,*) write(7,*) ' Coefficients for Shape functions ' do 390 i=1,na if (ns.eq.1) write(7,801) (a0(i,j),j=1,ns) if (ns.eq.2) write(7,802) (a0(i,j),j=1,ns) if (ns.eq.3) write(7,803) (a0(i,j),j=1,ns) if (ns.eq.4) write(7,804) (a0(i,j),j=1,ns) if (ns.eq.5) write(7,805) (a0(i,j),j=1,ns) if (ns.eq.6) write(7,806) (a0(i,j),j=1,ns) if (ns.eq.7) write(7,807) (a0(i,j),j=1,ns) if (ns.eq.8) write(7,808) (a0(i,j),j=1,ns) 390 continue write(7,*) write(7,*) ' Grid Points (xa,ya,za) ' write(7,*) ' xa ya za ' write(7,*) ' ------------------------------ ' do 400 k=1,nz1 do 401 j=1,ny1 do 402 i=1,nx1 write(7,813) xa(i),ya(j),za(k) 402 continue 401 continue 400 continue write(7,*) write(7,*) ' Number of iterations = ', it c c write file to with values of x,y,z,c(x,y,z),h(x,y,z),k'(x,y,z),res: c do 410 l=1,ns do 410 k1=1,nz do 410 k=1,nzp(k1) do 410 j1=1,ny mm = e((k1-1)*nx*ny+(j1-1)*nx+1,8) do 410 j=1,nyp(j1) nn = mm do 410 i1=1,nx nn = nn+(k-1)*nxp(i1)*nyp(j1)+(j-1)*nxp(i1) do 420 i=1,nxp(i1) nn = nn + 1 write(8,807) xx(i,i1),yy(j,j1),zz(k,k1),a(nn,l), c b(nn,l),c(nn,l),res(nn,l) 420 continue nn = nn+(nzp(k1)-k)*nxp(i1)*nyp(j1)+(nyp(j1)-j)* c nxp(i1) 410 continue c c write file to be used as next input file: c write(9,808) alpha,beta,delta,lam,phi,psi,rho,sig,theta do 430 i=1,ns write(9,803) ginv(i),bg(i),army(i) 430 continue do 440 i=1,ns write(9,807) (mcp(i,j),j=1,ns) 440 continue do 450 i=1,na write(9,807) (a0(i,j),j=1,ns) 450 continue do 460 i=1,nx write(9,811) xa(i),nxp(i) 460 continue write(9,801) xa(nx1) do 470 i=1,ny write(9,811) ya(i),nyp(i) 470 continue write(9,801) ya(ny1) do 480 i=1,nz write(9,811) za(i),nzp(i) 480 continue write(9,801) za(nz1) c c write file to be read into matlab: c write(10,801) alpha write(10,801) beta write(10,801) delta write(10,801) lam write(10,801) phi write(10,801) psi write(10,801) rho write(10,801) sig write(10,801) theta write(10,812) ns do 490 i=1,ns write(10,801) ginv(i) write(10,801) bg(i) write(10,801) army(i) 490 continue do 500 j=1,ns do 500 i=1,ns write(10,801) mcp(i,j) 500 continue write(10,812) nx1 write(10,812) ny1 write(10,812) nz1 do 510 i=1,nx1 write(10,801) xa(i) 510 continue do 520 i=1,ny1 write(10,801) ya(i) 520 continue do 530 i=1,nz1 write(10,801) za(i) 530 continue do 540 j=1,ns do 540 i=1,na write(10,801) a0(i,j) 540 continue c 801 format (1X, F12.6) 802 format (2(1X, F12.6)) 803 format (3(1X, F12.6)) 804 format (4(1X, F12.6)) 805 format (5(1X, F12.6)) 806 format (6(1X, F11.5)) 807 format (7(1X, F10.5)) 808 format (8(1X, F9.5)) 811 format (1X, F12.6, 4X, I4) 812 format (9X, I4) 813 format (4X, 3(1X, F12.6)) stop end