CFPP$ UNROLL (32) F c c PROGRAM: CUMC3DBi.F c c TITLE: implementation of the finite-element method for model c with capacity utilization and employment costs, in c ``The Effects of Government purchases on Employment c and Hours,'' (Toni Braun, Ellen McGrattan). c c PROBLEM: find c1^l(k,e,z), n^l(k,e,z), c0^l(k,e,z), and c h^l(k,e,z) that solve the functional equations: c c U_c(c1^l,1-h^l)(1-a^l) - beta sum_{j=1:n} pi(l,j) c * int U_c(c1^l(x~),1-h^l(x~))(1-a^j) c *(1-delta+z~*F_k(k~,n^j(x~),h^j(x~))) c *f(eps) d(eps) = 0 (1) c U(c1^l,1-h^l) - U(c0^l,1)-U_c(c1^l,1-h^l) c *(c1^l-c0^l-z*F(k,n^l,h^l)+2*alpha c *(n^l-e)-P'(n^l))*(1-a^l)-2*alpha*beta c *sum_{j=1:n} pi(l,j) int (1-a^j) c *U_c(c1^j(x~),1-h^j(x~))*(n^j(x~)-e~) c *f(eps) d(eps) = 0, (2) c U_c(c0^l,1) - U_c(c1^l,1-h^l) = 0, (3) c n^l*U_l(c1^l,1-h^l) - U_c(c1^l,1-h^l)*F_h(k,n^l,h^l)=0 (4) c c for l=1,..,ns, where x~=(k~,e~,z~) and c k~ = z*F(k,n^l(k,e,z),h^l(k,e,z))-n^l(k,e,z)* c c1^l(k,e,z)-(1-n^l(k,e,z))*c0^l(k,e,z) c +(1-delta)*k-(cg^l+zeta*log(z)) c -alpha(n^l(k,e,z)-e)^2-P(n^l(k,e,z)) c e~ = n^l(k,e,z) c z~ = z^psi*eps c c1(0,e,z) = 0 c c ALGORITHM: the Galerkin method with piecewise trilinear shape c functions N_a(x,y,z), e.g., c c1(x,y,z) = sum N_a(x,y,z)*c_a, c_a constant c n(x,y,z) = sum N_a(x,y,z)*n_a, n_a constant c (NOTE: LQ solution used as initial guess for c_a's,n_a's) c c MODEL: stochastic growth model with government capital, c varying hours and employment, and moving costs; c the planner chooses c1[t],c0[t],i[t],n[t],h[t] c to maximize c c E[ sum beta^t {n[t]*U(c1[t],1-h[t])+(1-n[t]) c t *U(c0[t],1)}*(1-a[t])|x[0]] c c subject to c c c[t]+i[t]+m[t]+p[t] = z[t]*F(k[t]+kg[t],n[t],h[t]) c -ig[t]-cg[t] c c[t] = n[t]*c1[t]+(1-n[t])*c0[t] c m[t] = M(n[t],n[t-1]) c p[t] = P(n[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 x[t] = (k[t],n[t-1],z[t]) c c FUNCTIONS: U(c,l) = [(c^gamma*l^(1-gamma))^(1-omega)-1]/(1-omega) c F(k,kg,n,h) = lam*z*(k+kg)^theta*n^(1-theta)*h c M(n,e) = alpha*(n-e)*(n-e) c P(n) = [(n-pa)*pb]^3+abs((n-pa)*pb)^3 c c TEST CASE: set alpha=0 and compare results to Kydland-Prescott c c INPUTS: (a) cumc3dai.inp: alpha,beta,delta,gamma,lam,omega, ... c pa,pb,pc,psi,sig,theta,zeta c ig(j),cg(j),a(j), j=1:ns c mcp(i,j) i=1:ns, j=1:ns c ssk,ssc,ssn,ssh c where alpha, etc. are parameters of utility/technology c ig's,cg's,a's are values for Markov chain c mcp's are the transition probabilities c ssk, etc. are initial guesses of the steady c for k, c, n, h (if 0, then program chooses) c c OUTPUTS: (a) cumc3d.dat: input file for cumc3d.m (Matlab file) c c OTHER (LINPACK) blas.f, dgeco.f, dgedi.f dgefa.f, dgesl.f c SUBROUTINES: (McGrattan) qgausl.f c c REFERENCES: Eischen, Jeff, (notes on the Finite Element Method) c Saad, Youcef (1993), ``Lecture notes: Scientific c Computation 8002,'' Univ. Minnesota c c SEE ALSO: cumc3da.f or cumc3db.f c c Ellen McGrattan, 10-24-93 c Revised, 10-24-93, ERM c parameter (na=126,nx=13,ny=2,nz=2,ns=4,nas=1008,nxm=3,nym=3,nzm=3, c ne1=10,nw=20000,nsv=4,ndv=4,nv=8,nt=1000) 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 a0(na,ns), a1(na,ns), b0(na,ns), b1(na,ns), zzt(ne1), c e1(ne1), w1(ne1), wgt(ne1), sum(34), wrk(ne1*ns,32), c mcp(ns,ns), ginv(ns), cg(ns), army(ns), lam, n0, n0t, c mu, wrk2(nas), foc(4), dfoc(4,4), del(8,8), c mom, amat(nsv,nsv), bmat(nsv,ndv), qmat(nsv,nsv), c rmat(ndv,ndv), wmat(nsv,ndv), aa(nsv,nsv), bb(nsv,nsv), c cc(nsv,nsv), dd(nsv,nsv), ee(nsv,nsv), qwr(nv,nv), c f0(ndv,nsv), f1(ndv,nsv), du(8), du2(8,8), ri(ndv,ndv), c rw(ndv,nsv), rb(ndv,ndv), det(2) integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), kk(ne1), c ind(ne1*ns), istat(nt) c open(unit=5, file='cumc3dbi.inp') open(unit=10, file='cumc3db.dat') c read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) gamma read(5,*) lam read(5,*) omega read(5,*) pa read(5,*) pb read(5,*) pc read(5,*) psi read(5,*) sig read(5,*) theta read(5,*) zeta read(5,*) itype do 1 i=1,ns read(5,*) ginv(i),cg(i),army(i) 1 continue do 5 i=1,ns read(5,*) (mcp(i,j),j=1,ns) 5 continue read(5,*) ssk,ssc,sse,ssh read(5,*) lqonly c c 3-dimensional domain (xa,ya,za): c nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 read(5,*) irgrid if (irgrid.eq.1) then do 7 i=1,nx1 read(5,*) xa(i) 7 continue do 8 i=1,ny1 read(5,*) ya(i) 8 continue do 9 i=1,nz1 read(5,*) za(i) 9 continue else do 11 i=1,nx1 xa(i) = 0.27d0+0.43d0*float(i)/float(nx1) 11 continue do 12 i=1,ny1 ya(i) = 0.2d0+0.45d0*float(i)/float(ny1) 12 continue do 13 i=1,nz1 za(i) = 0.8d0+0.3d0*float(i)/float(nz1) 13 continue endif c c initial decision functions (optional) c read(5,*) irfunc if (irfunc.eq.1) then do 14 i=1,na read(5,*) (a0(i,j),j=1,ns) 14 continue do 15 i=1,na read(5,*) (b0(i,j),j=1,ns) 15 continue endif do 10 i=1,nx nxp(i) = 2 10 continue do 20 i=1,ny nyp(i) = 2 20 continue do 30 i=1,nz nzp(i) = 2 30 continue delt1 = 1.d0-delta gamm1 = 1.d0-gamma ome1 = 1.d0-omega thet1 = 1.d0-theta cwgt = gamma*ome1 cwgt1 = cwgt-1.d0 hwgt = gamm1*ome1 hwgt1 = hwgt-1.d0 tiny = 0.000000001d0 pi = 3.14159265d0 crit = 0.0000001d0 maxit = 20 maxit0 = 20 if (lqonly.eq.1) maxit0 = 0 if (irfunc.eq.0) then do 40 i=1,8 do 40 j=1,8 del(i,j) =0.d0 40 continue if (ssk.lt.tiny) ssk = 0.4 if (ssc.lt.tiny) ssc = 0.1 if (sse.lt.tiny) sse = 0.4 if (ssh.lt.tiny) ssh = 0.42 do 50 l=1,ns del(1,1) = ssk*1.0e-5 del(2,2) = ssc*1.0e-5 del(3,3) = sse*1.0e-5 del(4,4) = ssh*1.0e-5 do 60 t=1,10 if (itype.eq.1) then pn = pa*sse+pb*sse*sse+pc dpn = pa+2.0*pb*sse endif if (itype.eq.2) then d = (sse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d dpn = 3.d0*d*(d+dsign(d,0.d0))*pb+pc endif if (itype.eq.3) then pn = (dexp(pa*sse)-1.d0)*pb+pc dpn = pa*pb*dexp(pa*sse) endif zkk = lam*ssk**theta c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) foc(1) = gamm1/gamma*sse**theta*ssc-zkk*(1.d0-ssh) foc(2) = ssc**cwgt*(1.d0-ssh)**hwgt/ome1-c0**cwgt/ome1- c gamma*ssc**cwgt1*(1.d0-ssh)**hwgt*(ssc-c0- c zkk*thet1*sse**(-theta)*ssh+dpn) foc(3) = sse*ssc+(1.d0-sse)*c0+delta*ssk+cg(l)- c zkk*ssh*sse**thet1+pn foc(4) = beta*(delt1+theta*zkk/ssk*sse**thet1*ssh)-1.d0 do 70 i=1,4 dssk = ssk + del(1,i) dssc = ssc + del(2,i) dsse = sse + del(3,i) dssh = ssh + del(4,i) zkk = lam*dssk**theta c0 = dssc*(1.d0-dssh)**(hwgt/cwgt1) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc dpn = pa+2.0*pb*dsse endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d dpn = 3.d0*d*(d+dsign(d,0.d0))*pb+pc endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc dpn = pa*pb*dexp(pa*dsse) endif foc1p = gamm1/gamma*dsse**theta*dssc-zkk*(1.d0-dssh) foc2p = dssc**cwgt*(1.d0-dssh)**hwgt/ome1-c0**cwgt/ome1- c gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh+dpn) foc3p = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1+pn foc4p = beta*(delt1+theta*zkk/dssk*dsse**thet1*dssh)-1.d0 dssk = ssk - del(1,i) dssc = ssc - del(2,i) dsse = sse - del(3,i) dssh = ssh - del(4,i) zkk = lam*dssk**theta c0 = dssc*(1.d0-dssh)**(hwgt/cwgt1) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc dpn = pa+2.0*pb*dsse endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d dpn = 3.d0*d*(d+dsign(d,0.d0))*pb+pc endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc dpn = pa*pb*dexp(pa*dsse) endif foc1m = gamm1/gamma*dsse**theta*dssc-zkk*(1.d0-dssh) foc2m = dssc**cwgt*(1.d0-dssh)**hwgt/ome1-c0**cwgt/ome1- c gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh+dpn) foc3m = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1+pn foc4m = beta*(delt1+theta*zkk/dssk*dsse**thet1*dssh)-1.d0 dfoc(1,i) = 0.5*(foc1p-foc1m)/del(i,i) dfoc(2,i) = 0.5*(foc2p-foc2m)/del(i,i) dfoc(3,i) = 0.5*(foc3p-foc3m)/del(i,i) dfoc(4,i) = 0.5*(foc4p-foc4m)/del(i,i) 70 continue call dgeco(dfoc,4,4,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'steady state values' call dgesl(dfoc,4,4,iwork,foc,0) ssk = ssk - foc(1) ssc = ssc - foc(2) sse = sse - foc(3) ssh = ssh - foc(4) 60 continue ssi = delta*ssk ssn = sse ssz = 0.d0 ssf = lam*ssk**theta*sse**thet1*ssh del(1,1) = ssk*1.0e-4 del(2,2) = ssn*1.0e-4 del(3,3) = 1.0e-5 del(4,4) = 1.0e-4 del(5,5) = ssc*1.0e-4 del(6,6) = ssi*1.0e-4 del(7,7) = sse*1.0e-4 del(8,8) = ssh*1.0e-4 if (itype.eq.1) then pn = pa*sse+pb*sse*sse+pc endif if (itype.eq.2) then d = (sse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*sse)-1.d0)*pb+pc endif c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) uu = (1.d0-army(l))*((sse*ssc**cwgt*(1.d0-ssh)**hwgt+ c (1.d0-sse)*c0**cwgt)/ome1) du(4) = 0.d0 do 80 i=1,8 du(i) = 0.d0 do 80 j=1,8 du2(i,j) = 0.d0 80 continue do 90 i=1,8 dssk = ssk + del(1,i) dssn = ssn + del(2,i) dssz = ssz + del(3,i) dssc = ssc + del(5,i) dssi = ssi + del(6,i) dsse = sse + del(7,i) dssh = ssh + del(8,i) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc endif c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)*(dsse-dssn)- c pn)/(1.d0-dsse) up = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) dssk = ssk - del(1,i) dssn = ssn - del(2,i) dssz = ssz - del(3,i) dssc = ssc - del(5,i) dssi = ssi - del(6,i) dsse = sse - del(7,i) dssh = ssh - del(8,i) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc endif c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)*(dsse-dssn)- c pn)/(1.d0-dsse) um = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) du(i)= 0.5*(up-um)/del(i,i) do 100 j=1,i dssk = ssk + del(1,i) + del(1,j) dssn = ssn + del(2,i) + del(2,j) dssz = ssz + del(3,i) + del(3,j) dssc = ssc + del(5,i) + del(5,j) dssi = ssi + del(6,i) + del(6,j) dsse = sse + del(7,i) + del(7,j) dssh = ssh + del(8,i) + del(8,j) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc endif c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn)-pn)/(1.d0-dsse) upp = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) dssk = ssk + del(1,i) - del(1,j) dssn = ssn + del(2,i) - del(2,j) dssz = ssz + del(3,i) - del(3,j) dssc = ssc + del(5,i) - del(5,j) dssi = ssi + del(6,i) - del(6,j) dsse = sse + del(7,i) - del(7,j) dssh = ssh + del(8,i) - del(8,j) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc endif c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn)-pn)/(1.d0-dsse) upm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) dssk = ssk - del(1,i) + del(1,j) dssn = ssn - del(2,i) + del(2,j) dssz = ssz - del(3,i) + del(3,j) dssc = ssc - del(5,i) + del(5,j) dssi = ssi - del(6,i) + del(6,j) dsse = sse - del(7,i) + del(7,j) dssh = ssh - del(8,i) + del(8,j) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc endif c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn)-pn)/(1.d0-dsse) ump = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) dssk = ssk - del(1,i) - del(1,j) dssn = ssn - del(2,i) - del(2,j) dssz = ssz - del(3,i) - del(3,j) dssc = ssc - del(5,i) - del(5,j) dssi = ssi - del(6,i) - del(6,j) dsse = sse - del(7,i) - del(7,j) dssh = ssh - del(8,i) - del(8,j) if (itype.eq.1) then pn = pa*dsse+pb*dsse*dsse+pc endif if (itype.eq.2) then d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*dsse)-1.d0)*pb+pc endif c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn)-pn)/(1.d0-dsse) umm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) du2(i,j) = 0.25*(upp-upm-ump+umm)/(del(i,i)*del(j,j)) du2(j,i) = du2(i,j) 100 continue 90 continue do 110 i=1,8 do 110 j=1,8 qwr(i,j) = 0.d0 110 continue do 120 j=1,8 qwr(1,j) = 0.5*du2(1,j) qwr(2,j) = 0.5*du2(2,j) qwr(3,j) = 0.5*du2(3,j) qwr(4,j) = 0.5*du2(4,j) qwr(5,j) = 0.5*du2(5,j) qwr(6,j) = 0.5*du2(6,j) qwr(7,j) = 0.5*du2(7,j) qwr(8,j) = 0.5*du2(8,j) 120 continue do 130 i=1,8 work(i) = 0.5*(du(i)-du2(i,1)*ssk-du2(i,2)*ssn-du2(i,5)* c ssc-du2(i,6)*ssi-du2(i,7)*sse-du2(i,8)*ssh) qwr(4,i) = qwr(4,i) + work(i) qwr(i,4) = qwr(i,4) + work(i) 130 continue qwr(4,4) = qwr(4,4) + uu-(work(1)+0.5*du(1))*ssk- c (work(2)+0.5*du(2))*ssn-(work(5)+0.5*du(5))*ssc- c (work(6)+0.5*du(6))*ssi-(work(7)+0.5*du(7))*sse- c (work(8)+0.5*du(8))*ssh do 140 i=1,4 do 140 j=1,4 qmat(i,j) = qwr(i,j) amat(i,j) = 0.d0 140 continue do 150 i=1,4 do 150 j=1,4 wmat(i,j) = qwr(i,j+4) bmat(i,j) = 0.d0 150 continue do 160 i=1,4 do 160 j=1,4 rmat(i,j) = qwr(i+4,j+4) 160 continue amat(1,1) = delt1 amat(3,3) = psi amat(4,4) = 1.d0 bmat(1,2) = 1.d0 bmat(2,3) = 1.d0 call double(beta,amat,bmat,qmat,rmat,wmat,4,4,icode,f1) if (icode.eq.-1) then write(*,*) 'ERROR: Parameters ncx,nsx too small in double.' stop endif n = 1 do 630 k=1,nz1 do 630 j=1,ny1 do 630 i=1,nx1 a0(n,l) = f1(1,1)*xa(i)+f1(1,2)*ya(j)+ c f1(1,3)*dlog(za(k))+f1(1,4) b0(n,l) = f1(3,1)*xa(i)+f1(3,2)*ya(j)+ c f1(3,3)*dlog(za(k))+f1(3,4) n = n+1 630 continue 50 continue endif c sum1 = 0.d0 sum2 = 0.d0 call qgausl(ne1,-2.88d0*sig,2.88d0*sig,e1,w1) do 640 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) 640 continue do 650 i=1,ne1 wgt(i) = wgt(i)/sum1 sum2 = sum2 + wgt(i) 650 continue do 660 i=1,nx call qgausl(nxp(i),xa(i),xa(i+1),x,u) do 660 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 660 continue do 670 i=1,ny call qgausl(nyp(i),ya(i),ya(i+1),y,v) do 670 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 670 continue do 680 i=1,nz call qgausl(nzp(i),za(i),za(i+1),z,w) do 680 j=1,nzp(i) zz(j,i) = z(j) wz(j,i) = w(j) 680 continue ne = nx*ny*nz nxy = nx*ny sumres = 1.d0 do 690 it=1,maxit0 do 700 i=1,nas f(i) = 0.d0 do 700 j=1,nas df(i,j) = 0.d0 700 continue do 710 l=1,ns army1 = 1.d0-army(l) do 720 n=1,ne iz = (n-1)/nxy+1 iy = (n-(iz-1)*nxy-1)/nx+1 ix = n-(iy-1)*nx-(iz-1)*nxy n1 = nxp(ix) n2 = nyp(iy) n3 = nzp(iz) x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 z1 = za(iz) z2 = za(iz+1) dz = z2-z1 nd1= ix+nx1*(iy-1)+nx1*ny1*(iz-1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 730 i=1,n1 x(i) = xx(i,ix) u(i) = wx(i,ix) 730 continue do 740 j=1,n2 y(j) = yy(j,iy) v(j) = wy(j,iy) 740 continue do 750 k=1,n3 z(k) = zz(k,iz) w(k) = wz(k,iz) zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz do 750 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 750 i=1,n1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c1 = 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 n0 = b0(nd1,l)*xp*yp*zp+b0(nd2,l)*xm*yp*zp+ c b0(nd3,l)*xm*ym*zp+b0(nd4,l)*xp*ym*zp+ c b0(nd5,l)*xp*yp*zm+b0(nd6,l)*xm*yp*zm+ c b0(nd7,l)*xm*ym*zm+b0(nd8,l)*xp*ym*zm h = 1-gamm1*c1*(n0/x(i))**theta/(z(k)*lam*gamma) c0 = c1*(1.d0-h)**(hwgt/cwgt1) h1 = 1.d0-h mu = gamma*c1**cwgt1*h1**hwgt cg1 = cg(l)+zeta*dlog(z(k)) if (itype.eq.1) then pn = pa*n0+pb*n0*n0+pc endif if (itype.eq.2) then d = (n0-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d endif if (itype.eq.3) then pn = (dexp(pa*n0)-1.d0)*pb+pc endif xt = delt1*x(i)+z(k)*lam*x(i)**theta*n0**thet1*h- c cg1-n0*c1-(1.d0-n0)*c0-alpha*(n0-y(j))* c (n0-y(j))-pn yt = n0 ii = 1 CFPP$ UNROLL do 760 m=1,nx j1 = idint( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 760 continue ii = max(1,ii-1) jj = 1 CFPP$ UNROLL do 770 m=1,ny j1 = idint( dmax1(yt-ya(m),0.d0)/(yt-ya(m))+.2) jj = jj + j1 770 continue jj = max(1,jj-1) CFPP$ UNROLL do 780 m=1,ne1 zzt(m) = dexp(psi*dlog(z(k))+e1(m)) kk(m) = 1 CFPP$ UNROLL do 790 mm=1,nz j1 = idint( dmax1(zzt(m)-za(mm),0.d0)/ c (zzt(m)-za(mm))+.2) kk(m)= kk(m) + j1 790 continue kk(m) = max(1,kk(m)-1) 780 continue do 800 m=1,34 sum(m) = 0.d0 800 continue do 810 m=1,nas work(m) = 0.d0 wrk2(m) = 0.d0 810 continue if (itype.eq.1) then pn = pa*n0+pb*n0*n0+pc dpn = pa+2.0*pb*n0 d2pn = 2.0*pb endif if (itype.eq.2) then d = (n0-pa)*pb pn = d*d*d+dabs(d*d*d)+pc*d dpn = 3.d0*d*(d+dsign(d,0.d0))*pb+pc d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb endif if (itype.eq.3) then pn = (dexp(pa*n0)-1.d0)*pb+pc dpn = pa*pb*dexp(pa*n0) d2pn = pa*pa*pb*dexp(pa*n0) endif dc0ca = -c0/c1*omega/cwgt1 dc0na = hwgt/cwgt1*c0*theta/n0 dhca = -h1/c1 dhna = -theta*h1/n0 dxca = z(k)*lam*x(i)**theta*n0**thet1*dhca-n0- c (1.d0-n0)*dc0ca dxna = z(k)*lam*x(i)**theta*n0**thet1*h*(thet1/ c n0+dhna/h)-c1+c0-(1.d0-n0)*dc0na-2.d0* c alpha*(n0-y(j))-dpn CFPP$ UNROLL do 820 m=1,ns CFPP$ UNROLL do 820 mm=1,ne1 nct = (m-1)*ne1+mm zt = zzt(mm) x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t y1t = ya(jj) y2t = ya(jj+1) dyt = y2t-y1t z1t = za(kk(mm)) z2t = za(kk(mm)+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 = (kk(mm)-1)*nx1*ny1+(jj-1)*nx1+ii 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 c1t = 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 n0t = b0(nd1t,m)*xpt*ypt*zpt+ c b0(nd2t,m)*xmt*ypt*zpt+ c b0(nd3t,m)*xmt*ymt*zpt+ c b0(nd4t,m)*xpt*ymt*zpt+ c b0(nd5t,m)*xpt*ypt*zmt+ c b0(nd6t,m)*xmt*ypt*zmt+ c b0(nd7t,m)*xmt*ymt*zmt+ c b0(nd8t,m)*xpt*ymt*zmt ht = 1.d0-gamm1*c1t*(n0t/xt)**theta/ c (gamma*zt*lam) dc1tx = (-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 dc1ty = (-a0(nd1t,m)*xpt*zpt-a0(nd2t,m)*xmt*zpt c +a0(nd3t,m)*xmt*zpt+a0(nd4t,m)*xpt*zpt c -a0(nd5t,m)*xpt*zmt-a0(nd6t,m)*xmt*zmt c +a0(nd7t,m)*xmt*zmt+a0(nd8t,m)*xpt*zmt) c /dyt dn0tx = (-b0(nd1t,m)*ypt*zpt+b0(nd2t,m)*ypt*zpt c +b0(nd3t,m)*ymt*zpt-b0(nd4t,m)*ymt*zpt c -b0(nd5t,m)*ypt*zmt+b0(nd6t,m)*ypt*zmt c +b0(nd7t,m)*ymt*zmt-b0(nd8t,m)*ymt*zmt) c /dxt dn0ty = (-b0(nd1t,m)*xpt*zpt-b0(nd2t,m)*xmt*zpt c +b0(nd3t,m)*xmt*zpt+b0(nd4t,m)*xpt*zpt c -b0(nd5t,m)*xpt*zmt-b0(nd6t,m)*xmt*zmt c +b0(nd7t,m)*xmt*zmt+b0(nd8t,m)*xpt*zmt) c /dyt out = beta*zt*lam*xt**theta*n0t**thet1*ht h1t = 1.d0-ht term1 = gamma*c1t**cwgt1*h1t**hwgt*(1.d0-army(m)) term2 = beta*delt1+theta*out/xt term3 = 2.d0*alpha*beta*(n0t-yt) dc1tca = dc1tx*dxca dn0tca = dn0tx*dxca dc1tna = dc1tx*dxna+dc1ty dn0tna = dn0tx*dxna+dn0ty dhtca =-h1t*(dc1tca/c1t+theta*dn0tca/n0t-theta* c dxca/xt) dhtna =-h1t*(dc1tna/c1t+theta*dn0tna/n0t-theta* c dxna/xt) dhtcb =-h1t/c1t dhtnb =-theta*h1t/n0t dt1ca = term1*(cwgt1*dc1tca/c1t-hwgt*dhtca/h1t) dt1cb = term1*(cwgt1/c1t-hwgt*dhtcb/h1t) dt1na = term1*(cwgt1*dc1tna/c1t-hwgt*dhtna/h1t) dt1nb =-term1*hwgt*dhtnb/h1t dt2ca = theta*out/xt*(-thet1*dxca/xt+thet1* c dn0tca/n0t+dhtca/ht) dt2cb = theta*out/xt*dhtcb/ht dt2na = theta*out/xt*(-thet1*dxna/xt+thet1* c dn0tna/n0t+dhtna/ht) dt2nb = theta*out/xt*(thet1/n0t+dhtnb/ht) dt3ca = 2.d0*alpha*beta*dn0tca dt3cb = 0.d0 dt3na = 2.d0*alpha*beta*(dn0tna-1.d0) dt3nb = 2.d0*alpha*beta sum(1) = sum(1) + mcp(l,m)*wgt(mm)*term1*term2 sum(2) = sum(2) + mcp(l,m)*wgt(mm)*term1*term3 drhsc = mcp(l,m)*wgt(mm)*(dt1ca*term2+ c term1*dt2ca) drhsn = mcp(l,m)*wgt(mm)*(dt1na*term2+ c term1*dt2na) c rhs of res1,res2 c d(res1) wrt a0(nd1),b0(nd1), .... sum(3) = sum(3) + drhsc*xp*yp*zp sum(4) = sum(4) + drhsn*xp*yp*zp sum(5) = sum(5) + drhsc*xm*yp*zp sum(6) = sum(6) + drhsn*xm*yp*zp sum(7) = sum(7) + drhsc*xm*ym*zp sum(8) = sum(8) + drhsn*xm*ym*zp sum(9) = sum(9) + drhsc*xp*ym*zp sum(10)= sum(10) + drhsn*xp*ym*zp sum(11)= sum(11) + drhsc*xp*yp*zm sum(12)= sum(12) + drhsn*xp*yp*zm sum(13)= sum(13) + drhsc*xm*yp*zm sum(14)= sum(14) + drhsn*xm*yp*zm sum(15)= sum(15) + drhsc*xm*ym*zm sum(16)= sum(16) + drhsn*xm*ym*zm sum(17)= sum(17) + drhsc*xp*ym*zm sum(18)= sum(18) + drhsn*xp*ym*zm drhsc = mcp(l,m)*wgt(mm)*(dt1ca*term3+ c term1*dt3ca) drhsn = mcp(l,m)*wgt(mm)*(dt1na*term3+ c term1*dt3na) c d(res2) wrt a0(nd1),b0(nd1), .... sum(19)= sum(19) + drhsc*xp*yp*zp sum(20)= sum(20) + drhsn*xp*yp*zp sum(21)= sum(21) + drhsc*xm*yp*zp sum(22)= sum(22) + drhsn*xm*yp*zp sum(23)= sum(23) + drhsc*xm*ym*zp sum(24)= sum(24) + drhsn*xm*ym*zp sum(25)= sum(25) + drhsc*xp*ym*zp sum(26)= sum(26) + drhsn*xp*ym*zp sum(27)= sum(27) + drhsc*xp*yp*zm sum(28)= sum(28) + drhsn*xp*yp*zm sum(29)= sum(29) + drhsc*xm*yp*zm sum(30)= sum(30) + drhsn*xm*yp*zm sum(31)= sum(31) + drhsc*xm*ym*zm sum(32)= sum(32) + drhsn*xm*ym*zm sum(33)= sum(33) + drhsc*xp*ym*zm sum(34)= sum(34) + drhsn*xp*ym*zm c drhsc = mcp(l,m)*wgt(mm)*(dt1cb*term2+ c term1*dt2cb) drhsn = mcp(l,m)*wgt(mm)*(dt1nb*term2+ c term1*dt2nb) c d(res1) wrt a0(nd1t),b0(nd1t), ... wrk(nct,1) = drhsc*xpt*ypt*zpt wrk(nct,2) = drhsn*xpt*ypt*zpt wrk(nct,3) = drhsc*xmt*ypt*zpt wrk(nct,4) = drhsn*xmt*ypt*zpt wrk(nct,5) = drhsc*xmt*ymt*zpt wrk(nct,6) = drhsn*xmt*ymt*zpt wrk(nct,7) = drhsc*xpt*ymt*zpt wrk(nct,8) = drhsn*xpt*ymt*zpt wrk(nct,9) = drhsc*xpt*ypt*zmt wrk(nct,10)= drhsn*xpt*ypt*zmt wrk(nct,11)= drhsc*xmt*ypt*zmt wrk(nct,12)= drhsn*xmt*ypt*zmt wrk(nct,13)= drhsc*xmt*ymt*zmt wrk(nct,14)= drhsn*xmt*ymt*zmt wrk(nct,15)= drhsc*xpt*ymt*zmt wrk(nct,16)= drhsn*xpt*ymt*zmt drhsc = mcp(l,m)*wgt(mm)*(dt1cb*term3+ c term1*dt3cb) drhsn = mcp(l,m)*wgt(mm)*(dt1nb*term3+ c term1*dt3nb) c d(res2) wrt a0(nd1t),b0(nd1t), ... wrk(nct,17)= drhsc*xpt*ypt*zpt wrk(nct,18)= drhsn*xpt*ypt*zpt wrk(nct,19)= drhsc*xmt*ypt*zpt wrk(nct,20)= drhsn*xmt*ypt*zpt wrk(nct,21)= drhsc*xmt*ymt*zpt wrk(nct,22)= drhsn*xmt*ymt*zpt wrk(nct,23)= drhsc*xpt*ymt*zpt wrk(nct,24)= drhsn*xpt*ymt*zpt wrk(nct,25)= drhsc*xpt*ypt*zmt wrk(nct,26)= drhsn*xpt*ypt*zmt wrk(nct,27)= drhsc*xmt*ypt*zmt wrk(nct,28)= drhsn*xmt*ypt*zmt wrk(nct,29)= drhsc*xmt*ymt*zmt wrk(nct,30)= drhsn*xmt*ymt*zmt wrk(nct,31)= drhsc*xpt*ymt*zmt wrk(nct,32)= drhsn*xpt*ymt*zmt 820 continue m = (l-1)*2*na i1 = m+2*nd1-1 i2 = i1+1 i3 = m+2*nd2-1 i4 = i3+1 i5 = m+2*nd3-1 i6 = i5+1 i7 = m+2*nd4-1 i8 = i7+1 i9 = m+2*nd5-1 i10 = i9+1 i11 = m+2*nd6-1 i12 = i11+1 i13 = m+2*nd7-1 i14 = i13+1 i15 = m+2*nd8-1 i16 = i15+1 dlhsc = -mu*army1*(cwgt1/c1-hwgt*dhca/h1) dlhsn = mu*army1*hwgt*dhna/h1 work(i1) = sum(3)+dlhsc*xp*yp*zp work(i2) = sum(4)+dlhsn*xp*yp*zp work(i3) = sum(5)+dlhsc*xm*yp*zp work(i4) = sum(6)+dlhsn*xm*yp*zp work(i5) = sum(7)+dlhsc*xm*ym*zp work(i6) = sum(8)+dlhsn*xm*ym*zp work(i7) = sum(9)+dlhsc*xp*ym*zp work(i8) = sum(10)+dlhsn*xp*ym*zp work(i9) = sum(11)+dlhsc*xp*yp*zm work(i10) = sum(12)+dlhsn*xp*yp*zm work(i11) = sum(13)+dlhsc*xm*yp*zm work(i12) = sum(14)+dlhsn*xm*yp*zm work(i13) = sum(15)+dlhsc*xm*ym*zm work(i14) = sum(16)+dlhsn*xm*ym*zm work(i15) = sum(17)+dlhsc*xp*ym*zm work(i16) = sum(18)+dlhsn*xp*ym*zm tem1 = c1-c0-z(k)*lam*(x(i)/n0)**theta*thet1*h+ c 2.d0*alpha*(n0-y(j))+dpn dlhsc =-army1*(mu-gamm1/gamma*mu*c1/h1*dhca- c gamma*c0**cwgt1*dc0ca-tem1*(cwgt1*mu/c1- c hwgt*mu/h1*dhca)-mu*(1.d0-dc0ca-z(k)* c lam*(x(i)/n0)**theta*thet1*dhca)) dlhsn =-army1*(-gamm1/gamma*mu*c1/h1*dhna-gamma* c c0**cwgt1*dc0na+tem1*hwgt*mu/h1* c dhna-mu*(d2pn-dc0na+z(k)*lam*(x(i)/n0)** c theta*thet1*h*(theta/n0-dhna/h)+2.d0* c alpha)) wrk2(i1) = sum(19)+dlhsc*xp*yp*zp wrk2(i2) = sum(20)+dlhsn*xp*yp*zp wrk2(i3) = sum(21)+dlhsc*xm*yp*zp wrk2(i4) = sum(22)+dlhsn*xm*yp*zp wrk2(i5) = sum(23)+dlhsc*xm*ym*zp wrk2(i6) = sum(24)+dlhsn*xm*ym*zp wrk2(i7) = sum(25)+dlhsc*xp*ym*zp wrk2(i8) = sum(26)+dlhsn*xp*ym*zp wrk2(i9) = sum(27)+dlhsc*xp*yp*zm wrk2(i10) = sum(28)+dlhsn*xp*yp*zm wrk2(i11) = sum(29)+dlhsc*xm*yp*zm wrk2(i12) = sum(30)+dlhsn*xm*yp*zm wrk2(i13) = sum(31)+dlhsc*xm*ym*zm wrk2(i14) = sum(32)+dlhsn*xm*ym*zm wrk2(i15) = sum(33)+dlhsc*xp*ym*zm wrk2(i16) = sum(34)+dlhsn*xp*ym*zm CFPP$ UNROLL do 830 m=1,ns CFPP$ UNROLL do 830 mm=1,ne1 nct = (m-1)*ne1+mm j1 = (m-1)*2*na+2*ind(nct)-1 j2 = j1+1 j3 = j1+2 j4 = j3+1 j5 = j3+2*nx1 j6 = j5+1 j7 = j1+2*nx1 j8 = j7+1 j9 = j1+2*nx1*ny1 j10 = j9+1 j11 = j3+2*nx1*ny1 j12 = j11+1 j13 = j11+2*nx1 j14 = j13+1 j15 = j9+2*nx1 j16 = j15+1 work(j1) = work(j1) + wrk(nct,1) work(j2) = work(j2) + wrk(nct,2) work(j3) = work(j3) + wrk(nct,3) work(j4) = work(j4) + wrk(nct,4) work(j5) = work(j5) + wrk(nct,5) work(j6) = work(j6) + wrk(nct,6) work(j7) = work(j7) + wrk(nct,7) work(j8) = work(j8) + wrk(nct,8) work(j9) = work(j9) + wrk(nct,9) work(j10) = work(j10) + wrk(nct,10) work(j11) = work(j11) + wrk(nct,11) work(j12) = work(j12) + wrk(nct,12) work(j13) = work(j13) + wrk(nct,13) work(j14) = work(j14) + wrk(nct,14) work(j15) = work(j15) + wrk(nct,15) work(j16) = work(j16) + wrk(nct,16) wrk2(j1) = wrk2(j1) + wrk(nct,17) wrk2(j2) = wrk2(j2) + wrk(nct,18) wrk2(j3) = wrk2(j3) + wrk(nct,19) wrk2(j4) = wrk2(j4) + wrk(nct,20) wrk2(j5) = wrk2(j5) + wrk(nct,21) wrk2(j6) = wrk2(j6) + wrk(nct,22) wrk2(j7) = wrk2(j7) + wrk(nct,23) wrk2(j8) = wrk2(j8) + wrk(nct,24) wrk2(j9) = wrk2(j9) + wrk(nct,25) wrk2(j10) = wrk2(j10) + wrk(nct,26) wrk2(j11) = wrk2(j11) + wrk(nct,27) wrk2(j12) = wrk2(j12) + wrk(nct,28) wrk2(j13) = wrk2(j13) + wrk(nct,29) wrk2(j14) = wrk2(j14) + wrk(nct,30) wrk2(j15) = wrk2(j15) + wrk(nct,31) wrk2(j16) = wrk2(j16) + wrk(nct,32) 830 continue res1 = (-mu*army1+sum(1))*u(i)*v(j)*w(k) res2 = ((-c1**cwgt*h1**hwgt/ome1+c0**cwgt/ c ome1+mu*(dpn+c1-c0-z(k)*lam*(x(i)/ c n0)**theta*thet1*h+2.d0*alpha*(n0- c y(j))))*army1+sum(2))*u(i)*v(j)*w(k) f(i1) = f(i1) + res1*xp*yp*zp f(i2) = f(i2) + res2*xp*yp*zp f(i3) = f(i3) + res1*xm*yp*zp f(i4) = f(i4) + res2*xm*yp*zp f(i5) = f(i5) + res1*xm*ym*zp f(i6) = f(i6) + res2*xm*ym*zp f(i7) = f(i7) + res1*xp*ym*zp f(i8) = f(i8) + res2*xp*ym*zp f(i9) = f(i9) + res1*xp*yp*zm f(i10) = f(i10) + res2*xp*yp*zm f(i11) = f(i11) + res1*xm*yp*zm f(i12) = f(i12) + res2*xm*yp*zm f(i13) = f(i13) + res1*xm*ym*zm f(i14) = f(i14) + res2*xm*ym*zm f(i15) = f(i15) + res1*xp*ym*zm f(i16) = f(i16) + res2*xp*ym*zm c tem3 = u(i)*v(j)*w(k) do 840 mm=1,nas tem1 = tem3*work(mm) tem2 = tem3*wrk2(mm) df(i1,mm) = df(i1,mm) + tem1*xp*yp*zp df(i2,mm) = df(i2,mm) + tem2*xp*yp*zp df(i3,mm) = df(i3,mm) + tem1*xm*yp*zp df(i4,mm) = df(i4,mm) + tem2*xm*yp*zp df(i5,mm) = df(i5,mm) + tem1*xm*ym*zp df(i6,mm) = df(i6,mm) + tem2*xm*ym*zp df(i7,mm) = df(i7,mm) + tem1*xp*ym*zp df(i8,mm) = df(i8,mm) + tem2*xp*ym*zp df(i9,mm) = df(i9,mm) + tem1*xp*yp*zm df(i10,mm) = df(i10,mm) + tem2*xp*yp*zm df(i11,mm) = df(i11,mm) + tem1*xm*yp*zm df(i12,mm) = df(i12,mm) + tem2*xm*yp*zm df(i13,mm) = df(i13,mm) + tem1*xm*ym*zm df(i14,mm) = df(i14,mm) + tem2*xm*ym*zm df(i15,mm) = df(i15,mm) + tem1*xp*ym*zm df(i16,mm) = df(i16,mm) + tem2*xp*ym*zm 840 continue 750 continue 720 continue 710 continue call dgeco(df,nas,nas,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when solving df*x=f' call dgesl(df,nas,nas,iwork,f,0) sum3 = 0.d0 do 850 j=1,ns do 850 i=1,na i1 = (j-1)*na*2+2*i-1 tem1 = f(i1) tem2 = f(i1+1) a1(i,j) = a0(i,j) - tem1 b1(i,j) = b0(i,j) - tem2 sum3 = sum3 + tem1*tem1+tem2*tem2 850 continue sumres = dsqrt(sum3)/float(nas) write(*,*) ' at iteration ', it, ' the residual is ', sumres do 860 j=1,ns do 860 i=1,na a0(i,j) = a1(i,j) b0(i,j) = b1(i,j) 860 continue if (sumres.lt.crit) go to 999 690 continue 999 continue c c print results c 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) gamma write(10,801) lam write(10,801) omega write(10,801) pa write(10,801) pb write(10,801) pc write(10,801) psi write(10,801) sig write(10,801) theta write(10,801) zeta write(10,812) itype write(10,812) ns do 870 i=1,ns write(10,801) ginv(i) write(10,801) cg(i) write(10,801) army(i) 870 continue do 880 j=1,ns do 880 i=1,ns write(10,801) mcp(i,j) 880 continue write(10,812) nx1 write(10,812) ny1 write(10,812) nz1 do 890 i=1,nx1 write(10,801) xa(i) 890 continue do 900 i=1,ny1 write(10,801) ya(i) 900 continue do 910 i=1,nz1 write(10,801) za(i) 910 continue do 920 j=1,ns do 920 i=1,na write(10,801) a0(i,j) 920 continue do 930 j=1,ns do 930 i=1,na write(10,801) b0(i,j) 930 continue c 801 format (4X, E12.6) 802 format (2(4X, F12.5)) 803 format (3(4X, F12.5)) 804 format (4(4X, F12.5)) 805 format (5(4X, F12.5)) 806 format (6(2X, F11.5)) 807 format (7(2X, F10.5)) 808 format (8(2X, F9.5)) 811 format (4X, E12.6, 4X, I4) 812 format (4X, I4) 813 format (4X, 3(1X, E9.3)) 1001 format (A27,4(F8.3,1X)) 1002 format (A27,E10.5,1X,2(F8.3,1X)) stop end