subroutine opacder(d,t,opcs,dlklr,dlklt) * * opacities, but from Lagrange re-interpolator * Call for each chemical species. If pure, use derivatives from * within the routine. Otherwise, use 5 point differencing. * See opac3 for further details. * implicit double precision(a-h,o-z) common/je/je common/comp/amhyhe,amheca,x(4) common/opcswch/first logical first real*8 op(4),dkrx(4),ytb(4),xtb(4) real*8 ta(30),tb(30),ty(30,30,3) * * read in the opacities if we haven't already * if ( first ) then first = .false. read(25,*) iaend,jbend read(25,705) (ta(li),li=1,iaend) read(25,705) (tb(li),li=1,jbend) 705 format(8f10.5) 710 format(10f8.4) do 1000 lk=1,3 do 1000 li=1,iaend read(25,710) (ty(li,lj,lk),lj=1,jbend) 1000 continue endif * * locate our place on the grid. * ider=1 if(km.lt.1)then i = 9 j = 7 km = 0 endif 65 continue km=km+1 * * set temp = xa and dens = xb * xa = t xb = d if( xa .lt. ta(1) ) go to 50 if( xa .gt. ta(iaend)) go to 50 1 if ( xa - ta(i)) 4, 6, 2 2 if (xa - ta(i+1)) 6, 5, 3 3 i = i + 2 go to 1 4 i = i - 1 go to 1 5 i = i + 1 6 continue if( xb .lt. tb(1) ) go to 50 if( xb .gt. tb(jbend)) go to 50 11 if (xb - tb(j)) 14, 16, 12 12 if(xb-tb(j+1)) 16, 15, 13 13 j = j + 2 go to 11 14 j = j - 1 go to 11 15 j = j + 1 16 continue itop=0 it10=i-1 if(it10 .ge. 1) go to 150 it10=it10+1 itop=1 150 continue if(ty(it10,j,je).gt.-7.99 .or. ty(it10,j+1,je).gt.-7.99) goto 160 it10=it10+1 itop=1 if(it10 .gt.i) go to 400 go to 150 160 it1f=it10+3 if(it1f.le.iaend) goto 170 if(itop .ne. 0) go to 400 it10=it10-1 goto 150 170 continue if(ty(it1f,j,je).gt.-7.99 .or. ty(it1f,j+1,je).gt.-7.99)goto 180 if(itop .ne. 0) go to 400 it10 = it10-1 if(it10.le. 0) go to 400 if(it1f .le. i+1) go to 400 go to 150 180 continue do 250 int=1,4 it1=it10+int-1 if(ty(it1,j,je).le.-7.99 .or. ty(it1,j+1,je).le.-7.99) go to 200 js1=j-1 if(js1 .le. 0) js1=js1+1 if(ty(it1,js1,je) .le. -7.99) js1=js1+1 if(js1+3 .gt. jbend) js1=js1-1 if(ty(it1,js1+3,je) .le. -7.99) js1=js1-1 do 190 k=1,4 k1=js1+k-1 * * add conductive opac. to rad. opac. * ocond = ocon(tb(k1),ta(it1)) ytb(k)=dlog10(1./(10.**(-ty(it1,k1,je))+10.**(-ocond))) 190 xtb(k)=tb(k1) call fh4(xb,ytb,xtb,op(int),dkrx(int),ider) go to 250 200 if(ty(it1,j,je) .le. -7.99) js1=j+1 if(ty(it1,j+1,je) .le. -7.99) js1=j-1 * * add conductive opac. to rad. opac. * ocond1 = ocon(tb(js1+1),ta(it1)) ocond2 = ocon(tb(js1),ta(it1)) op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1))) op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2))) dkrx(int)=(op1-op2)/(tb(js1+1)-tb(js1)) op(int)=op1+dkrx(int)*(xb-tb(js1+1)) 250 continue do 260 k=1,4 260 xtb(k)=ta(it10+k-1) call fh4(xa,op,xtb,opc,dlklt,ider) opcs=opc if(ider .ne. 1)then return endif call fh4(xa,dkrx,xtb,dlklr,dumfh,0) return 400 continue do 500 k=1,2 it1=i+k-1 if(ty(it1,j,je).le.-7.99 .and. ty(it1,j+1,je).le.-7.99) go to 450 js1=j if(ty(it1,j,je).le.-7.99) js1=j+1 if(ty(it1,j+1,je).le.-7.99) js1=j-1 frxs=(xb-tb(js1+1))/(tb(js1+1)-tb(js1)) * * add conductive opac. to rad. opac. * ocond1 = ocon(tb(js1+1),ta(it1)) ocond2 = ocon(tb(js1),ta(it1)) op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1))) op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2))) op(k)=op1+(op1-op2)*frxs if(ider .ne. 1) go to 500 * * add conductive opac. to rad. opac. * dop1=(op1-op2)/(tb(js1+1)-tb(js1)) dop2=dop1 if(js1-1 .le. 0) go to 420 if(ty(it1,js1-1,je) .gt. -7.99) then * * add conductive opac. to rad. opac. * ocond1 = ocon(tb(js1-1),ta(it1)) ocond2 = ocon(tb(js1),ta(it1)) op1=dlog10(1./(10.**(-ty(it1,js1-1,je))+10.**(-ocond1))) op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2))) dop1=0.5*(dop1+(op2-op1)/(tb(js1)-tb(js1-1))) endif 420 if(js1+2 .gt. jbend) go to 430 if(ty(it1,js1+2,je).gt.-7.99) then * * add conductive opac. to rad. opac. * ocond1 = ocon(tb(js1+1),ta(it1)) ocond2 = ocon(tb(js1+2),ta(it1)) op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1))) op2=dlog10(1./(10.**(-ty(it1,js1+2,je))+10.**(-ocond2))) dop2=0.5*(dop2+(op3-op1)/(tb(js1+2)-tb(js1+1))) endif 430 dkrx(k)=dop2+(dop2-dop1)/(tb(js1+1)-tb(js1))*(xb-tb(js1+1)) go to 500 450 continue jc=1 jco=0 js1=j if(j.gt. 11) jc=-1 455 do 460 ljs=1,jbend js1=js1+jc if(js1.ge.jbend .or. js1.le.1) go to 465 if(ty(it1,js1,je) .gt. -7.99) go to 470 460 continue 465 if(jco .ne. 0) go to 50 jco=jco+1 jc=-jc go to 455 470 if(jc .eq. -1) js1=js1-1 * * add conductive opac. to rad. opac. * ocond1 = ocon(tb(js1+1),ta(it1)) ocond2 = ocon(tb(js1),ta(it1)) op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1))) op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2))) dkrx(k)=(op1-op2)/(tb(js1+1)-tb(js1)) op(k)=op1+dkrx(k)*(xb-tb(js1+1)) 500 continue dlklt=(op(2)-op(1))/(ta(i+1)-ta(i)) opcs=op(2)+dlklt*(xa-ta(i+1)) if(ider .eq. 0) return dlklr=dkrx(2)+(dkrx(2)-dkrx(1))*(xa-ta(i+1))/(ta(i+1)-ta(i)) return 50 continue stop end ************************************************************************