program crystal ! !*********************************************************************** ! !! CRYSTAL is the main program for the crystal simulation. ! ! ! Modified: ! ! 22 March 2002 ! ! Author: ! ! John Burkardt ! implicit none ! integer, parameter :: liv = 60 integer, parameter :: npar = 1 ! integer, parameter :: lv = 77+npar*(npar+17)/2 ! real cost external cryfun real d(npar) external dummy integer i integer ido integer iopt integer ipar(1) integer iv(liv) integer nfun real par(npar) real rpar(1) real tarray(2) real temp real time1 real time2 real v(lv) ! call timestamp ( ) call etime ( tarray ) time1 = tarray(1) + tarray(2) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL:' write ( *, '(a)' ) ' Version of 25 January 1996' iopt = 0 par(1) = 1713.0E+00 ! ! If IOPT = 0, just call CRYFUN once. ! if ( iopt == 0 ) then call cryfun ( npar, par, nfun, cost, ipar, rpar, dummy ) ! ! If IOPT = 1, then carry out an optimization. ! else d(1:npar) = 1.0E+00 ! ! Set the 611 work arrays to default values. ! ido = 2 call deflt ( ido, iv, liv, lv, v ) ! ! Tell SMSNO that we have called DEFLT already. ! iv(1) = 12 ! ! Set the maximum number of iterations. ! iv(18) = 40 call smsno ( npar, d, par, cryfun, iv, liv, lv, v, ipar, rpar, dummy ) end if call etime ( tarray ) time2 = tarray(1) + tarray(2) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL:' write ( *, '(a,g14.6,a)' )' Total elapsed CPU time = ', & time2 - time1 ,' seconds,' write ( *, '(a,g14.6,a)' )' = ', & ( time2 - time1 ) / 60.0E+00, 'minutes.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine dummy ! !*********************************************************************** ! !! DUMMY is a dummy subroutine needed as formal input to CRYFUN. ! ! ! Modified: ! ! 22 March 2002 ! ! Author: ! ! John Burkardt ! implicit none ! return end subroutine cryfun ( npar, par, nfun, cost, ipar, rpar, dummy ) ! !*********************************************************************** ! !! CRYFUN evaluates the cost function for the minimizing software. ! ! implicit none ! integer, parameter :: maxbot = 20 integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer npar integer, parameter :: ns = 10 ! real ae1(ni,nj) real ae2(ni,nj) real ak1(ni,nj) real ak2(ni,nj) real area(ni,nj) real areal real areas real areat real b real b1jbl(nj) real b2jbl(nj) real birad real bo real cappa real cd real ce1 real ce2 real cfo real cinc real cinco real cmu real cost real cvn real delt real dtm external dummy real epsad real epsil real ewall real f(ni,nj,ns) real fcsl real fjeta(ni,nj) real fjksi(ni,nj) real fks real fksl real fma real fmax(ns) real fn(ni,nj,ns) real fnu real fo(ni,nj,ns) real fr real frsl real gamt(ni,nj) real grash real hamag real heta(ni,nj) real hf real hksi(ni,nj) integer i integer icost integer icrys integer inturb integer ipar(1) integer iplot integer ipref integer iprint integer iter integer izone integer j integer jcrys integer jpref integer k integer l0 integer l1 integer last integer lastt logical lblk(ns) logical lconv logical lortho logical lsolve(ns) integer m0 integer m1 integer mode integer nbot integer, save :: ncall = 0 integer ndt integer nf integer nfun integer np integer npc integer nsolve(ns) integer ntimes(ns) real orth real par(npar) real pr real r(ni,nj) real ra real rdtm real re real recb real rect real relax(nk) real res(ns) real rho(ni,nj) real rhocon real rpar(1) real rpr real rueta(ni,nj) real ruksi(ni,nj) real sige real sigk real sigma real sigt real smax real smooth real ssum real stel real stes real t1 real t2 real t3 real tal real tanca real tanca2 real tas real tend real tf real tinit character ( len = 25 ) title(ns) real tnow real tw real vave real vol(ni,nj) real x(ni,nj) real xbot(maxbot) real xc(ni,nj) real xlen real y(ni,nj) real ybot(maxbot) real yc(ni,nj) real ylen ! ncall = ncall+1 write ( *, * ) 'CRYFUN: PAR(1) = ',par(1) ! ! Initialize the data. ! call inidat(ae1,ae2,ak1,ak2,area,b,birad,bo,cappa,cd,ce1, & ce2,cfo,cinc,cmu,cost,cvn,delt,dtm,epsad,epsil,ewall,f,fcsl,fks, & fksl,fma,fn,fnu,fo,fr,frsl,gamt,grash,hamag,heta,hf, & hksi,icost,icrys,inturb,iplot,ipref,iprint,jcrys, & jpref,l0,l1,last,lastt,lblk,lortho,lsolve,m0,m1,maxbot,mode, & nbot,ndt,ni,nj,nk,np,npar,npc,ns,nsolve,ntimes,orth, & par,pr,ra,rdtm,re,recb,rect,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk, & sigma,sigt,smooth,stel,stes,tal,tanca,tanca2,tas,tend,tf, & tinit,title,tnow,tw,vave,vol,x,xbot,xc,xlen,y,ybot,yc,ylen) ! ! Print out data. ! if ( ncall == 1 ) then call prdat(b,birad,bo,cappa,cfo,cvn,delt,dtm,fcsl,fksl,fma, & fnu,fr,frsl,grash,icost,icrys,inturb,iprint, & jcrys,l0,last,lastt,m0,mode,nbot,ns,nsolve,ntimes,orth, & pr,ra,rdtm,recb,rect,rhocon,smooth,stel,stes,tanca,tanca2, & tend,tf,tinit,title,tw,vave,xlen,ylen) end if ! ! Generate the initial grid XC, YC. ! call inigrd(icrys,jcrys,l0,m0,nbot,xbot,xc,xlen,ybot,yc,ylen) ! ! Adapt the grid. ! do k = 1, 10 call adapt ( cvn, epsad, icrys, iprint, jcrys, l1, m1, nbot, & orth, smooth, xbot, xc, ybot, yc ) end do ! ! Once XC and YC are determined, compute the control volume areas. ! call doarea ( area, areal, areas, areat, icrys, jcrys, l0, m0, & ni, nj, xc, yc ) ! ! Each time iteration begins at this point. ! 10 continue ! ! Begin the iterative solution of the state equations in the ! solid zone. ! izone = 1 l1 = l0 m1 = jcrys ! ! Only the temperature (variable 5) needs to be solved for. ! lsolve(1) = .false. lsolve(2) = .false. lsolve(3) = .false. lsolve(4) = .false. lsolve(5) = .true. lsolve(6) = .false. lsolve(7) = .false. lsolve(8) = .false. lsolve(9) = .false. lsolve(10) = .false. do iter = 0, last lconv = .true. ! ! Set the X and Y coordinates of the primary nodes from XC, YC, ! the coordinates of the corner nodes. ! ! Note that this is only done for the left portion of the region! ! call setx ( l1, m1, ni, nj, x, xc, y, yc ) ! ! Compute various geometric quantities required for computing ! derivatives. ! call setgeo(ae1,ae2,ak1,ak2,heta,hksi,l1,m1,mode,ni,nj, & r,vol,x,xc,y,yc) ! ! Estimate pressure and momentum at control volume interfaces. ! if ( ndt == 0 ) then if ( iter == 0 ) then call initl(fjeta,fjksi,heta,hksi,l1,m1,ni,nj,f(1,1,3),rho, & rueta,ruksi,f(1,1,1),f(1,1,2),x,y) end if end if call setup(ae1,ae2,ak1,ak2,b1jbl, & b2jbl,birad,cappa,cd,ce1,ce2,cmu,epsil,ewall,f, & fcsl,fjeta,fjksi,fksl,fma,fmax,fn,fo,frsl,gamt,grash,hamag, & heta,hksi,inturb,icrys,iter,izone,jcrys,l1,lblk,lconv,lortho, & lsolve,m1,mode,nf,np,npc,nsolve,ntimes,pr,r,rdtm,re,recb, & rect,res,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk,sigt, & smax,ssum,stel,stes,tal,tas,tf,tw,vol,x,xc,y,yc) if ( iprint > 0 ) then call output(cfo,iter,izone,res,smax,ssum,f(1,1,5), & tnow,f(1,1,1),f(1,1,6)) end if if ( lconv .and. iter >= 5 ) go to 20 end do if ( iprint > 0 ) then write ( *, * ) ' ' write ( *, * ) 'CRYFUN - Warning!' write ( *, * ) ' The solid iteration has not converged' write ( *, * ) ' after ',iter,' iterations.' write ( *, * ) ' ' write ( *, * ) ' Solid norm and max relative change:' write ( *, * ) ' ' do i = 1,ns if ( lsolve(i) ) then fmax(i) = maxval ( f(1:l0,1:m0,i) ) write(*,'(i3,2x,a25,2g14.6)') i, title(i), fmax(i), res(i) end if end do end if ! ! Zone 2 (LIQUID) calculation ! 20 continue do i = 1,icrys do j = 1,jcrys f(i,j,5) = fo(i,j,5) f(i,j,9) = fo(i,j,9) end do end do izone = 2 l1 = icrys m1 = m0 rueta(1:icrys,2) = 0.0E+00 rueta(1:icrys,m0) = 0.0E+00 rueta(1,1:m0) = 0.0E+00 rueta(icrys,1:m0) = 0.0E+00 ruksi(1:icrys,1) = 0.0E+00 ruksi(1:icrys,m0) = 0.0E+00 ruksi(2,1:m0) = 0.0E+00 ruksi(icrys,1:m0) = 0.0E+00 if ( inturb == 0 ) then lsolve(1) = .true. lsolve(2) = .true. lsolve(3) = .true. lsolve(4) = .true. lsolve(5) = .true. lsolve(6) = .false. lsolve(7) = .false. lsolve(8) = .false. lsolve(9) = .false. lsolve(10) = .false. else lsolve(1) = .true. lsolve(2) = .true. lsolve(3) = .true. lsolve(4) = .true. lsolve(5) = .true. lsolve(6) = .false. lsolve(7) = .true. lsolve(8) = .true. lsolve(9) = .false. lsolve(10) = .false. end if do iter = 1,last lconv = .true. ! ! Set the turbulent viscosity. ! if ( inturb == 0 ) then do i = 2,l1-1 do j = 2,m1-1 gamt(i,j) = 0.0E+00 end do end do else if ( inturb == 1 ) then do i = 2,l1-1 do j = 2,m1-1 gamt(i,j) = (1.0-relax(12))*gamt(i,j) & +relax(12)*cmu*rho(i,j)*f(i,j,7)**2/f(i,j,8) end do end do end if ! ! Set the X, Y values. ! call setx(l1,m1,ni,nj,x,xc,y,yc) ! ! Compute various geometric quantities required for computing ! derivatives. ! call setgeo(ae1,ae2,ak1,ak2,heta,hksi,l1,m1,mode,ni,nj, & r,vol,x,xc,y,yc) ! ! Set the right hand side of certain flux boundary conditions. ! call setup(ae1,ae2,ak1,ak2,b1jbl, & b2jbl,birad,cappa,cd,ce1,ce2,cmu,epsil,ewall,f, & fcsl,fjeta,fjksi,fksl,fma,fmax,fn,fo,frsl,gamt,grash,hamag, & heta,hksi,inturb,icrys,iter,izone,jcrys,l1,lblk,lconv, & lortho, & lsolve,m1,mode,nf,np,npc,nsolve,ntimes,pr,r,rdtm,re,recb, & rect,res,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk,sigt, & smax,ssum,stel,stes,tal,tas,tf,tw,vol,x,xc,y,yc) ! ! Compute the stream function PSI. ! f(2,2,10) = 0.0 do i = 2, icrys if ( i > 2 ) then f(i,2,10) = f(i-1,2,10) - rueta(i-1,2) * ae1(i-1,2) & - 0.5E+00 * ( ruksi(i-1,1) + ruksi(i,1) ) * ae2(i-1,2) end if do j = 3, m1-1 t1 = (rueta(i,j-1)+rueta(i,j)) t2 = (rueta(i-1,j-1)+rueta(i-1,j)) t3 = 0.5*(hksi(i,j-1)*t2+hksi(i-1,j-1)*t1)/(hksi(i,j-1) & +hksi(i-1,j-1)) f(i,j,10) = f(i,j-1,10)+ruksi(i,j-1)*ak1(i,j-1)-t3*ak2(i,j-1) end do end do f(1,1,10) = f(2,2,10) do j = 2, m0 f(1,j,10) = f(2,j,10) end do do i = 2, l0 f(i,1,10) = f(i,2,10) end do ! ! Optional printout. ! if ( iprint > 0 ) then call output(cfo,iter,izone,res,smax,ssum,f(1,1,5), & tnow,f(1,1,1),f(1,1,6)) end if if ( lconv ) go to 30 end do if ( iprint > 0 ) then write ( *, * ) ' ' write ( *, * ) 'CRYFUN - Warning!' write ( *, * ) ' The liquid iteration has not converged' write ( *, * ) ' after ',iter,' iterations.' write ( *, * ) ' ' write ( *, * ) ' Liquid norms and max relative change:' write ( *, * ) ' ' do i = 1, ns if ( lsolve(i) ) then fmax(i) = maxval ( f(1:l0,1:m0,i) ) write(*,'(i3,2x,a25,2g14.6)') i, title(i), fmax(i), res(i) end if end do end if ! ! We have computed the state variables for the current time. ! 30 continue ! ! Evaluate the cost function integrand at the current time. ! cinco = cinc call setcst ( area, cinc, icost, icrys, m0, ni, nj, f(1,1,5), tf, & f(1,1,1), f(1,1,2), vave ) ! ! Update the total cost, by adding the estimated contribution ! from the current time interval. ! ! TEMPORARY ! if ( ndt > 0 ) then ! cost = cost + 0.5 * delt *( cinco + cinc ) cost = cost + 0.5 * dtm * ( cinco + cinc ) end if ! ! If we've reached the end time, then write out final data, ! possibly save a restart file, and stop. ! if ( ndt >= lastt ) then if ( iprint > 0 ) then call pmod(ipref,jcrys,jpref,l1,m1,f(1,1,3)) end if cost = cost / ( sqrt ( areal ) * ( tend - tinit ) ) ! ! If requested, write out plot data. ! Sadly, it is necessary to call SETX to generate the X, Y arrays ! for the entire region. Previous calls only generate them for ! a portion of the region. ! if ( iplot == 1 ) then call setx ( l0, m0, ni, nj, x, xc, y, yc ) call rswrit ( cost, f(1,1,9), gamt, icrys, jcrys, l0, m0, nbot, & f(1,1,3), f(1,1,4), f(1,1,10), & rueta, ruksi, f(1,1,5), f(1,1,8), f(1,1,7), tnow, f(1,1,1), & f(1,1,2), f(1,1,6), x, xbot, xc, y, yc, ybot ) end if write ( *, * ) 'CRYFUN: COST = ',cost return end if ! ! Update the number of steps, and the current time. ! ndt = ndt+1 tnow = tinit+ndt*dtm ! ! Save a copy of the current data. ! do i = 1,l0 do j = 1,m0 do k = 1,ns fo(i,j,k) = f(i,j,k) end do end do end do call movgrd(b1jbl,b2jbl,bo,delt,fksl,fr,frsl,icrys,iprint,jcrys, & l0,m0,f(1,1,3),pr,re,stel,tanca,tanca2,xc,xlen,yc) l1 = l0 m1 = m0 call adapt(cvn,epsad,icrys,iprint,jcrys,l1,m1,nbot, & orth,smooth,xbot,xc,ybot,yc) ! ! Once XC and YC are determined, compute the control volume areas. ! call doarea(area,areal,areas,areat,icrys,jcrys,l0,m0,ni,nj,xc,yc) go to 10 end subroutine adapt ( cvn, epsad, icrys, iprint, jcrys, l1, m1, nbot, & orth, smooth, xbot, xc, ybot, yc ) ! !*********************************************************************** ! !! ADAPT executes MAGG, the multizone adaptive grid generation. ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! real a1 real a11 real a12 real a2 real a21 real a22 real a3 real b1 real b2 real b3 real c1 real c2 real c3 real cvn real det real dot real dx real dxc(ni,nj) real dxmax real dy real dyc(ni,nj) real epsad real gn real gt real gx real gy integer i integer icrys integer iprint integer j integer jcrys integer l1 integer m1 integer nin real orth real pb(ni,nj) real pc1 real pc2 real pd(ni,nj) real ratio real res1 real res2 real rmax real s1 real s2 real scale real smooth real ssmot real ssweg real vmag2 real wn real wt real wx real wy real xbot(nbot) real xc(ni,nj) real xdisp real xij real xjac real xn real xnn real xpp real xt real xtn real xtt real ybot(nbot) real yc(ni,nj) real ydisp real yij real yn real ynn real ypp real yt real ytn real ytt ! dxmax = 0.0E+00 dxc(1:l1,1:m1) = 0.0E+00 dyc(1:l1,1:m1) = 0.0E+00 do i = 2,l1 do j = 2,m1 if ( i < icrys ) then pc1 = 0.95 * ((i-(2+icrys)/2)/10.0)**2 + 0.05E+00 else pc1 = 0.8 * ((i-l1)/12.0)**2 + 0.2E+00 end if if ( j <= jcrys ) then pc2 = 0.9 * ((j-(2+jcrys)/2)/10.0)**2 + 0.05 else pc2 = 0.9 * ((j-(m1+jcrys)/2)/10.0)**2 + 0.05 end if pb(i,j) = pc1 * pc2 pd(i,j) = pc1 * pc2 end do end do ! ! This DO loop iteration is an iterative solution of a linear system. ! do nin = 1,100 ssmot = 0.0E+00 ssweg = 0.0E+00 rmax = 0.0E+00 if ( iprint > 0 ) then if ( nin == 1 ) then write ( *, * ) ' ' write ( *, * ) 'ADAPT: Iter DXMAX XC(20,20) YC(20,20)' // & ' XC(3,M1) YC(3,M1) RMAX' end if if ( mod(nin,50) == 0 ) then write(*,'(6x,i5,6g10.3)') & nin,dxmax,xc(20,20),yc(20,20),xc(3,m1),yc(3,m1),rmax end if end if do j = 3, m1-1 do i = 3, l1-1 xt = 0.5 * ( xc(i+1,j) - xc(i-1,j) ) xn = 0.5 * ( xc(i,j+1) - xc(i,j-1) ) yt = 0.5 * ( yc(i+1,j) - yc(i-1,j) ) yn = 0.5 * ( yc(i,j+1) - yc(i,j-1) ) xjac = xt*yn-xn*yt ! ! Refer to table 2.1, page 14, of Hui Zhang's thesis, for ! a description of these coefficients. ! a1 = -2.0*smooth*(xt*yt+xn*yn)*(xn*xn+yn*yn)/xjac**3 a2 = 4.0*smooth*(xt*yt+xn*yn)*(xt*xn+yt*yn)/xjac**3 a3 = -2.0*smooth*(xt*yt+xn*yn)*(xt*xt+yt*yt)/xjac**3 b1 = 2.0*smooth*(yt*yt+yn*yn)*(xn*xn+yn*yn)/xjac**3 b2 = -4.0*smooth*(yt*yt+yn*yn)*(xt*xn+yt*yn)/xjac**3 b3 = 2.0*smooth*(yt*yt+yn*yn)*(xt*xt+yt*yt)/xjac**3 c1 = 2.0*smooth*(xt*xt+xn*xn)*(xn*xn+yn*yn)/xjac**3 c2 = -4.0*smooth*(xt*xt+xn*xn)*(xt*xn+yt*yn)/xjac**3 c3 = 2.0*smooth*(xt*xt+xn*xn)*(xt*xt+yt*yt)/xjac**3 a1 = a1+orth*pd(i,j)*xn*yn a2 = a2+orth*pd(i,j)*(xt*yn+xn*yt) a3 = a3+orth*pd(i,j)*xt*yt b1 = b1+orth*pd(i,j)*xn*xn b2 = b2+2.0*orth*pd(i,j)*(2.0*xt*xn+yt*yn) b3 = b3+orth*pd(i,j)*xt*xt c1 = c1+orth*pd(i,j)*yn*yn c2 = c2+2.0*orth*pd(i,j)*(xt*xn+2.0*yt*yn) c3 = c3+orth*pd(i,j)*yt*yt a1 = a1-2.0*cvn*pb(i,j)*xn*yn a2 = a2+2.0*cvn*pb(i,j)*(xt*yn+xn*yt) a3 = a3-2.0*cvn*pb(i,j)*xt*yt b1 = b1+2.0*cvn*pb(i,j)*yn*yn b2 = b2-4.0*cvn*pb(i,j)*yt*yn b3 = b3+2.0*cvn*pb(i,j)*yt*yt c1 = c1+2.0*cvn*pb(i,j)*xn*xn c2 = c2-4.0*cvn*pb(i,j)*xt*xn c3 = c3+2.0*cvn*pb(i,j)*xt*xt ! ! Now compute terms from the derivatives of the weight functions. ! gt = 0.5*(pd(i+1,j)-pd(i-1,j)) gn = 0.5*(pd(i,j+1)-pd(i,j-1)) gx = orth*0.5*(gt*yn-gn*yt)*(xt*xn+yt*yn)**2/xjac gy = orth*0.5*(gn*xt-gt*xn)*(xt*xn+yt*yn)**2/xjac wt = 0.5*(pb(i+1,j)-pb(i-1,j)) wn = 0.5*(pb(i,j+1)-pb(i,j-1)) wx = cvn*xjac*(yn*wt-yt*wn) wy = cvn*xjac*(xt*wn-xn*wt) ! ! (RES1, RES2) is the residual that should be driven to zero. ! xtt = xc(i+1,j)-2.0*xc(i,j)+xc(i-1,j) xnn = xc(i,j+1)-2.0*xc(i,j)+xc(i,j-1) ytt = yc(i+1,j)-2.0*yc(i,j)+yc(i-1,j) ynn = yc(i,j+1)-2.0*yc(i,j)+yc(i,j-1) xtn = 0.25*(xc(i+1,j+1)+xc(i-1,j-1)-xc(i-1,j+1)-xc(i+1,j-1)) ytn = 0.25*(yc(i+1,j+1)+yc(i-1,j-1)-yc(i-1,j+1)-yc(i+1,j-1)) res1 = b1*xtt+b2*xtn+b3*xnn+a1*ytt+a2*ytn+a3*ynn+gx+wx res2 = a1*xtt+a2*xtn+a3*xnn+c1*ytt+c2*ytn+c3*ynn+gy+wy if ( abs(res1)+abs(res2) > rmax ) then rmax = abs(res1)+abs(res2) end if a11 = -2.0*(b1+b3) a12 = -2.0*(a1+a3) a22 = -2.0*(c1+c3) a21 = a12 det = a11*a22-a12*a21 ! ! (DX, DY) is the pointwise solution of the linear system. ! dx = (res2*a21-res1*a22)/det dy = (res1*a12-res2*a11)/det ! ! Now look at how movement of XC(I,J), YC(I,J) affects the neighbors ! ! (I+1,J-1) (I+1,J) (I+1,J+1) ! ! (I, J-1) (I, J) (I, J+1) ! ! (I-1,J-1) (I-1,J) (I-1,J+1) ! ! to determine the linear factor SCALE that will multiply (DX,DY). ! ratio = 0.25 xij = xc(i,j) yij = yc(i,j) xdisp = xc(i+1,j)-xij ydisp = yc(i+1,j)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i+1,j+1)-xij ydisp = yc(i+1,j+1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i,j+1)-xij ydisp = yc(i,j+1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i-1,j+1)-xij ydisp = yc(i-1,j+1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i-1,j)-xij ydisp = yc(i-1,j)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i-1,j-1)-xij ydisp = yc(i-1,j-1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i,j-1)-xij ydisp = yc(i,j-1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) xdisp = xc(i+1,j-1)-xij ydisp = yc(i+1,j-1)-yij vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) scale = min(0.5d0,0.25d0/ratio) dxc(i,j) = scale*dx dyc(i,j) = scale*dy ssmot = ssmot+(xt*xt+xn*xn+yt*yt+yn*yn)/xjac ssweg = ssweg+cvn*pb(i,j)*xjac**2 end do end do ! ! Move the points. ! do j = 3,m1-1 do i = 3,l1-1 if ( i /= icrys ) then xc(i,j) = xc(i,j)+dxc(i,j) end if if ( j /= jcrys ) then yc(i,j) = yc(i,j)+dyc(i,j) end if end do end do do i = icrys+1,l1-1 yc(i,jcrys) = 0.4+0.02*sin(xc(i,jcrys)*20.0*3.14159) end do ! ! Handle the nodes along the bottom row, from 10 positions to the right ! of the crystal, to the right hand wall. ! do j = jcrys+10,m1-1 ratio = 0.25 call findp(j,xc(3,j),yc(3,j),xpp,ypp,xc,yc) dx = xpp-xc(2,j) dy = ypp-yc(2,j) xdisp = xc(2,j+1)-xc(2,j) ydisp = yc(2,j+1)-yc(2,j) vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) s1 = vmag2 xdisp = xc(2,j-1)-xc(2,j) ydisp = yc(2,j-1)-yc(2,j) vmag2 = xdisp*xdisp+ydisp*ydisp dot = dx*xdisp+dy*ydisp ratio = max(dot/vmag2,ratio) s2 = max(vmag2,s1) scale = min(0.3d0,0.25d0/ratio) s1 = (scale*dx)**2+(scale*dy)**2 if ( s1 < s2 ) then xc(2,j) = xc(2,j)+scale*dx yc(2,j) = yc(2,j)+scale*dy end if call cubic(nbot,yc(2,j),ybot,xc(2,j),xbot) yc(l1,j) = yc(l1-1,j) end do ! ! Handle the nodes along the bottom row, from the left axis ! of symmetry, to 9 positions beyond the crystal. ! do j = 2,jcrys+9 yc(2,j) = yc(3,j) call cubic(nbot,yc(2,j),ybot,xc(2,j),xbot) yc(l1,j) = yc(l1-1,j) end do do i = 3,l1 xc(i,2) = xc(i,3) xc(i,m1) = xc(i,m1-1) if ( xc(i,m1-1) <= (xc(2,m1)+0.002*float(i-2)) ) then xc(i,m1) = xc(2,m1)+0.002*float(i-2) end if end do ! ! Compute the maximum movement of all corner nodes. ! dxmax = abs(dxc(3,3)) do j = 3,m1-1 do i = 3,l1-1 dxmax = max(dxmax,abs(dxc(i,j))) dxmax = max(dxmax,abs(dyc(i,j))) end do end do ! ! Copy values to the dummy corner nodes with I = 1 or J=1. ! xc(1,1) = xc(2,2) yc(1,1) = yc(2,2) do j = 2,m1 xc(1,j) = xc(2,j) yc(1,j) = yc(2,j) end do do i = 2,l1 xc(i,1) = xc(i,2) yc(i,1) = yc(i,2) end do ! ! Check for convergence. ! if ( dxmax <= epsad)return end do return end subroutine cubic ( nbot, ynew, ybot, xnew, xbot ) ! !*********************************************************************** ! !! CUBIC constructs a cubic spline through data. ! ! ! Discussion: ! ! CUBIC is given NBOT points (YBOT,XBOT) which lie along the curve ! defining the bottom of the crucible. ! ! CUBIC is also given a value YNEW which lies between YBOT(1) and ! YBOT(NBOT). ! ! CUBIC constructs a cubic spline through the data, and evaluates it at ! YNEW, returning the value as XNEW. ! ! Modified: ! ! 22 March 2002 ! ! Author: ! ! John Burkardt ! implicit none ! integer nbot integer, parameter :: ni = 64 ! real dome integer i integer ist real p2 real p3 real shift real t1 real t2 real tm(ni) real tt(ni) real x1 real x2 real xbot(nbot) real xl(ni) real xnew real y1 real y2 real ybot(nbot) real yl(ni) real ymin real ynew ! if ( ynew < ybot(1) .or. ybot(nbot) < ynew ) then write ( *, * ) ' ' write ( *, * ) 'CUBIC - Fatal error!' write ( *, * ) ' YNEW = ',ynew write ( *, * ) ' outside of range YBOT(1) = ',ybot(1) write ( *, * ) ' through YBOT(NBOT) = ',ybot(nbot) stop end if xl(1:ni) = 0.0E+00 yl(1:ni) = 0.0E+00 tt(1:ni) = 0.0E+00 tm(1:ni) = 0.0E+00 yl(1) = ybot(1)+(ybot(1)-ybot(3)) yl(2) = ybot(2)+(ybot(1)-ybot(3)) do i = 3,nbot+2 yl(i) = ybot(i-2) xl(i) = xbot(i-2) end do yl(nbot+3) = yl(nbot+1)+(yl(nbot+2)-yl(nbot)) yl(nbot+4) = yl(nbot+2)+(yl(nbot+2)-yl(nbot)) shift = (xl(3)-xl(4))/(yl(3)-yl(4))-(xl(4)-xl(5))/(yl(4)-yl(5)) xl(2) = xl(3)+(yl(2)-yl(3))*(shift+(xl(3)-xl(4))/(yl(3)-yl(4))) xl(1) = xl(2)+(yl(1)-yl(2))*(shift+(xl(2)-xl(3))/(yl(2)-yl(3))) shift = (xl(nbot+2)-xl(nbot+1))/(yl(nbot+2)-yl(nbot+1)) & -(xl(nbot+1)-xl(nbot))/(yl(nbot+1)-yl(nbot)) xl(nbot+3) = xl(nbot+2)+(yl(nbot+3)-yl(nbot+2))* & (shift+(xl(nbot+2)-xl(nbot+1))/(yl(nbot+2)-yl(nbot+1))) xl(nbot+4) = xl(nbot+3)+(yl(nbot+4)-yl(nbot+3))* & (shift+(xl(nbot+3)-xl(nbot+2))/(yl(nbot+3)-yl(nbot+2))) do i = 1, nbot+3 tm(i) = (xl(i+1)-xl(i)) / (yl(i+1)-yl(i)) end do do i = 3, nbot+2 dome = abs ( tm(i+1) - tm(i) ) + abs ( tm(i-1) + tm(i-2) ) if ( dome < 1.0e-10) then tt(i) = 0.0E+00 else tt(i) = (abs(tm(i+1)-tm(i))*tm(i-1) & +abs(tm(i-1)-tm(i-2))*tm(i))/dome end if end do ! ! Find the node YL(IST) which is nearest to YNEW. ! ymin = abs(yl(3)-ynew) ist = 3 do i = 3,nbot+2 if ( abs(yl(i)-ynew) < ymin ) then ist = i ymin = abs(yl(i)-ynew) end if end do if ( (yl(ist)-ynew) > 0.0E+00 ) then y1 = yl(ist-1) x1 = xl(ist-1) y2 = yl(ist) x2 = xl(ist) t1 = tt(ist-1) t2 = tt(ist) else y1 = yl(ist) x1 = xl(ist) y2 = yl(ist+1) x2 = xl(ist+1) t1 = tt(ist) t2 = tt(ist+1) end if ! ! Evaluate the spline at YNEW. ! p2 = (3.0*(x2-x1)/(y2-y1)-2.0*t1-t2)/(y2-y1) p3 = (t1+t2-2.0*(x2-x1)/(y2-y1))/(y2-y1)**2 xnew = x1+t1*(ynew-y1)+p2*(ynew-y1)**2+p3*(ynew-y1)**3 return end subroutine diflow ( acof, diff, flow ) ! !*********************************************************************** ! !! DIFLOW computes the convection-diffusion coefficient. ! ! ! Discussion: ! ! DIFLOW computes ACOF, given FLOW, the mass velocity RHO*U, ! and DIFF, the value of GAMMA/DELX. ! ! Several schemes are available, but currently the power law is used. ! ! See Patankar, chapter 5. ! implicit none ! real acof real diff real flow ! acof = diff if ( diff == 0.0E+00 ) then return end if ! ! Power Law. ! ! Define ! FE = rho*u ! DE = gamma/delx ! Then the coefficient is ! ! AE = DE * Max(0, (1-0.1*Abs(FE)/DE)**5) + Max(0,-FE) ! acof = diff * max(1.0d-10,(1.0d0-0.1*abs(flow/diff))**5) if ( flow < 0.0 ) acof = acof - flow return end subroutine doarea ( area, areal, areas, areat, icrys, jcrys, l0, m0, & ni, nj, xc, yc ) ! !*********************************************************************** ! !! DOAREA computes the area of the control volumes. ! ! ! Discussion: ! ! DOAREA is given (XC,YC), the locations of the "corners" of the ! control volumes, and calculates the area of each control volume. ! ! DOAREA uses the fact that the area of a polygon can be computed ! by ! ! AREA = 0.5 * SUM (I=1 to N) X(I) * (Y(I+1)-Y(I-1)) ! ! where Y(0) should be replaced by Y(N), and Y(N+1) by Y(1). ! ! Using this formula, AREA may come out negative, depending on ! whether the nodes are given in clockwise or counterclockwise order. ! implicit none ! integer ni integer nj ! real area(ni,nj) real areal real areas real areat integer i integer icrys integer j integer jcrys integer l0 integer m0 integer nbad real xc(ni,nj) real yc(ni,nj) ! do i = 1,l0 do j = 1,m0 if ( i == 1 .or. i == l0 .or. j == 1 .or. j == m0 ) then area(i,j) = 0.0 else area(i,j) = 0.5*( & xc(i, j)* (yc(i+1,j) -yc(i, j+1)) & +xc(i+1,j)* (yc(i+1,j+1)-yc(i, j)) & +xc(i+1,j+1)*(yc(i, j+1)-yc(i+1,j)) & +xc(i, j+1)*(yc(i, j) -yc(i+1,j+1))) end if end do end do ! ! Check for illegal zero length sides. ! nbad = 0 do i = 2, l0-1 do j = 2, m0-1 if ( xc(i,j) == xc(i+1,j) .and. yc(i,j) == yc(i+1,j) ) then nbad = nbad+1 write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' Zero length side for cell I,J:',i,j write ( *, * ) 'XC,YC(I,J) = XC,YC(I+1,J)=',xc(i,j),yc(i,j) else if ( xc(i+1,j) == xc(i+1,j+1) .and. yc(i+1,j) == yc(i+1,j+1) ) then nbad = nbad+1 write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' Zero length side for cell I,J:',i,j write ( *, * ) 'XC,YC(I+1,J) = XC,YC(I+1,J+1)=',xc(i+1,j),yc(i+1,j) else if ( xc(i+1,j+1) == xc(i,j+1) .and. yc(i+1,j+1) == yc(i,j+1) ) then nbad = nbad+1 write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' Zero length side for cell I,J:',i,j write ( *, * ) 'XC,YC(I+1,J+1) = XC,YC(I,J+1)=',xc(i+1,j+1),yc(i+1,j+1) else if ( xc(i,j+1) == xc(i,j) .and. yc(i,j+1) == yc(i,j) ) then nbad = nbad+1 write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' Zero length side for cell I,J:',i,j write ( *, * ) 'XC,YC(I,J+1) = XC,YC(I,J)=',xc(i,j+1),yc(i,j+1) end if end do end do if ( nbad > 0 ) then write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' A total of ',nbad,' zero length cell sides.' stop end if ! ! Check for illegal zero area cells. ! nbad = 0 do i = 2, l0-1 do j = 2, m0-1 if ( area(i,j) == 0.0 ) then nbad = nbad+1 write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' Zero area for cell I = ',i,' J=',j end if end do end do if ( nbad > 0 ) then write ( *, * ) ' ' write ( *, * ) 'SETARE - Fatal error!' write ( *, * ) ' A total of ',nbad,' null cells.' stop end if ! ! Compute liquid, solid, and total areas. ! areal = 0.0 do i = 1,icrys-1 do j = 1,m0 areal = areal+area(i,j) end do end do areas = 0.0 do i = icrys,l0 do j = 1,jcrys areas = areas+area(i,j) end do end do areat = sum ( area(1:l0,1:m0) ) return end subroutine findp ( j, xold, yold, xnew, ynew, xc, yc ) ! !*********************************************************************** ! !! FINDP finds the boundary nodes for a Neumann BC. ! ! ! Discussion: ! ! A Neumann boundary condition is applied to the derivative of a quantity. ! ! FINDP is given a point (XOLD,YOLD). ! ! FINDP returns a pair of values (XNEW,YNEW), which represent ??? ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! real a1 real ai real aj real alp real b1 real bet real c1 real dels real denm real dxb1 real dxb2 real dyb1 real dyb2 integer j real r0 real r1 real r11 real r2 real r22 real r33 real slpi real xc(ni,nj) real xnew real xnew1 real xnew2 real xold real yc(ni,nj) real yn1 real yn2 real ynew real yold ! ! Consider the three points P, Q and R, which have a constant ! KSI coordinate, and an increasing ETA coordinate: ! ! P = ( XC(2,J-1), YC(2,J-1) ) ! Q = ( XC(2,J), YC(2,J) ) ! R = ( XC(2,J+1), YC(2,J+1) ) ! ! ! ^ R = [2,J+1] ! | ! E ! T Q = [2,J] ! A ! | ! | P = [2,J-1] ! | ! +------KSI-------------> ! ! The slope of the line from P to Q is (YQ-YP)/(XQ-XP), and ! the slope of the line from Q to R is (YR-YQ)/(XR-XQ). ! ! If these slopes are close enough, we can assume the three points ! lie on a straight line. So look at the size of ! DENM = (YQ-YP)*(XR-XQ)-(YR-YQ)*(XQ-XP). ! dxb1 = xc(2,j)-xc(2,j-1) dxb2 = xc(2,j+1)-xc(2,j) dyb1 = yc(2,j)-yc(2,j-1) dyb2 = yc(2,j+1)-yc(2,j) denm = dyb2*dxb1-dyb1*dxb2 ! if ( abs(denm) < 0.0001E+00 ) then if ( abs(dxb1) >= 0.001 ) then slpi = -dyb1 / dxb1 aj = yc(2,j) + slpi*xc(2,j) ai = xold - slpi*yold ynew = (aj-ai*slpi) / (1.0+slpi*slpi) xnew = slpi*ynew+ai else ynew = yold xnew = xc(2,j) end if ! ! Use 3 points to find a circle. ! (X-ALP)**2+(Y-BET)**2 = R0 ! ! ! Find the cross point between circle and straight line ! x = ai+slpi*y get a1*y**2-2*b1*y+c1=0. ! else r11 = xc(2,j-1)**2+yc(2,j-1)**2 r22 = xc(2,j)**2+yc(2,j)**2 r33 = xc(2,j+1)**2+yc(2,j+1)**2 r1 = 0.5*(r22-r11) r2 = 0.5*(r33-r22) alp = (r1*dyb2-r2*dyb1)/denm bet = (r2*dxb1-r1*dxb2)/denm r0 = (alp-xc(2,j))**2+(bet-yc(2,j))**2 slpi = (xold-alp)/(yold-bet) ai = alp-slpi*bet a1 = 1.0 + slpi*slpi b1 = slpi * (alp-ai)+bet c1 = (alp-ai)**2 + bet*bet-r0 dels = sqrt ( b1*b1 - a1 * c1 ) yn1 = (b1+dels) / a1 yn2 = (b1-dels) / a1 xnew1 = slpi*yn1+ai xnew2 = slpi*yn2+ai r1 = (xnew1-xc(2,j))**2+(yn1-yc(2,j))**2 r2 = (xnew2-xc(2,j))**2+(yn2-yc(2,j))**2 if ( r1 < r2 ) then ynew = yn1 else ynew = yn2 end if xnew = slpi*ynew+ai end if return end subroutine flux ( ae1, ae2, aim, aip, ajm, ajp, ak1, ak2, ap, b1jbl, & b2jbl, cappa, cd, cmu, con, epsil, f, fjeta, fjksi, fmax, fn, fo, & gam, heta, hksi, icrys, isol, iter, izone, jcrys, l1, lblk, lconv, & lortho, m1, mode, nf, nsolve, ntimes, r, relax, res, rueta, ruksi ) ! !*********************************************************************** ! !! FLUX computes the value of the flux of various quantities. ! ! ! Discussion: ! ! It looks like if ISOL is 0, you just go through the big loop once. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! real acof real ae1(ni,nj) real ae2(ni,nj) real aim(ni,nj) real aip(ni,nj) real ajm(ni,nj) real ajp(ni,nj) real ak1(ni,nj) real ak2(ni,nj) real ap(ni,nj) real ap0(ni,nj) real b1jbl(nj) real b2jbl(nj) real b3jbl(ni) real cappa real cd real cmu real con(ni,nj) real con0(ni,nj) real diff real epsil real error real f(ni,nj,ns) real fjbi1(nj,ns) real fjbj1(ni,ns) real fjbl1(nj,ns) real fjbm1(ni,ns) real fjeta(ni,nj) real fjksi(ni,nj) real flow real fmax(ns) real fn(ni,nj,ns) real fo(ni,nj,ns) real gam(ni,nj) real heta(ni,nj) real hksi(ni,nj) integer i integer ii integer icrys integer isol integer iter integer izone integer j integer jcrys integer jj integer l1 logical lblk(ns) logical lconv logical lortho integer m1 integer mode integer nf integer nsolve(ns) integer nt integer ntimes(ns) real pt(nmaxij) real qt(nmaxij) real r(ni,nj) real relax(nk) real res(ns) real rueta(ni,nj) real ruksi(ni,nj) real t1 real t2 real t3 real t4 real tem1 real tem2 real xje real xjn real xjw real xjs real xme real xmn real xmw real xms ! ! Save copies of CON and AP. ! con0(1:l1,1:m1) = con(1:l1,1:m1) ap0(1:l1,1:m1) = ap(1:l1,1:m1) ! ! I think ITER = 0 occurs only for Temperature in the Solid zone... ! if ( iter == 0 ) then if ( isol /= 0 ) then call solve1(aim,aip,ajm,ajp,ap,cappa,cd,cmu,con,f,fo, & heta,hksi,icrys,izone,jcrys,l1,m1,nf) call solve2(aim,aip,ajm,ajp,ap,con,f,l1,lblk,m1,nf,nsolve) end if return end if do nt = 1,ntimes(nf) ! ! Find the flux on the first step only. ! if ( nt == 1 ) then do j = 1,m1 diff = gam(2,j)/(0.5*hksi(2,j)) flow = -ruksi(2,j) call diflow(acof,diff,flow) qt(2) = ruksi(2,j)*f(2,j,nf)+acof*(f(1,j,nf)-f(2,j,nf)) do i = 3,l1-1 jj = j if ( j == 1) jj = 2 if ( j == m1) jj = m1-1 diff = gam(i,jj)*gam(i-1,jj)/(0.5*hksi(i,jj)*gam(i-1,jj) & +0.5*hksi(i-1,jj)*gam(i,jj)) flow = -ruksi(i,j) call diflow(acof,diff,flow) qt(i) = ruksi(i,j)*f(i,j,nf)+acof*(f(i-1,j,nf)-f(i,j,nf)) end do diff = gam(l1-1,j)/(0.5*hksi(l1-1,j)) flow = ruksi(l1,j) call diflow(acof,diff,flow) qt(l1) = ruksi(l1,j)*f(l1-1,j,nf)+acof*(f(l1-1,j,nf)-f(l1,j,nf)) do i = 2,l1 fjksi(i,j) = qt(i) end do end do ! ! Along J. ! do i = 1,l1 diff = gam(i,2)/(0.5*heta(i,2)) flow = -rueta(i,2) call diflow(acof,diff,flow) qt(2) = rueta(i,2)*f(i,2,nf)+acof*(f(i,1,nf)-f(i,2,nf)) do j = 3,m1-1 if ( i == 1 ) then ii = 2 else if ( i == l1 ) then ii = l1-1 else ii = i end if diff = gam(ii,j)*gam(ii,j-1)/(0.5*heta(ii,j)*gam(ii,j-1) & +0.5*heta(ii,j-1)*gam(ii,j)) flow = -rueta(i,j) call diflow(acof,diff,flow) qt(j) = rueta(i,j)*f(i,j,nf)+acof*(f(i,j-1,nf)-f(i,j,nf)) end do diff = gam(i,m1-1)/(0.5*heta(i,m1-1)) flow = rueta(i,m1) call diflow(acof,diff,flow) qt(m1) = rueta(i,m1)*f(i,m1-1,nf)+acof*(f(i,m1-1,nf)-f(i,m1,nf)) do j = 2,m1 fjeta(i,j) = qt(j) end do end do ! ! Finish jksi and jeta calculation. ! Construct boundary flux for j type boundary condition. ! Boundary condition con(1,j),*,*,* go into jksi and jeta ! ak2/ak1,* is project on xi - direction from eta - direction ! variable, in program let ak2 = 0.0, because sym. orth. ! do j = 2,m1-1 if ( gam(1,j) == 0.0 ) then fjksi(2,j) = (con(1,j)*heta(1,j)+0.5*ak2(2,j) & *(fjeta(1,j)+fjeta(1,j+1)))/ak1(2,j) end if if ( gam(l1,j) == 0.0 ) then fjksi(l1,j) = (con(l1,j)*heta(l1,j)+0.5*ak2(l1,j) & *(fjeta(l1,j)+fjeta(l1,j+1)))/ak1(l1,j) end if end do do i = 2,l1-1 if ( gam(i,1) == 0.0 ) then fjeta(i,2) = (con(i,1)*hksi(i,1)+0.5*ae2(i,2) & *(fjksi(i,1)+fjksi(i+1,1)))/ae1(i,2) end if if ( gam(i,m1) == 0.0 ) then fjeta(i,m1) = (con(i,m1)*hksi(i,m1)+0.5*ae2(i,m1) & *(fjksi(i,m1)+fjksi(i+1,m1)))/ae1(i,m1) end if end do ! ! (2-31) and (2-42) find out vector J cdot vector n or J_n in bound. ! do i = 2,l1-1 t1 = 0.5*(fjksi(i,1)+fjksi(i+1,1)) t2 = 0.5*(fjksi(i,m1)+fjksi(i+1,m1)) fjbj1(i,nf) = (fjeta(i,2)*ae1(i,2)-t1*ae2(i,2))/hksi(i,1) fjbm1(i,nf) = (fjeta(i,m1)*ae1(i,m1)-t2*ae2(i,m1))/hksi(i,m1) if ( mode == 1 ) then fjbj1(i,nf) = fjbj1(i,nf)/r(i,2) fjbm1(i,nf) = fjbm1(i,nf)/r(i,m1) end if end do do j = 2,m1-1 t1 = 0.5*(fjeta(1,j)+fjeta(1,j+1)) t2 = 0.5*(fjeta(l1,j)+fjeta(l1,j+1)) fjbi1(j,nf) = (fjksi(2,j)*ak1(2,j)-t1*ak2(2,j))/heta(1,j) fjbl1(j,nf) = (fjksi(l1,j)*ak1(l1,j) & -t2*ak2(l1,j))/heta(l1,j) if ( mode == 1 ) then fjbi1(j,nf) = fjbi1(j,nf)/r(2,j) fjbl1(j,nf) = fjbl1(j,nf)/r(l1,j) end if end do ! ! This code is only for the temperature variable. ! if ( nf == 5 ) then do j = 2,jcrys-1 t1 = gam(icrys-1,j)/hksi(icrys-1,j)*(f(icrys-1,j,5)-f(icrys,j,5)) b1jbl(j) = t1*ak1(icrys,j)/heta(icrys,j) if ( mode == 1 ) then b1jbl(j) = b1jbl(j)/r(icrys,j) end if end do b1jbl(jcrys) = b1jbl(jcrys+1) b1jbl(m1) = b1jbl(m1-1) do j = 2,jcrys-1 t1 = gam(icrys+1,j)/hksi(icrys+1,j)*(f(icrys,j,5)-f(icrys+1,j,5)) b2jbl(j) = t1*ak1(icrys,j)/heta(icrys,j) if ( mode == 1 ) then b2jbl(j) = b2jbl(j)/r(icrys,j) end if end do b2jbl(jcrys) = b2jbl(jcrys+1) b2jbl(m1) = b2jbl(m1-1) do j = jcrys+1,m1-1 t1 = gam(icrys-1,j)/hksi(icrys-1,j)*(f(icrys-1,j,5)-f(icrys,j,5)) t2 = gam(icrys-1,j)/heta(icrys,j)*(f(icrys,j+1,5)-f(icrys,j,5)) b3jbl(j) = (t1*ak1(icrys,j)-t2*ak2(icrys,j))/heta(icrys,j) if ( mode == 1 ) then b3jbl(j) = b3jbl(j)/r(icrys,j) end if end do b3jbl(jcrys) = b3jbl(jcrys+1) b3jbl(m1) = b3jbl(m1-1) end if else ! ! This code is done if this is NOT the first iteration. ! ! correct j from known phi and j ! correct jksi ! do j = 2,m1-1 if ( gam(1,j) /= 0.0 ) then fjksi(2,j) = fjksi(2,j)+aim(2,j)/ak1(2,j)* & (f(1,j,nf)-f(2,j,nf))+ruksi(2,j)*f(2,j,nf) end if do i = 3,l1 if ( gam(l1,j) /= 0.0 .or. i.ne.l1 ) then fjksi(i,j) = fjksi(i,j)+aip(i-1,j)/ak1(i,j)*(f(i-1,j,nf) & -f(i,j,nf))+ruksi(i,j)*f(i-1,j,nf) end if end do end do ! ! Correct FJETA. ! do i = 2,l1-1 if ( gam(i,1) /= 0.0 ) then fjeta(i,2) = fjeta(i,2)+ajm(i,2)/ae1(i,2)* & (f(i,1,nf)-f(i,2,nf))+rueta(i,2)*f(i,2,nf) end if do j = 3,m1 if ( gam(i,m1) /= 0.0 .or. j.ne.m1 ) then fjeta(i,j) = fjeta(i,j)+ajp(i,j-1)/ae1(i,j)* & (f(i,j-1,nf)-f(i,j,nf))+rueta(i,j)*f(i,j-1,nf) end if end do end do end if ! ! End of "Is this the first iteration or not" block. ! ! Restore old values of CON and AP. ! do i = 1,l1 do j = 1,m1 con(i,j) = con0(i,j) ap(i,j) = ap0(i,j) end do end do if ( .not.lortho ) then do j = 2,m1-1 do i = 1,l1 pt(i) = 0.5*(fjeta(i,j)+fjeta(i,j+1)) qt(i) = 0.5*(rueta(i,j)+rueta(i,j+1)) end do do i = 2,l1-1 t1 = hksi(i,j)/(hksi(i+1,j)+hksi(i,j)) t2 = hksi(i+1,j)/(hksi(i+1,j)+hksi(i,j)) t3 = hksi(i,j)/(hksi(i-1,j)+hksi(i,j)) t4 = hksi(i-1,j)/(hksi(i-1,j)+hksi(i,j)) xje = t2*pt(i)+t1*pt(i+1) xme = t2*qt(i)+t1*qt(i+1) xjw = t4*pt(i)+t3*pt(i-1) xmw = t4*qt(i)+t3*qt(i-1) con(i,j) = con(i,j)+(xje-f(i,j,nf)*xme)*ak2(i+1,j) & -(xjw-f(i,j,nf)*xmw)*ak2(i,j) end do end do do i = 2,l1-1 do j = 1,m1 pt(j) = 0.5*(fjksi(i,j)+fjksi(i+1,j)) qt(j) = 0.5*(ruksi(i,j)+ruksi(i+1,j)) end do do j = 2,m1-1 tem1 = heta(i,j+1)+heta(i,j) tem2 = heta(i,j-1)+heta(i,j) t1 = heta(i,j)/tem1 t2 = heta(i,j+1)/tem1 t3 = heta(i,j)/tem2 t4 = heta(i,j-1)/tem2 xjn = t2*pt(j)+t1*pt(j+1) xmn = t2*qt(j)+t1*qt(j+1) xjs = t4*pt(j)+t3*pt(j-1) xms = t4*qt(j)+t3*qt(j-1) con(i,j) = con(i,j)+(xjn-f(i,j,nf)*xmn)*ae2(i,j+1) & -(xjs-f(i,j,nf)*xms)*ae2(i,j) end do end do end if ! ! calculate jhat ! step 5 use equation (2-97) and (2-102), (2-103), (2-104) ! calculate boundary data from J and other purpose for b_sp. ! from here to end. pc = jksi hat or jeta hat ! ! jksi hat ! do j = 2,m1-1 f(2,j,4) = 0.0 if ( gam(1,j) == 0.0 ) then t1 = fjksi(2,j)-ruksi(2,j)*f(2,j,nf) diff = gam(2,j)/(0.5*hksi(2,j)) flow = -ruksi(2,j) call diflow(acof,diff,flow) if ( acof /= 0.0 ) then f(1,j,nf) = f(2,j,nf)+t1/acof else f(1,j,nf) = f(2,j,nf) end if f(2,j,4) = fjksi(2,j) end if f(l1,j,4) = 0.0 if ( gam(l1,j) == 0.0 ) then t1 = fjksi(l1,j)-ruksi(l1,j)*f(l1-1,j,nf) diff = gam(l1-1,j)/(0.5*hksi(l1-1,j)) flow = ruksi(l1,j) call diflow(acof,diff,flow) if ( acof /= 0.0 ) then f(l1,j,nf) = f(l1-1,j,nf)-t1/acof else f(l1,j,nf) = f(l1-1,j,nf) end if f(l1,j,4) = fjksi(l1,j) end if do i = 3,l1-1 f(i,j,4) = 0.0E+00 end do end do do j = 2,m1-1 do i = 2,l1 fjksi(i,j) = f(i,j,4) end do end do ! ! jeta hat ! do i = 2,l1-1 f(i,2,4) = 0.0 if ( gam(i,1) == 0.0 ) then t1 = fjeta(i,2)-rueta(i,2)*f(i,2,nf) diff = gam(i,2)/(0.5*heta(i,2)) flow = -rueta(i,2) call diflow(acof,diff,flow) if ( acof /= 0.0 ) then f(i,1,nf) = f(i,2,nf)+t1/acof else f(i,1,nf) = f(i,2,nf) end if f(i,2,4) = fjeta(i,2) end if f(i,m1,4) = 0.0 if ( gam(i,m1) == 0.0 ) then t1 = fjeta(i,m1)-rueta(i,m1)*f(i,m1-1,nf) diff = gam(i,m1-1)/(0.5*heta(i,m1-1)) flow = rueta(i,m1) call diflow(acof,diff,flow) if ( acof /= 0.0 ) then f(i,m1,nf) = f(i,m1-1,nf)-t1/acof else f(i,m1,nf) = f(i,m1-1,nf) end if f(i,m1,4) = fjeta(i,m1) end if do j = 3,m1-1 f(i,j,4) = 0.0 end do end do do i = 2,l1-1 do j = 2,m1 fjeta(i,j) = f(i,j,4) end do end do ! ! Calculate and solve phy equation ! (2-106) and (2-107) b = b_s + b_no + b_sp, t1=b_sp ! do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)+fjksi(i,j)*ak1(i,j)-fjksi(i+1,j)*ak1(i+1,j) & +fjeta(i,j)*ae1(i,j)-fjeta(i,j+1)*ae1(i,j+1) end do end do if ( isol == 0) return do i = 2,l1-1 do j = 2,m1-1 ap(i,j) = ap(i,j)/relax(nf) con(i,j) = con(i,j)+(1.0-relax(nf))*ap(i,j)*f(i,j,nf) end do end do call solve1(aim,aip,ajm,ajp,ap,cappa,cd,cmu,con,f,fo, & heta,hksi,icrys,izone,jcrys,l1,m1,nf) call solve2(aim,aip,ajm,ajp,ap,con,f,l1,lblk,m1,nf,nsolve) ! ! BEGIN NEW. ! ! Compute the largest solution component. ! fmax(nf) = 0.0 do i = 1,l1 do j = 1,m1 fmax(nf) = max(fmax(nf),abs(f(i,j,nf))) end do end do ! ! Compare current and previous iterates. ! if ( nt > 1 ) then res(nf) = 0.0 do i = 1,l1 do j = 1,m1 error = abs(f(i,j,nf)-fn(i,j,nf)) if ( fmax(nf) > 0.0 ) then error = error/fmax(nf) end if res(nf) = max(res(nf),error) end do end do end if ! ! Save the new solution in FN. ! do i = 1,l1 do j = 1,m1 fn(i,j,nf) = f(i,j,nf) end do end do if ( res(nf) < epsil .and. nt > 1 ) then ! ! write ( *, * ) 'FLUX - Convergence for variable ',nf,' on step ',nt ! write ( *, * ) ' RES(NF) = ',res(nf),' FMAX(NF)=',fmax(nf) ! return end if end do ! ! End of big iteration ! ! Compute the largest solution component. ! ! fmax(nf) = 0.0 ! do i = 1,l1 ! do j = 1,m1 ! fmax(nf) = max(fmax(nf),abs(f(i,j,nf))) ! end do ! end do ! ! Compute the largest relative change in the solution. ! ! res(nf) = 0.0 ! do i = 1,l1 ! do j = 1,m1 ! ! error = abs(f(i,j,nf)-fn(i,j,nf)) ! ! if ( fmax(nf) > 0.0 ) then ! error = error/fmax(nf) ! end if ! ! res(nf) = max(res(nf),error) ! ! end do ! end do ! ! Decide whether to declare convergence or not. ! if ( res(nf) > epsil ) then lconv = .false. end if ! ! Save the new solution in FN. ! ! do i = 1,l1 ! do j = 1,m1 ! fn(i,j,nf) = f(i,j,nf) ! end do ! end do ! nt = ntimes(nf) ! write ( *, * ) 'FLUX - No convergence for variable ',nf,' on step ',nt ! write ( *, * ) ' RES(NF) = ',res(nf),' FMAX(NF)=',fmax(nf) return end subroutine gamsor ( ap, birad, ce1, ce2, cmu, con, ewall, f, fcsl, fksl, & fma, fo, frsl, gam, gamt, grash, hamag, heta, hksi, icrys, inturb, & izone, jcrys, l1, lsolve, m1, mode, nf, pr, r, rdtm, re, recb, rect, & rho, rhocon, rpr, sige, sigk, sigt, stel, stes, tal, tas, tf, tw, & x, xc, y, yc ) ! !*********************************************************************** ! !! GAMSOR sets the coefficients for the transport problems. ! ! ! Discussion: ! ! These coefficients are associated with U, V, T, W, TK, TE and E. ! ! Note that GAM(I,1) = 0 means the gradient is zero at wall J=1. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! real ap(ni,nj) real ar real birad real ce1 real ce2 real cm4 real cmu real con(ni,nj) real dcdr real depdr real depdx real dudx(ni,nj) real dudy(ni,nj) real dvdx(ni,nj) real dvdy(ni,nj) real dwdx(ni,nj) real dwdy(ni,nj) real ewall real f(ni,nj,ns) real fcsl real fksl real fma real fo(ni,nj,ns) real frsl real gam(ni,nj) real gamt(ni,nj) real gdudx(ni,nj) real gdudy(ni,nj) real gdvdx(ni,nj) real gdvdy(ni,nj) real grash real hamag real heta(ni,nj) real hksi(ni,nj) integer i integer icrys integer inturb integer izone integer j integer jcrys integer l1 logical lsolve(ns) integer m1 integer mode integer modelm integer nf real p11 real p12 real p21 real p22 real p31 real p32 real p41 real p42 real pbig real pr real prod(ni,nj) real r(ni,nj) real rdtm real re real recb real rect real rho(ni,nj) real rhocon real rpr real sige real sigk real sigt real stel real stes real tal real tas real tf real tw real x(ni,nj) real xc(ni,nj) real xplus(nmaxij) real y(ni,nj) real yc(ni,nj) real yplus(nmaxij) ! ap(1:l1,1:m1) = - rdtm * rho(1:l1,1:m1) con(1:l1,1:m1) = 0.0E+00 gam(1:l1,1:m1) = gamt(1:l1,1:m1) + rhocon / re ! ! Horizontal velocity. ! if ( nf == 1 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,1)*rho(i,j)*rdtm+rho(i,j)*f(i,j,5)*grash/re**2 end do end do if ( inturb == 1 ) then do i = 1,l1 do j = 1,m1 gdudx(i,j) = 0.0 gdudy(i,j) = 0.0 gdvdx(i,j) = 0.0 gdvdy(i,j) = 0.0 end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,f(1,1,1),dudx(i,j),dudy(i,j),xc,yc) call gradnt(heta,hksi,i,j,f(1,1,2),dvdx(i,j),dvdy(i,j),xc,yc) gdudx(i,j) = gamt(i,j)*dudx(i,j) gdudy(i,j) = gamt(i,j)*dudy(i,j) gdvdx(i,j) = gamt(i,j)*dvdx(i,j) gdvdy(i,j) = gamt(i,j)*dvdy(i,j) end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,gdudx,p11,p12,xc,yc) call gradnt(heta,hksi,i,j,gdudy,p21,p22,xc,yc) call gradnt(heta,hksi,i,j,gdvdx,p31,p32,xc,yc) call gradnt(heta,hksi,i,j,gdvdy,p41,p42,xc,yc) con(i,j) = con(i,j)+p11+p32 end do end do cm4 = cmu**0.25 if ( mode == 0 ) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5*heta(i,2)/rhocon if ( yplus(i) <= 11.63 ) then gam(i,1) = rhocon/re else gam(i,1) = yplus(i)/(2.5*log(ewall*yplus(i))) end if end do end if do i = 2,l1-1 yplus(i) = re*rho(i,m1-1)*sqrt(f(i,m1-1,7))*cm4*0.5*heta(i,m1-1)/rhocon if ( yplus(i) <= 11.63 ) then gam(i,m1) = rhocon/re else gam(i,m1) = yplus(i)/(2.5*log(ewall*yplus(i))) end if end do do j = 2,m1-1 xplus(j) = re*rho(2,j)*sqrt(f(2,j,7))*cm4*0.5*hksi(2,j)/rhocon if ( xplus(j) <= 11.63 ) then gam(1,j) = rhocon/re else gam(1,j) = xplus(j)/(2.5*log(ewall*xplus(j))) end if end do do j = 2,m1-1 xplus(j) = re*rho(l1-1,j)*sqrt(f(l1-1,j,7))*cm4*0.5*hksi(l1-1,j)/rhocon if ( xplus(j) <= 11.63 ) then gam(l1,j) = rhocon/re else gam(l1,j) = xplus(j)/(2.5*log(ewall*xplus(j))) end if end do end if do i = 2, l1-1 gam(i,1) = 0.0E+00 end do do j = 2, m1-1 f(1,j,1) = 0.0 f(l1,j,1) = 0.0 end do f(l1,1,1) = 0.0 do j = jcrys+1, m1-1 f(l1,j,1) = fo(l1,j,2)*(xc(l1,j+1)-xc(l1,j))/(yc(l1,j+1)-yc(l1,j)) end do f(l1,jcrys+1,1) = 0.0 ! ! Vertical velocity. ! else if ( nf == 2 ) then do i = 2, l1-1 do j = 2, m1-1 con(i,j) = fo(i,j,2)*rho(i,j)*rdtm ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re if ( mode == 1 ) then con(i,j) = con(i,j)+rho(i,j)*f(i,j,6)**2/r(i,j)**3 ap(i,j) = ap(i,j)-rho(i,j)/r(i,j)**2 end if end do end do if ( inturb == 1 ) then do i = 1, l1 do j = 1, m1 gdudx(i,j) = 0.0 gdudy(i,j) = 0.0 gdvdx(i,j) = 0.0 gdvdy(i,j) = 0.0 end do end do do i = 2, l1-1 do j = 2, m1-1 call gradnt(heta,hksi,i,j,f(1,1,1),dudx(i,j),dudy(i,j),xc,yc) call gradnt(heta,hksi,i,j,f(1,1,2),dvdx(i,j),dvdy(i,j),xc,yc) gdudx(i,j) = gamt(i,j)*dudx(i,j) gdudy(i,j) = gamt(i,j)*dudy(i,j) gdvdx(i,j) = gamt(i,j)*dvdx(i,j) gdvdy(i,j) = gamt(i,j)*dvdy(i,j) end do end do do i = 2, l1-1 do j = 2, m1-1 call gradnt(heta,hksi,i,j,gdudx,p11,p12,xc,yc) call gradnt(heta,hksi,i,j,gdudy,p21,p22,xc,yc) call gradnt(heta,hksi,i,j,gdvdx,p31,p32,xc,yc) call gradnt(heta,hksi,i,j,gdvdy,p41,p42,xc,yc) con(i,j) = con(i,j)+p21+p42 end do end do cm4 = cmu**0.25E+00 if ( mode == 0 ) then do i = 2, l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5*heta(i,2)/rhocon gam(i,1) = rhocon/re if ( yplus(i) > 11.63 ) then gam(i,1) = yplus(i) / (2.5*log(ewall*yplus(i))) end if end do end if do i = 2, l1-1 yplus(i) = re * rho(i,m1-1) * sqrt ( f(i,m1-1,7) ) * cm4 & * 0.5 * heta(i,m1-1) / rhocon gam(i,m1) = rhocon / re if ( yplus(i) > 11.63 ) then gam(i,m1) = yplus(i) / (2.5*log(ewall*yplus(i))) end if end do do j = 2, m1-1 xplus(j) = re*rho(2,j)*sqrt(f(2,j,7))*cm4*0.5*hksi(2,j)/rhocon gam(1,j) = rhocon / re if ( xplus(j) > 11.63 ) then gam(1,j) = xplus(j) / (2.5*log(ewall*xplus(j))) end if end do do j = 2, m1-1 xplus(j) = re*rho(l1-1,j)*sqrt(f(l1-1,j,7))*cm4*0.5*hksi(l1-1,j)/rhocon gam(l1,j) = rhocon/re if ( xplus(j) > 11.63 ) then gam(l1,j) = xplus(j)/(2.5*log(ewall*xplus(j))) end if end do end if do i = 2,l1-1 f(i,1,2) = 0.0E+00 end do do j = jcrys+1,m1-1 dcdr = (f(l1,j+1,5)-f(l1,j,5))/(y(l1,j+1)-y(l1,j)) f(l1,j,2) = f(l1-1,j,2)+(xc(l1,j)-xc(l1-1,j))*dcdr*(fma/(re*pr)) end do f(l1,jcrys+1,2) = 0.0E+00 do j = 2,jcrys f(1,j,2) = 0.0 f(l1,j,2) = 0.0 end do f(l1,m1-1,2) = 0.0 f(l1,m1,2) = 0.0 ! ! Temperature computations in the solid zone. ! else if ( nf == 5 .and. izone == 1 ) then do i = 1,l1 do j = 1,m1 gam(i,j) = rpr*fksl/fcsl/frsl end do end do do i = icrys+1,l1 gam(i,1) = 0.0 end do do i = 2, l1-1 do j = 2, m1-1 con(i,j) = fo(i,j,5)*rho(i,j)*rdtm end do end do do j = 2, jcrys f(l1,j,5) = -stes / stel end do do i = icrys+1, l1 f(i,m1,5) = f(i,m1-1,5)-birad*((f(i,m1,5)*(tw-tf)+tf)**4-tas**4) & *(y(i,m1)-y(i,m1-1))/100000000.0 end do ! ! Temperature calculations in the liquid zone. ! else if ( nf == 5 .and. izone == 2 ) then do i = 1, l1 do j = 1, m1 gam(i,j) = rpr + gamt(i,j) / sigt end do end do do i = 2, l1-1 do j = 2, m1-1 con(i,j) = fo(i,j,5)*rho(i,j)*rdtm end do end do if ( inturb == 1 ) then cm4 = cmu**0.25 pbig = 9.24*((pr/sigt)**0.75-1.0)*(1.0+0.28*exp(-0.007*pr/sigt)) if ( mode == 0) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5*heta(i,2)/rhocon gam(i,1) = rpr if ( yplus(i) > 11.63 ) then gam(i,1) = max(rhocon/re, & yplus(i)/(sigt*(2.5*log(ewall*yplus(i))+pbig))) end if end do end if do i = 2, l1-1 yplus(i) = re*rho(i,m1-1)*sqrt(f(i,m1-1,7))*cm4*0.5*heta(i,m1-1)/rhocon gam(i,m1) = rpr if ( yplus(i) > 11.63 ) then gam(i,m1) = max(rhocon/re,yplus(i)/ & (sigt*(2.5*log(ewall*yplus(i))+pbig))) end if end do do j = 2, m1-1 xplus(j) = re*rho(2,j)*sqrt(f(2,j,7))*cm4*0.5*hksi(2,j)/rhocon gam(1,j) = rpr if ( xplus(j) > 11.63 ) then gam(1,j) = max(rhocon/re, & xplus(j)/(sigt*(2.5*log(ewall*xplus(j))+pbig))) end if end do do j = 2, m1-1 xplus(j) = re*rho(l1-1,j)*sqrt(f(l1-1,j,7)) & *cm4*0.5*hksi(l1-1,j) / rhocon gam(l1,j) = rpr if ( xplus(j) > 11.63 ) then gam(l1,j) = max(rhocon/re, & xplus(j)/(sigt*(2.5*log(ewall*xplus(j))+pbig))) end if end do end if do j = 2, jcrys f(l1,j,5) = 0.0 end do do j = jcrys+1, m1 f(l1,j,5) = f(l1-1,j,5)-0.1*birad*((f(l1,j,5)*(tw-tf)+tf)**4-tal**4) & *(x(l1,j)-x(l1-1,j))/100000000.0 end do do j = 1, m1 f(1,j,5) = 0.5+0.5*yc(2,j) end do do i = 1, l1 f(i,m1,5) = 1.0 end do do i = 1, l1 gam(i,1) = 0.0 end do ! ! Angular momentum. ! else if ( nf == 6 ) then ! ! modelm = 1. based on rectangular de/dx of Sabhapathy and Salcudean ! modelm = 2. based on curvilinear de/dx of Sabhapathy and Salcudean. ! modelm = 3. same as 2 but different numerical treatment for source ! terms based on strong magnetic. 1,2 is good for weak Ha ! 3 is good for moderate. ! modelm = 4. based on strong magnetic. rotation is eliminated. ! modelm = 5. is a test. ! if ( hamag == 0.0 ) then modelm = 0.0 else if ( hamag > 0.0 .and. hamag <= 20.0 ) then modelm = 2 else if ( hamag > 20.0 .and. hamag <= 40.0 ) then modelm = 5 else if ( hamag > 40.0 ) then modelm = 4 end if do i = 2, l1-1 do j = 2, m1-1 con(i,j) = fo(i,j,6)*rho(i,j)*rdtm if ( modelm == 1 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)* & (f(i+1,j,9)-f(i,j,9))/(x(i+1,j)-x(i,j)) else if ( modelm == 2 ) then call gradnt(heta,hksi,i,j,f(1,1,9),depdx,depdr,xc,yc) if ( j >= 3 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)*depdx end if else if ( modelm == 3 ) then call gradnt(heta,hksi,i,j,f(1,1,9),depdx,depdr,xc,yc) if ( j >= 3 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)*depdx & +rho(i,j)*hamag**2/re*f(i,j,6) ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re end if else if ( modelm == 4 ) then ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re else if ( modelm == 5 ) then call gradnt(heta,hksi,i,j,f(1,1,9),depdx,depdr,xc,yc) if ( j >= 20 ) then con(i,j) = con(i,j)-rho(i,j)*hamag**2/re/r(i,j)*depdx end if if ( j <= 20 ) then ap(i,j) = ap(i,j)-rho(i,j)*hamag**2/re end if end if ! ! -2/r*dOmega/dr /re ! ar = 2.0/r(i,j)/(y(i,j)-y(i,j-1))/re con(i,j) = con(i,j)+rho(i,j)*ar*f(i,j-1,6) ap(i,j) = ap(i,j)-rho(i,j)*ar end do end do do i = 1, l1 f(i,m1,6) = recb f(i,1,6) = 0.0 end do do j = jcrys+1, m1 gam(l1,j) = 0.0 end do do j = 1, m1 f(1,j,6) = recb*r(2,j)**2 if ( j > jcrys ) then f(l1,j,6) = f(l1-1,j,6) else f(l1,j,6) = rect*r(l1,j)**2 end if end do ! ! Turbulent kinetic energy. ! else if ( nf == 7 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,7)*rdtm gam(i,j) = rhocon/re+gamt(i,j)/sigk end do end do do i = 2,l1-1 do j = 2,m1-1 call gradnt(heta,hksi,i,j,f(1,1,1),dudx(i,j),dudy(i,j),xc,yc) call gradnt(heta,hksi,i,j,f(1,1,2),dvdx(i,j),dvdy(i,j),xc,yc) prod(i,j) = gamt(i,j)*(2.0*(dudx(i,j)**2+dvdy(i,j)**2) & +(dudy(i,j)+dvdx(i,j))**2) if ( mode == 1 ) then prod(i,j) = prod(i,j)+gamt(i,j)*2.0*(f(i,j,2)/r(i,j))**2 end if if ( lsolve(6) ) then call gradnt(heta,hksi,i,j,f(1,1,6),dwdx(i,j),dwdy(i,j),xc,yc) prod(i,j) = prod(i,j)+gamt(i,j)*(dwdx(i,j)**2 & +(dwdy(i,j)-2.0*f(i,j,6)/r(i,j))**2)/r(i,j)**2 end if con(i,j) = con(i,j)+prod(i,j) ap(i,j) = ap(i,j)-cmu*rho(i,j)**2*abs(f(i,j,7))/(gamt(i,j)+1.e-5) end do end do ! ! Turbulent dissipation. ! else if ( nf == 8 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = fo(i,j,8)*rdtm gam(i,j) = rhocon/re+gamt(i,j)/sige end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)+ce1*cmu*rho(i,j)*abs(f(i,j,7))*prod(i,j) & /((gamt(i,j)+1.e-5)) ap(i,j) = ap(i,j)-ce2*cmu*rho(i,j)**2*abs(f(i,j,7)) & /(gamt(i,j)+1.e-5) end do end do end if return end subroutine gradnt ( heta, hksi, i, j, phi, dphidx, dphidy, xc, yc ) ! !*********************************************************************** ! !! GRADNT calculates gradients at the primary nodes. ! ! ! Discussion: ! ! GRADNT is given the value of a state quantity PHI at the primary nodes, ! the KSI and ETA spacing between primary nodes, and the coordinates of ! the corner nodes that define the control volumes. ! ! GRADNT calculates the gradient (dPHI/dX, dPHI/dY) at the primary nodes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! real detadx real detady real dpdeta real dpdksi real dphidx real dphidy real dxdeta real dxdksi real dksidx real dksidy real dydeta real dydksi real heta(ni,nj) real hksi(ni,nj) integer i integer j real phi(ni,nj) real phia real phib real phic real phid real t1 real t2 real volume real xc(ni,nj) real xd real xm real xp real xu real yc(i,j) real yd real ym real yp real yu ! ! Interpolate the value of PHI at the cell volume interfaces ! from its value at primary nodes. ! ! ^ (I,J+1) ! | ! | PHI(B) ! E ! T (I-1,J) PHI(C) (I,J) PHI(A) (I+1,J) ! A ! | PHI(D) ! | ! | (I,J-1) ! | ! +-------------------KSI-------------> ! t1 = hksi(i,j)/(hksi(i,j)+hksi(i+1,j)) t2 = hksi(i+1,j)/(hksi(i,j)+hksi(i+1,j)) phia = t1*phi(i+1,j)+t2*phi(i,j) t1 = heta(i,j)/(heta(i,j)+heta(i,j+1)) t2 = heta(i,j+1)/(heta(i,j)+heta(i,j+1)) phib = t1*phi(i,j+1)+t2*phi(i,j) t1 = hksi(i,j)/(hksi(i,j)+hksi(i-1,j)) t2 = hksi(i-1,j)/(hksi(i,j)+hksi(i-1,j)) phic = t1*phi(i-1,j)+t2*phi(i,j) t1 = heta(i,j)/(heta(i,j)+heta(i,j-1)) t2 = heta(i,j-1)/(heta(i,j)+heta(i,j-1)) phid = t1*phi(i,j-1)+t2*phi(i,j) ! ! Now subtract opposing values, assuming a nominal spacing of ! 1 for KSI and ETA, to get dPHI/dKSI and dPHI/dETA. ! dpdksi = phia-phic dpdeta = phib-phid ! ! Now, using the coordinates of the "corner nodes", locate the midpoints ! of the control volume interfaces that surround primary node (I,J). ! ! ^ ! | [I,J+1] Up [I+1,J+1] ! E ! T Minus (I,J) Plus ! A ! | [I,J] Down [I+1,J] ! | ! +-----XSI----> ! xp = 0.5 * (xc(i+1,j+1)+xc(i+1,j)) yp = 0.5 * (yc(i+1,j+1)+yc(i+1,j)) xm = 0.5 * (xc(i,j+1)+xc(i,j)) ym = 0.5 * (yc(i,j+1)+yc(i,j)) xu = 0.5 * (xc(i,j+1)+xc(i+1,j+1)) yu = 0.5 * (yc(i,j+1)+yc(i+1,j+1)) xd = 0.5 * (xc(i,j)+xc(i+1,j)) yd = 0.5 * (yc(i,j)+yc(i+1,j)) ! ! From the coordinates of the control volume wall centers, ! estimate dX/dKSI, dX/dETA, dY/dKSI and dY/dETA at the primary node, ! again assuming a spacing of 1 for KSI and ETA between control ! volume wall centers. ! dxdksi = xp-xm dxdeta = xu-xd dydksi = yp-ym dydeta = yu-yd ! ! Compute the determinant of the Jacobian (XSI,ETA)-->(X,Y). ! ! J(XSI,ETA) = ( dX/dKSI dX/dETA) ! ( dY/dKSI dY/dETA) ! volume = dxdksi*dydeta-dxdeta*dydksi ! ! Now compute the elements of the inverse Jacobian. ! dksidx = dydeta/volume dksidy = -dxdeta/volume detadx = -dydksi/volume detady = dxdksi/volume ! ! Now finally compute dPHI/dX and dPHI/dY using the chain rule. ! dphidx = dpdksi*dksidx+dpdeta*detadx dphidy = dpdeta*detady+dpdksi*dksidy return end subroutine inidat ( ae1, ae2, ak1, ak2, area, b, birad, bo, cappa, cd, ce1, & ce2, cfo, cinc, cmu, cost, cvn, delt, dtm, epsad, epsil, ewall, f, fcsl, & fks,fksl,fma,fn,fnu,fo,fr,frsl,gamt,grash,hamag,heta,hf,hksi, & icost,icrys,inturb,iplot,ipref,iprint,jcrys, & jpref,l0,l1,last,lastt,lblk,lortho,lsolve,m0,m1,maxbot,mode, & nbot,ndt,ni,nj,nk,np,npar,npc,ns,nsolve,ntimes,orth,par,pr,ra,rdtm, & re,recb,rect,relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk, & sigma,sigt,smooth,stel,stes,tal,tanca,tanca2,tas,tend,tf,tinit, & title,tnow,tw,vave,vol,x,xbot,xc,xlen,y,ybot,yc,ylen) ! !*********************************************************************** ! !! INIDAT sets the initial values of certain data. ! ! ! Modified: ! ! 22 March 2002 ! ! Parameters: ! ! Output, integer IPLOT, controls restart file creation. ! 0, do not create a restart file. ! 1, create a restart file. ! implicit none ! integer maxbot integer ni integer nj integer nk integer npar integer ns ! real ae1(ni,nj) real ae2(ni,nj) real ak1(ni,nj) real ak2(ni,nj) real area(ni,nj) real b real birad real bo real cappa real cd real ce1 real ce2 real cfo real cinc real cl real cmu real cost real cs real cvn real delt real dtm real epsad real epsil real ewall real f(ni,nj,ns) real fcsl real fkl real fks real fksl real fma real fn(ni,nj,ns) real fnu real fo(ni,nj,ns) real fr real frsl real gamt(ni,nj) real grash real grav real hamag real heta(ni,nj) real hf real hksi(ni,nj) integer i integer icost integer icrys integer inturb integer iplot integer ipref integer iprint integer j integer jcrys integer jpref integer k integer l0 integer l1 integer last integer lastt logical lblk(ns) logical lortho logical lsolve(ns) integer m0 integer m1 integer mode integer nbot integer ndt integer np integer npc integer nsolve(ns) integer ntimes(ns) real orth real par(npar) real pr real ra real rdtm real re real recb real rect real relax(nk) real rho(ni,nj) real rhocon real rhol real rhos real rpr real rueta(ni,nj) real ruksi(ni,nj) real sige real sigk real sigma real sigt real smooth real stel real stes real tal real tanca real tanca2 real tas real tend real tf real tin real tinit character ( len = 25 ) title(ns) real tnow real tw real vave real vol(ni,nj) real x(ni,nj) real xbot(maxbot) real xc(ni,nj) real xlen real y(ni,nj) real ybot(maxbot) real yc(ni,nj) real ylen ! ! Set things. ! ae1(1:ni,1:nj) = 0.0 ae2(1:ni,1:nj) = 0.0 ak1(1:ni,1:nj) = 0.0 ak2(1:ni,1:nj) = 0.0 area(1:ni,1:nj) = 0.0 b = 0.1778 cappa = 0.4187 cd = 1.0 ce1 = 1.44 ce2 = 1.92 cinc = 0.0 cl = 1000.0 cmu = 0.09 cost = 0.0 cs = 1000.0 cvn = 200000.0 dtm = 100.0 epsad = 0.0001 epsil = 0.0001 ewall = 9.793 f(1:ni,1:nj,1:ns) = 0.0 fkl = 64.0 fks = 22.0 fma = -1000.0 fn(1:ni,1:nj,1:ns) = 0.0 fnu = 3.0e-7 fo(1:ni,1:nj,1:ns) = 0.0 fr = 1.0e-10 grash = 10000000.0 grav = 9.81 hamag = 0.0 heta(1:ni,1:nj) = 0.0 hf = 1800000.0 hksi(1:ni,1:nj) = 0.0 ! ! ICOST = 1, U**2+V**2 ! 2, (T-TF)**2 ! 3, (VAVE-SQRT(U**2+V**2)) ! icost = 1 ! ! TEMPORARY ! icost = 3 icrys = 42 inturb = 0 iplot = 1 ipref = 11 ! ! IPRINT = 0, don't print very much. ! 1, print intermediate information. ! iprint = 0 jcrys = 22 jpref = 30 l0 = 62 l1 = 62 ! ! Interstep iteration number. ! Use LAST = 40 for more accurate results. ! Use LAST = 10 for quicker results. ! last = 40 ! ! Set the number of timesteps. ! ! For tests, set LASTT = 2. ! For a real run, try LASTT = 200. ! lastt = 2 lblk(1:4) = .false. lblk(5:10) = .true. lortho = .false. lsolve(1:ns) = .false. m0 = 42 m1 = 42 mode = 1 ndt = 0 np = 3 npc = 4 ! ! Number of linear iterations. ! nsolve(1) = 3 nsolve(2) = 3 nsolve(3) = 1 nsolve(4) = 1 nsolve(5) = 3 nsolve(6) = 3 nsolve(7) = 3 nsolve(8) = 3 nsolve(9) = 1 nsolve(10) = 1 ! ! Number of nonlinear iterations. ! ntimes(1) = 5 ntimes(2) = 5 ntimes(3) = 1 ntimes(4) = 1 ntimes(5) = 3 ntimes(6) = 3 ntimes(7) = 3 ntimes(8) = 3 ntimes(9) = 1 ntimes(10) = 1 orth = 50000.0E+00 pr = 0.015E+00 rdtm = 0.0E+00 re = 1.0E+00 recb = 0.0E+00 rect = 0.0E+00 relax(1) = 0.3E+00 relax(2) = 0.3E+00 relax(3) = 0.8E+00 relax(4) = 1.0E+00 relax(5) = 0.7E+00 relax(6) = 0.2E+00 relax(7) = 0.5E+00 relax(8) = 0.5E+00 relax(9) = 1.0E+00 relax(10) = 1.0E+00 relax(11) = 1.0E+00 relax(12) = 0.6E+00 relax(13) = 1.0E+00 relax(14) = 1.0E+00 rhocon = 1.0E+00 rhol = 2490.0E+00 rhos = 2490.0E+00 rueta(1:ni,1:nj) = 0.0E+00 ruksi(1:ni,1:nj) = 0.0E+00 sige = 1.3E+00 sigk = 1.0E+00 sigma = 0.72E+00 sigt = 0.9E+00 smooth = 0.0001E+00 tal = 1523.0E+00 tanca = 1.0E+00 tanca2 = 1.0E+00 tas = 1523.0E+00 tf = 1683.0E+00 tin = 1523.0E+00 tinit = 0.0E+00 title(1) = 'U velocity' title(2) = 'V velocity' title(3) = 'Pressure' title(4) = 'Corrected pressure' title(5) = 'Temperature' title(6) = 'Rotational velocity' title(7) = 'Turbulent dissipation' title(8) = 'Turbulent energy' title(9) = 'Magnetic stream function' title(10) = 'Stream function' ! ! TEMPORARY ! tw = 1713.0E+00 tw = par(1) vave = 1850.0E+00 vol(1:ni,1:nj) = 0.0E+00 x(1:ni,1:nj) = 0.0E+00 xc(1:ni,1:nj) = 0.0E+00 xlen = 1.2E+00 y(1:ni,1:nj) = 0.0E+00 yc(1:ni,1:nj) = 0.0 ylen = 1.0 ! ! Set things that depend on things. ! birad = 0.25*5.67*1.5*0.0254/fks/(tw-tf)/re bo = (rhol*grav*b**2)/sigma cfo = fnu/(b*b) if ( inturb == 1 ) then f(1:ni,1:nj,7) = 0.005E+00 end if fcsl = cs/cl fksl = fks/fkl frsl = rhos/rhol ra = grash*pr rho(1:ni,1:nj) = rhocon rpr = rhocon/(pr*re) stel = cl*(tw-tf)/hf stes = cs*(tf-tin)/hf tend = tinit+dtm*lastt tnow = tinit ! ! Set things that depend on things that depend on things. ! delt = cfo*dtm if ( inturb == 1 ) then f(1:ni,1:nj,8) = f(1:ni,1:nj,7)**1.5 / 0.006E+00 gamt(1:ni,1:nj) = 0.006 * cmu * rho(1:ni,1:nj) * f(1:ni,1:nj,7)**0.5 end if ! ! Set the crucible shape. ! xbot(1) = 0.0 ybot(1) = 0.0 xbot(2) = 0.015 ybot(2) = 0.3 xbot(3) = 0.03 ybot(3) = 0.5 xbot(4) = 0.046 ybot(4) = 0.6 xbot(5) = 0.07 ybot(5) = 0.7 xbot(6) = 0.10 ybot(6) = 0.8 xbot(7) = 0.14 ybot(7) = 0.9 xbot(8) = 0.18 ybot(8) = 0.95 xbot(9) = 0.25 ybot(9) = 1.0 nbot = 9 return end subroutine inigrd ( icrys, jcrys, l0, m0, nbot, xbot, xc, xlen, ybot, & yc, ylen ) ! !*********************************************************************** ! !! INIGRD makes an initial assignment of the grid points XC, YC. ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! real free integer i integer icrys integer j integer jcrys integer l0 integer licrys integer ll1 integer llst1 integer m0 real t1 real t2 real wave real xbot(nbot) real xc(ni,nj) real xlen real ybot(nbot) real yc(ni,nj) real ylen ! ll1 = icrys/2+1 llst1 = jcrys/2+1 licrys = (jcrys+m0)/2+1 free = 0.5E+00 * xlen ! ! Set the Y coordinates ! do i = 2,l0 if ( i < ll1 ) then t2 = free * (0.5*(dble(i-2)/dble(ll1-2))**1.5) else if ( i < icrys ) then t2 = free * (1.0-0.5*(dble(icrys-i)/dble(icrys-ll1))**1.5) else t2 = free + (xlen-free)*(dble(i-icrys)/dble(l0-icrys))**1.2 end if do j = 2, m0 wave = 0.4 + 0.02*sin(t2*20.0*3.14159) if ( j < llst1 ) then t1 = wave/2.0*(dble(j-2)/dble(llst1-2))**1.5 else if ( j < jcrys ) then t1 = wave-wave/2.0*(dble(jcrys-j)/dble(jcrys-llst1))**1.5 else if ( j < licrys ) then t1 = 0.4+0.3*(dble(j-jcrys)/dble(licrys-jcrys))**1.5 else t1 = 1.0-0.3*(dble(m0-j)/dble(m0-licrys))**1.5 end if yc(i,j) = t1*ylen end do end do ! ! Set the XC coordinates along the bottom boundary. ! xc(2,2) = xbot(1) do j = 3, m0-1 call cubic(nbot,yc(2,j),ybot,xc(2,j),xbot) end do xc(2,m0) = xbot(nbot) ! ! Now assign the XC coordinates of the other nodes. ! Here, we explicitly assume that the free surface and crystal ! boundaries occur at a coordinate value of 0.6. ! do j = 2, m0 do i = 3, l0 if ( i < ll1 ) then t1 = 0.5*(dble(i-2)/dble(ll1-2))**1.5 xc(i,j) = (1.0-t1)*xc(2,j)+t1*free else if ( i <= icrys ) then t1 = 1.0E+00 - 0.5E+00 * (dble(icrys-i)/dble(icrys-ll1))**1.5E+00 xc(i,j) = (1.0-t1)*xc(2,j)+t1*free else t1 = (dble(i-icrys)/dble(l0-icrys))**1.2E+00 if ( j < jcrys ) then xc(i,j) = xc(icrys,j)+t1*((xlen-free)+0.5*yc(l0,jcrys)**2 & -0.5*yc(l0,j)**2) else xc(i,j) = xc(icrys,j)+t1*(xlen-free) end if end if end do end do ! ! Copy values to dummy nodes with I = 1 or J=1. ! xc(1,1) = xc(2,2) yc(1,1) = yc(2,2) do j = 2,m0 xc(1,j) = xc(2,j) yc(1,j) = yc(2,j) end do do i = 2,l0 xc(i,1) = xc(i,2) yc(i,1) = yc(i,2) end do return end subroutine initl ( fjeta, fjksi, heta, hksi, l1, m1, ni, nj, p, rho, & rueta, ruksi, u, v, x, y ) ! !*********************************************************************** ! !! INITL estimates quantities at control volume interfaces. ! ! ! Discussion: ! ! INITL estimates the values of P, and of RHO*dU/dKSI and RHO*dU/dETA ! at the interfaces of the control volumes by averaging values associated ! with primary nodes. ! implicit none ! integer ni integer nj ! real fjeta(ni,nj) real fjksi(ni,nj) real heta(ni,nj) real hksi(ni,nj) integer i integer j integer l1 integer m1 real p(ni,nj) real rho(ni,nj) real rhov real rueta(ni,nj) real ruksi(ni,nj) real tt real u(ni,nj) real ucs real ulen real v(ni,nj) real vcs real x(ni,nj) real xcomp real y(ni,nj) real ycomp ! do i = 2,l1 do j = 1,m1 if ( hksi(i,j)+hksi(i-1,j) /= 0.0E+00 ) then tt = hksi(i,j)/(hksi(i,j)+hksi(i-1,j)) else write ( *, * ) 'HKSI(I,J),HKSI(I-1,J) = 0, I=',i,' J=',j tt = 0.0E+00 end if ucs = tt*u(i-1,j)+(1.0-tt)*u(i,j) vcs = tt*v(i-1,j)+(1.0-tt)*v(i,j) fjksi(i,j) = tt*p(i-1,j)+(1.0-tt)*p(i,j) rhov = tt*rho(i-1,j)+(1.0-tt)*rho(i,j) xcomp = x(i,j)-x(i-1,j) ycomp = y(i,j)-y(i-1,j) ulen = sqrt(xcomp**2+ycomp**2) if ( ulen > 0.0E+00 ) then ruksi(i,j) = rhov*(ucs*xcomp+vcs*ycomp)/ulen else ruksi(i,j) = 0.0E+00 end if end do end do do i = 1,l1 do j = 2,m1 if ( heta(i,j)+heta(i,j-1) /= 0.0 ) then tt = heta(i,j)/(heta(i,j)+heta(i,j-1)) else tt = 0.0 write ( *, * ) 'HETA(I,J),HETA(I,J-1) = 0, I=',i,' J=',j end if ucs = tt*u(i,j-1)+(1.0-tt)*u(i,j) vcs = tt*v(i,j-1)+(1.0-tt)*v(i,j) fjeta(i,j) = tt*p(i,j-1)+(1.0-tt)*p(i,j) rhov = tt*rho(i,j-1)+(1.0-tt)*rho(i,j) xcomp = x(i,j)-x(i,j-1) ycomp = y(i,j)-y(i,j-1) ulen = sqrt(xcomp**2+ycomp**2) if ( ulen > 0.0E+00 ) then rueta(i,j) = rhov*(ucs*xcomp+vcs*ycomp)/ulen else rueta(i,j) = 0.0E+00 end if end do end do return end subroutine movgrd ( b1jbl, b2jbl, bo, delt, fksl, fr, frsl, icrys, & iprint, jcrys, l0, m0, p, pr, re, stel, tanca, tanca2, xc, xlen, yc ) ! !*********************************************************************** ! !! MOVGRD calculates the new position of the interface and free surface. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! real ac real as real b1jbl(nj) real b2jbl(nj) real bo real coeff real coff real d2hdy2(nj) real delt real dhdy(nj) real dhdym(nj) real dhdyp(nj) real dlen real dx1 real dxm1 real dy1 real dym1 real fksl real flow(ni) real flomax real flomin real fr real frsl real high(ni) real hilev real himax real himin integer i integer icrys integer interl integer iprint integer j integer jcrys integer k integer l0 integer m0 real omega real omega2 real p(ni,nj) real pcsurf real pr real pres(ni) real pullv real re real rr(nj) real stel real tanca real tanca2 real xb(ni) real xc(ni,nj) real xlen real yc(ni,nj) real zb1(ni) ! do j = 2,m0 zb1(j) = xc(icrys,j) end do do k = 1,5 do j = jcrys+1,m0-1 dhdy(j) = (xc(icrys,j+1)-xc(icrys,j))/(yc(icrys,j+1)-yc(icrys,j)) dhdym(j) = (xc(icrys,j)-xc(icrys,j-1))/(yc(icrys,j)-yc(icrys,j-1)) dhdyp(j) = (xc(icrys,j+1)-xc(icrys,j))/(yc(icrys,j+1)-yc(icrys,j)) end do do j = jcrys+1,m0-1 d2hdy2(j) = (dhdyp(j)-dhdym(j))/(yc(icrys,j+1)-yc(icrys,j)) rr(j) = d2hdy2(j)/(1.0+(dhdy(j)*dhdy(j)))**1.5 & +dhdy(j)/yc(icrys,j)/(1.0+(dhdy(j)*dhdy(j)))**0.5 end do pres(jcrys+1) = 0.0E+00 do j = jcrys+1,m0-2 pres(j+1) = pres(j)+(p(icrys-1,j+1)-p(icrys-1,j))/yc(icrys-1,j+1) end do himax = 0.0E+00 himin = 0.0E+00 interl = (jcrys+m0)/2 do j = jcrys+1,m0-1 pcsurf = pres(j)-pres(interl) high(j) = fr*pcsurf+rr(j)/bo end do hilev = high(interl) omega = 0.0E+00 do j = jcrys+1,m0-1 high(j) = omega*((high(j)-hilev)-(xc(icrys,j)-xc(icrys,interl))) himax = max(himax,high(j)) himin = min(himin,high(j)) xc(icrys,j) = xc(icrys,j)+high(j) end do if ( omega /= 0.0E+00 ) then xc(icrys,jcrys) = xc(icrys,jcrys+1)+ & tanca*(yc(icrys,jcrys+1)-yc(icrys,jcrys)) xc(icrys,m0) = xc(icrys,m0-1) + & tanca2*(yc(icrys,m0)-yc(icrys,m0-1)) end if flomax = 0.0E+00 flomin = 0.0E+00 if ( iprint > 0 ) then if ( k == 5 ) then write ( *, * ) ' ' write ( *, * ) 'MOVGRD:' write ( *, * ) ' Himin = ',himin write ( *, * ) ' Himax = ',himax write ( *, * ) ' Omega = ',omega write ( *, * ) ' ' write(*,'(7f10.4)')high(jcrys+1),high(jcrys+2),high(m0-10), & high(m0-8),high(m0-4),high(m0-2),high(m0-1) write(*,'(7f9.4)')xc(icrys,jcrys),xc(icrys,jcrys+1), & xc(icrys,jcrys+4),xc(icrys,m0-8),xc(icrys,m0-4), & xc(icrys,m0-1),xc(icrys,m0) write ( *, * ) ' ' write ( *, * ) 'MOVGRD:' write ( *, * ) ' Free surface movement' end if end if end do omega2 = 0.0E+00 do j = 2,jcrys-1 dxm1 = xc(icrys,j+1)-xc(icrys,j) dym1 = yc(icrys,j+1)-yc(icrys,j) dlen = sqrt(dxm1*dxm1+dym1*dym1) as = -dxm1/dlen ac = dym1/dlen dx1 = delt*frsl*(b1jbl(j)-fksl*b2jbl(j))*stel/(re*pr)*ac dy1 = delt*frsl*(b1jbl(j)-fksl*b2jbl(j))*stel/(re*pr)*as if ( j == jcrys-1 ) then flow(j) = dx1 else flow(j) = dx1+dy1**2/dx1 end if end do ! ! What is PULLV? ! pullv = flow(jcrys-1) do j = 3,jcrys-1 flomax = max(flomax,(flow(j)-pullv)) flomin = min(flomin,(flow(j)-pullv)) end do omega2 = min(omega2, 0.003/(flomax+1.e-6)) omega2 = min(omega2, 0.003/(-flomin+1.e-6)) do j = 3,jcrys-1 flow(j) = omega2*(flow(j)-pullv)+(xc(icrys,jcrys)-zb1(jcrys)) xb(j) = xc(icrys,j)+flow(j) xc(icrys,j) = xb(j) end do xc(icrys,2) = xc(icrys,3) if ( iprint > 0 ) then write ( *, * ) ' ' write ( *, * ) 'MOVGRD:' write ( *, * ) ' Flomin = ',flomin write ( *, * ) ' Flomax = ',flomax write ( *, * ) ' Omega = ',omega2 write ( *, * ) ' ' write(*,'(2x,7f9.4)')xb(jcrys-1),xb(jcrys-2),xb(8),xb(6),xb(4),xb(3),xb(2) write(*,'(2x,7f9.4)')b1jbl(jcrys-1),-fksl*b2jbl(jcrys-1), & b1jbl(jcrys-2),-fksl*b2jbl(jcrys-2),b1jbl(jcrys-5),-fksl*b2jbl(jcrys-5) write(*,'(2x,7f9.4)')b1jbl(jcrys-1)-fksl*b2jbl(jcrys-1), & b1jbl(jcrys-2)-fksl*b2jbl(jcrys-2),b1jbl(jcrys-5)-fksl*b2jbl(jcrys-5), & b1jbl(jcrys-10)-fksl*b2jbl(jcrys-10),b1jbl(jcrys-15)-fksl*b2jbl(jcrys-15), & b1jbl(jcrys-17)-fksl*b2jbl(jcrys-17),b1jbl(jcrys-19)-fksl*b2jbl(jcrys-19) write(*,'(2x,7f9.4)')flow(21),flow(19),flow(17),flow(15), & flow(13),flow(11),flow(10) write(*,'(2x,7f9.4)')flow(9),flow(8),flow(7),flow(6),flow(5),flow(4),flow(3) write ( *, * ) ' ' write ( *, * ) 'MOVGRD:' write ( *, * ) ' Solid-liquid interface movement' write ( *, * ) ' ' end if do j = 2,m0 coeff = (xc(2,j)-xc(icrys,j))/(xc(2,j)-zb1(j)) coff = (xlen-xc(icrys,j))/(xlen-zb1(j)) do i = 2,icrys-1 xc(i,j) = xc(2,j)-(xc(2,j)-xc(i,j))*coeff end do do i = icrys+1,l0 xc(i,j) = xlen-(xlen-xc(i,j))*coff if ( j <= jcrys .and. i == l0 ) then xc(l0,j) = xlen+0.5E+00*yc(l0,jcrys)**2-0.5E+00*yc(l0,j)**2 end if end do end do return end subroutine output ( cfo, iter, izone, res, smax, ssum, t, tnow, u, w ) ! !*********************************************************************** ! !! OUTPUT prints information about the current solution. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! real cfo integer iter integer izone real res(ns) real smax real ssum real t(ni,nj) real tnow real u(ni,nj) real w(ni,nj) ! ! Solid zone output. ! if ( izone == 1 ) then if ( iter == 0 ) then write ( *, * ) ' ' write ( *, * ) 'OUTPUT:' write ( *, * ) ' Current time = ',tnow write ( *, * ) ' = ',tnow*cfo,' seconds.' write ( *, * ) ' Solid zone results.' write ( *, * ) ' Iter Res(5) T(50,21) T(50,22)' write ( *, * ) ' ' end if write(*,'(i4,2x,10g12.4)')iter,res(5),t(50,21),t(50,22) ! ! Liquid zone output. ! else if ( izone == 2 ) then if ( iter == 1 ) then write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'OUTPUT:' write ( *, * ) ' Liquid zone results.' write ( *, * ) ' SMAX is the maximum local mass imbalance.' write ( *, * ) ' SSUM is the total mass imbalance.' write ( *, * ) ' ' write ( *, * ) ' Iter SMAX SSUM U(5,5) W(5,5) T(5,5)' write ( *, * ) ' ' end if write(*,'(i4,2x,10g12.4)')iter,smax,ssum,u(5,5),w(5,5),t(5,5) end if return end subroutine pmod ( ipref, jcrys, jpref, l1, m1, p ) ! !*********************************************************************** ! !! PMOD extends pressure values from primary nodes to corner nodes. ! ! ! Discussion: ! ! PMOD carries out some simple operations to extend the value of the ! pressure to primary nodes at the corners, and to normalize the ! pressure by subtracting off its value at the reference point. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! integer i integer ipref integer j integer jcrys integer jpref integer l1 integer m1 real p(ni,nj) real pref ! ! Extrapolate to get pressures on the boundary corners. ! p(1,1) = p(2,1)+p(1,2)-p(2,2) p(l1,1) = p(l1-1,1)+p(l1,2)-p(l1-1,2) p(1,jcrys) = p(2,jcrys)+p(1,jcrys-1)-p(2,jcrys-1) p(l1,jcrys) = p(l1-1,jcrys)+p(l1,jcrys-1)-p(l1-1,jcrys-1) ! ! Subtract off the reference pressure. ! pref = p(ipref,jpref) do j = 1,m1 do i = 1,l1 p(i,j) = p(i,j) - pref end do end do return end subroutine prdat ( b, birad, bo, cappa, cfo, cvn, delt, dtm, fcsl, fksl, & fma, fnu, fr, frsl, grash, icost, icrys, inturb, iprint, & jcrys, l0, last, lastt, m0, mode, nbot, ns, nsolve, ntimes, orth, & pr, ra, rdtm, recb, rect, rhocon, smooth, stel, stes, tanca, tanca2, tend, & tf, tinit, title, tw, vave, xlen, ylen ) ! !*********************************************************************** ! !! PRDAT prints out the initial values of certain data. ! implicit none ! integer ns ! real b real birad real bo real cappa real cfo real cvn real delt real dtm real fcsl real fksl real fma real fnu real fr real frsl real grash integer i integer icost integer icrys integer inturb integer iprint integer jcrys integer l0 integer last integer lastt integer m0 integer mode integer nbot integer nsolve(ns) integer ntimes(ns) real orth real pr real ra real rdtm real recb real rect real rhocon real smooth real stel real stes real tanca real tanca2 real tend real tf real tinit character ( len = 25 ) title(ns) real tw real vave real xlen real ylen ! write ( *, * ) ' ' write ( *, * ) 'PRDAT:' write ( *, * ) ' ' write ( *, * ) ' B, Crucible radius = ',b write ( *, * ) ' BIRAD, thermal coefficient = ',birad write ( *, * ) ' BO, the Bond number = ',bo write ( *, * ) ' CAPPA, wall function constant = ',cappa write ( *, * ) ' CFO, time scale = ',cfo write ( *, * ) ' CVN, ADAPT cell volume weight = ',cvn write ( *, * ) ' DELT, time step in seconds = ',delt write ( *, * ) ' DTM, dimensionless time step = ',dtm write ( *, * ) ' FCSL = CS/CL = ',fcsl write ( *, * ) ' FKSL = FKS/FKL = ',fksl write ( *, * ) ' FMA, Marangoni number FMA = ',fma write ( *, * ) ' FR, Froude number = ',fr write ( *, * ) ' FRSL = RHOS/RHOL = ',frsl write ( *, * ) ' GRASH, Grashof number = ',grash write ( *, * ) ' ICOST, cost functional = ',icost write ( *, * ) ' ICRYS, maximum I of crystal = ',icrys write ( *, * ) ' INTURB, turbulence option = ',inturb write ( *, * ) ' IPRINT, printing option = ',iprint write ( *, * ) ' JCRYS, maximum J of crystal = ',jcrys write ( *, * ) ' L0, number of I nodes = ',l0 write ( *, * ) ' LAST, number of zone iterations on ' write ( *, * ) ' each time step = ',last write ( *, * ) ' LASTT, number of time steps = ',lastt write ( *, * ) ' M0, number of J nodes = ',m0 write ( *, * ) ' MODE, 0 cartesian, 1 axisymmetric = ',mode write ( *, * ) ' NBOT, number of boundary points = ',nbot write ( *, * ) ' ORTH, ADAPT orthogonality weight = ',orth write ( *, * ) ' PR, Prandtl number PR = ',pr write ( *, * ) ' RA, Rayleigh number = ',ra write ( *, * ) ' RBD = FMA/RA = ',fma/ra write ( *, * ) ' RDTM, = 1/DTM or 0 = ',rdtm write ( *, * ) ' RECB, crucible Reynolds number = ',recb write ( *, * ) ' RECT, crystal Reynolds number = ',rect write ( *, * ) ' RHOCON, density constant = ',rhocon write ( *, * ) ' SMOOTH, ADAPT smoothness weight = ',smooth write ( *, * ) ' STEL, liquid Stefan number STEL = ',stel write ( *, * ) ' STES, solid Stefan number STES = ',stes write ( *, * ) ' TANCA = ',tanca write ( *, * ) ' TANCA2 = ',tanca2 write ( *, * ) ' TEND, dimensionless end time = ',tend write ( *, * ) ' TF, crystal melting temperature = ',tf write ( *, * ) ' TINIT, dimensionless start time = ',tinit write ( *, * ) ' TW, the wall temperature = ',tw write ( *, * ) ' VAVE, desired average velocity = ',vave write ( *, * ) ' WEBER = FR*BO = ',fr*bo write ( *, * ) ' XLEN = problem region length = ',xlen write ( *, * ) ' YLEN = problem region height = ',ylen write ( *, * ) ' Characteristic velocity FNU/B = ',fnu/b write ( *, * ) ' ' write ( *, * ) 'Variable Nonlinear Linear' write ( *, * ) ' Iterations Iterations' write ( *, * ) ' ' do i = 1,ns write ( *, * ) title(i),ntimes(i),nsolve(i) end do return end subroutine resid ( aim, aip, ajm, ajp, ap, con, f, l1, m1, nf ) ! !*********************************************************************** ! !! RESID computes the linear equation residual at interior primary nodes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! real aim(ni,nj) real aip(ni,nj) real ajm(ni,nj) real ajp(ni,nj) real ap(ni,nj) real con(ni,nj) real conmax real diamax real f(ni,nj,ns) integer i integer imax integer iprint integer j integer jmax integer l1 integer m1 integer nf real offmax real ratio real ratmin real res real resmax real sum real xmax ! iprint = 0 conmax = 0.0E+00 diamax = 0.0E+00 offmax = 0.0E+00 ratmin = 100000.0E+00 resmax = 0.0E+00 xmax = 0.0E+00 imax = 0 jmax = 0 do i = 2,l1-1 do j = 2,m1-1 res = con(i,j)-ap(i,j)*f(i,j,nf)+aim(i,j)*f(i-1,j,nf) & +aip(i,j)*f(i+1,j,nf)+ajm(i,j)*f(i,j-1,nf)+ajp(i,j)*f(i,j+1,nf) if ( abs(res) >= resmax ) then resmax = abs(res) imax = i jmax = j end if if ( abs(con(i,j)) >= conmax ) then conmax = abs(con(i,j)) end if if ( abs(f(i,j,nf)) >= xmax ) then xmax = abs(f(i,j,nf)) end if if ( abs(ap(i,j)) >= diamax ) then diamax = abs(ap(i,j)) end if sum = abs(aim(i,j))+abs(aip(i,j))+abs(ajm(i,j))+abs(ajp(i,j)) if ( sum >= offmax ) then offmax = sum end if if ( sum /= 0.0E+00 ) then ratio = abs(ap(i,j))/sum if ( ratio < ratmin ) then ratmin = ratio end if end if end do end do if ( iprint == 1 ) then write ( *, * ) ' ' write ( *, * ) 'RESID' write ( *, * ) ' Maxixum residual for variable ',nf,' is ',resmax write ( *, * ) ' at I,J = ',imax,jmax i = imax j = jmax write ( *, * ) ' ' write ( *, * ) ' A(I,J)*X(I,J) = ',ap(i,j),f(i,j,nf) write ( *, * ) ' A(I-1,J)*X(I-1,J) = ',aim(i,j),f(i-1,j,nf) write ( *, * ) ' A(I+1,J)*X(I+1,J) = ',aip(i,j),f(i+1,j,nf) write ( *, * ) ' A(I,J-1)*X(I,J-1) = ',ajm(i,j),f(i,j-1,nf) write ( *, * ) ' A(I,J+1)*X(I,J+1) = ',ajp(i,j),f(i,j+1,nf) write ( *, * ) ' CON(I,J) = ',con(i,j) write ( *, * ) ' ' write ( *, * ) ' Norm of variable is ',xmax write ( *, * ) ' Norm of RHS is ',conmax write ( *, * ) ' Norm of diagonal is ',diamax write ( *, * ) ' Norm of off-diag is ',offmax write ( *, * ) ' Min diag/off-diag ',ratmin end if return end subroutine rswrit ( cost, e, gamt, icrys, jcrys, l0, m0, nbot, p, pc, & psi, rueta, ruksi, t, te, tk, tnow, u, v, w, x, xbot, xc, y, ybot, yc ) ! !*********************************************************************** ! !! RSWRIT writes out restart information. ! ! ! Modified: ! ! 22 March 2002 ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! real cost real e(ni,nj) character ( len = 80 ) :: file_name = 'crystal_rs.txt' real gamt(ni,nj) integer i integer icrys integer j integer jcrys integer l0 integer m0 real p(ni,nj) real pc(ni,nj) real psi(ni,nj) real rueta(ni,nj) real ruksi(ni,nj) real t(ni,nj) real te(ni,nj) real tk(ni,nj) real tnow real u(ni,nj) real v(ni,nj) real w(ni,nj) real x(ni,nj) real xbot(nbot) real xc(ni,nj) real y(ni,nj) real ybot(nbot) real yc(ni,nj) ! write ( *, '(a)' )' ' write ( *, '(a)' )'RSWRIT:' write ( *, '(a)' )' Writing restart information to ' // trim ( file_name ) write ( *, '(a)' )' ' open ( unit = 10, file = file_name, status = 'replace' ) write ( 10, '(e12.5)' ) cost write ( 10, '(i6)' ) l0 write ( 10, '(i6)' ) jcrys write ( 10, '(i6)' ) icrys write ( 10, '(i6)' ) m0 write ( 10, '(i6)' ) nbot write ( 10, '(e12.5)' ) tnow write ( 10, '(6e12.5)' ) ((e(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((gamt(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((p(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((pc(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((psi(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((rueta(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((ruksi(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((t(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((te(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((tk(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((u(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((v(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((w(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((x(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) xbot(1:nbot) write ( 10, '(6e12.5)' ) ((xc(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ((y(i,j),i = 1,l0),j=1,m0) write ( 10, '(6e12.5)' ) ybot(1:nbot) write ( 10, '(6e12.5)' ) ((yc(i,j),i = 1,l0),j=1,m0) close ( unit = 10 ) return end subroutine setcst ( area, cinc, icost, icrys, m0, ni, nj, t, tf, u, & v, vave ) ! !*********************************************************************** ! !! SETCST computes a portion of the cost functional. ! ! ! Discussion: ! ! SETCST computes a portion on the cost functional, which is currently ! the integral over time and space of the norm of the fluid velocity. ! ! The integral to be computed for a particular time is: ! ! ICOST = 1: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! (U(T,X,Y)**2 + V(T,X,Y)**2) dX dY ) ! ! ICOST = 2: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! (T(T,X,Y)-TF)**2 dX dY ) ! ! ICOST = 3: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! ( Vave - SQRT(U(T,X,Y)**2 + V(T,X,Y)**2)) dX dY ) ! ! and the total cost is ! ! COST = Integral (T=TINIT to TEND) F(T) dT/ (TEND-TINIT) ! ! Currently, a crude approximation is made to the integrals. ! For instance, the area of the flow region, LIQUID(T), is computed ! by summing up the estimated areas of each control volume. These ! areas, in turn, are estimated by computing the jacobian at the ! center of the control volume. ! ! +-----V-----+ ! | | ! | | ! U N U ! | | ! | | ! +-----V-----+ ! ! Modified: ! ! 22 March 2002 ! ! Author: ! ! John Burkardt ! implicit none ! integer ni integer nj ! real area(ni,nj) real cinc integer i integer icost integer icrys integer j integer m0 real t(ni,nj) real tf real u(ni,nj) real v(ni,nj) real vave ! if ( icost == 1 ) then cinc = 0.0 do i = 1,icrys-1 do j = 1,m0 cinc = cinc+(u(i,j)**2+v(i,j)**2)*area(i,j) end do end do cinc = sqrt(cinc) else if ( icost == 2 ) then cinc = 0.0 do i = 1,icrys-1 do j = 1,m0 cinc = cinc+(t(i,j)-tf)**2*area(i,j) end do end do cinc = sqrt(cinc) else if ( icost == 3 ) then cinc = 0.0E+00 do i = 1,icrys-1 do j = 1,m0 cinc = cinc+area(i,j)*(vave-sqrt(u(i,j)**2+v(i,j)**2))**2 end do end do cinc = sqrt(cinc) end if return end subroutine setgeo ( ae1, ae2, ak1, ak2, heta, hksi, l1, m1, mode, ni, & nj, r, vol, x, xc, y, yc ) ! !*********************************************************************** ! !! SETGEO calculates various geometric quantities. ! ! ! Discussion: ! ! SETGEO is given (XC,YC), the locations of the "corners" of the ! control volumes, and calculates various related geometric parameters, ! including: ! ! X, Y the position of the primary nodes, ! HETA and HKSI, dH/dETA and dH/dKSI, ! the Jacobian, ! AK1, AE1, dALPHA/dKSI, dALPHA/dETA, ! AK2, AE2, dBETA/dKSI, dBETA/dETA, ! R, a weight factor for cylindrical geometries. ! VOL, the volume of the control volumes. ! implicit none ! integer ni integer nj ! real ae1(ni,nj) real ae2(ni,nj) real ak1(ni,nj) real ak2(ni,nj) real dxdeta real dxdksi real dydeta real dydksi real heta(ni,nj) real hksi(ni,nj) integer i integer j integer l1 integer m1 integer mode real r(ni,nj) real t1 real t2 real t3 real vol(ni,nj) real x(ni,nj) real xc(ni,nj) real xd real xjacb real xm real xp real xu real y(ni,nj) real yc(ni,nj) real yd real ym real yp real yu ! ! Calculate the local scale factors HETA and HKSI, and ! the control volumes VOL. do i = 2,l1-1 do j = 2,m1-1 ! ! Compute the coordinates of the mid-side nodes. ! ! ^ ! | [I,J+1] Up [I+1,J+1] ! E ! T Minus (I,J) Plus ! A ! | [I,J] Down [I+1,J] ! | ! +-----XSI----> ! xp = 0.5 * (xc(i+1,j+1)+xc(i+1,j)) yp = 0.5 * (yc(i+1,j+1)+yc(i+1,j)) xm = 0.5 * (xc(i,j+1)+xc(i,j)) ym = 0.5 * (yc(i,j+1)+yc(i,j)) xu = 0.5 * (xc(i,j+1)+xc(i+1,j+1)) yu = 0.5 * (yc(i,j+1)+yc(i+1,j+1)) xd = 0.5 * (xc(i,j)+xc(i+1,j)) yd = 0.5 * (yc(i,j)+yc(i+1,j)) ! ! The mid-side nodes differ by Delta KSI or Delta ETA = 1. ! Use finite differences to estimate derivatives like dX/dKSI. ! dxdksi = xp-xm dxdeta = xu-xd dydksi = yp-ym dydeta = yu-yd ! ! Now compute the length of the line segments that connect ! opposing mid-side nodes. ! hksi(i,j) = sqrt(dxdksi**2+dydksi**2) heta(i,j) = sqrt(dxdeta**2+dydeta**2) ! ! Now compute the area of the control volume, which is just ! the determinant of the Jacobian matrix. ! vol(i,j) = dxdksi*dydeta-dxdeta*dydksi end do end do ! ! Take care of values along the borders J = 1 and J=M1. ! do i = 1,l1 heta(i,1) = 0.0 heta(i,m1) = 0.0 vol(i,1) = 0.0 vol(i,m1) = 0.0 if ( i /= 1 .and. i /= l1 ) then hksi(i,1) = sqrt((xc(i+1,2)-xc(i,2))**2+(yc(i+1,2)-yc(i,2))**2) hksi(i,m1) = sqrt((xc(i+1,m1)-xc(i,m1))**2+(yc(i+1,m1)-yc(i,m1))**2) else hksi(i,1) = 0.0 hksi(i,m1) = 0.0 end if end do ! ! Take care of values along the borders I = 1 and I=L1. ! do j = 1,m1 hksi(1,j) = 0.0 hksi(l1,j) = 0.0 vol(1,j) = 0.0 vol(l1,j) = 0.0 if ( j /= 1 .and. j /= m1 ) then heta(1,j) = sqrt((xc(2,j)-xc(2,j+1))**2+(yc(2,j)-yc(2,j+1))**2) heta(l1,j) = sqrt((xc(l1,j)-xc(l1,j+1))**2+(yc(l1,j)-yc(l1,j+1))**2) else heta(1,j) = 0.0 heta(l1,j) = 0.0 end if end do ! ! Calculate dALPHA/dKSI, dBETA/dKSI, dALPHA/dETA, dBETA/dETA, ! the areas on the control-volume faces. ! do j = 2,m1-1 do i = 2,l1 dxdeta = xc(i,j+1)-xc(i,j) dydeta = yc(i,j+1)-yc(i,j) dxdksi = x(i,j)-x(i-1,j) dydksi = y(i,j)-y(i-1,j) if ( i == 2 .or. i == l1 ) then dxdksi = dxdksi*2.0 dydksi = dydksi*2.0 end if t1 = dxdeta**2+dydeta**2 t2 = dxdksi**2+dydksi**2 t3 = dxdksi*dxdeta+dydeta*dydksi xjacb = dxdksi*dydeta-dxdeta*dydksi if ( xjacb /= 0.0 ) then ak1(i,j) = sqrt(t2)*t1/xjacb ak2(i,j) = sqrt(t1)*t3/xjacb else ak1(i,j) = 0.0 ak2(i,j) = 0.0 end if end do end do do i = 2,l1-1 do j = 2,m1 dxdeta = x(i,j)-x(i,j-1) dydeta = y(i,j)-y(i,j-1) dxdksi = xc(i+1,j)-xc(i,j) dydksi = yc(i+1,j)-yc(i,j) if ( j == 2 .or. j == m1 ) then dxdeta = dxdeta*2.0 dydeta = dydeta*2.0 end if t1 = dxdeta**2+dydeta**2 t2 = dxdksi**2+dydksi**2 t3 = dxdksi*dxdeta+dydeta*dydksi xjacb = dxdksi*dydeta-dxdeta*dydksi if ( xjacb /= 0.0 ) then ae1(i,j) = sqrt(t1)*t2/xjacb ae2(i,j) = sqrt(t2)*t3/xjacb else ae1(i,j) = 0.0 ae2(i,j) = 0.0 end if end do end do ! ! Set the cylindrical geometry weight factor. ! if ( mode == 0 ) then r(1:l1,1:m1) = 1.0E+00 else if ( mode == 1 ) then do i = 1,l1 do j = 1,m1 r(i,j) = y(i,j) ak1(i,j) = ak1(i,j)*r(i,j) ak2(i,j) = ak2(i,j)*r(i,j) ae1(i,j) = ae1(i,j)*r(i,j) ae2(i,j) = ae2(i,j)*r(i,j) vol(i,j) = vol(i,j)*r(i,j) end do end do end if return end subroutine setup ( ae1, ae2, ak1, ak2, b1jbl, b2jbl, birad, cappa, & cd,ce1,ce2,cmu,epsil,ewall,f,fcsl,fjeta,fjksi,fksl,fma, & fmax,fn,fo,frsl,gamt,grash,hamag,heta,hksi,inturb,icrys, & iter,izone,jcrys,l1,lblk,lconv,lortho,lsolve,m1, & mode,nf,np,npc,nsolve,ntimes,pr,r,rdtm,re,recb,rect,res, & relax,rho,rhocon,rpr,rueta,ruksi,sige,sigk,sigt,smax,ssum, & stel, stes, tal, tas, tf, tw, vol, x, xc, y, yc ) ! !*********************************************************************** ! !! SETUP calculates the coefficients of the equations, and solves them. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! real a1 real a2 real acof real ae1(ni,nj) real ae2(ni,nj) real aim(ni,nj) real aip(ni,nj) real ajm(ni,nj) real ajp(ni,nj) real ak1(ni,nj) real ak2(ni,nj) real ap(ni,nj) real apr real b1jbl(nj) real b2jbl(nj) real birad real biu(nj) real biv(nj) real bju(ni) real bjv(ni) real blu(nj) real blv(nj) real bmu(ni) real bmv(ni) real bpi real bpi1 real bpj real bpj1 real bpm real cappa real cd real ce1 real ce2 real cmu real cofu(ni,nj,5) real cofv(ni,nj,5) real con(ni,nj) real conu(ni,nj) real conv(ni,nj) real denom real diff real dpeta real dpksi real dxdeta real dxdksi real dydeta real dydksi real em real ep real epsil real ewall real f(ni,nj,ns) real fcsl real fjeta(ni,nj) real fjksi(ni,nj) real fksl real flow real fma real fmax(ns) real fn(ni,nj,ns) real fo(ni,ni,ns) real frc real frsl real gam(ni,nj) real gamt(ni,nj) real grash real hamag real heta(ni,nj) real hetap(ni,nj) real hksi(ni,nj) real hksip(ni,nj) real hutop(ni,nj) real hvtop(ni,nj) integer i integer icrys integer inturb integer isol integer iter integer izone integer j integer jcrys integer l1 logical lblk(ns) logical lconv logical lortho logical lsolve(ns) integer m1 integer mode integer n integer nf integer np integer npc integer nsolve(ns) integer ntimes(ns) real peta(ni,nj) real pksi(ni,nj) real pr real qt(nmaxij) real r(ni,nj) real rdtm real re real recb real rect real relax(nk) real res(ns) real rho(ni,nj) real rhocon real rpr real rueij real rueij1 real ruet real rueta(ni,nj) real ruki1j real rukij real ruksi(ni,nj) real rukt real sige real sigk real sigt real smax real soor real ssum real stel real stes real t1 real t2 real t3 real t4 real tal real tas real tem1 real tem2 real temp real tf real tmp(ni,nj) real tw real vol(ni,nj) real x(ni,nj) real xc(ni,nj) real xd real xl real xr real xu real y(ni,nj) real yc(ni,nj) real yd real yl real yr real yu ! do n = 1,ns if ( lsolve(n) ) then nf = n ! ! N = 1, U VELOCITY. ! if ( n == 1 ) then call gamsor(ap,birad,ce1,ce2,cmu,con,ewall,f,fcsl,fksl, & fma,fo,frsl,gam,gamt,grash,hamag,heta,hksi,icrys,inturb, & izone,jcrys,l1,lsolve,m1,mode,nf,pr,r,rdtm,re,recb,rect, & rho,rhocon,rpr,sige,sigk,sigt,stel,stes,tal,tas,tf,tw,x, & xc,y,yc) do j = 2,m1-1 diff = 2.0*gam(2,j)/hksi(2,j) flow = -ruksi(2,j) call diflow(acof,diff,flow) aim(2,j) = acof*ak1(2,j) diff = 2.0*gam(l1-1,j)/hksi(l1-1,j) flow = ruksi(l1,j) call diflow(acof,diff,flow) aip(l1-1,j) = acof*ak1(l1,j) do i = 2,l1-2 diff = gam(i,j)*gam(i+1,j)/(0.5*hksi(i,j)*gam(i+1,j) & +0.5*hksi(i+1,j)*gam(i,j)) flow = ruksi(i+1,j) call diflow(acof,diff,flow) aip(i,j) = acof*ak1(i+1,j) aim(i+1,j) = (acof+ruksi(i+1,j))*ak1(i+1,j) end do end do do i = 2,l1-1 diff = 2.0*gam(i,2)/heta(i,2) flow = -rueta(i,2) call diflow(acof,diff,flow) ajm(i,2) = acof*ae1(i,2) diff = 2.0*gam(i,m1-1)/heta(i,m1-1) flow = rueta(i,m1) call diflow(acof,diff,flow) ajp(i,m1-1) = acof*ae1(i,m1) do j = 2,m1-2 diff = gam(i,j)*gam(i,j+1)/(0.5*heta(i,j)*gam(i,j+1) & +0.5*heta(i,j+1)*gam(i,j)) flow = rueta(i,j+1) call diflow(acof,diff,flow) ajp(i,j) = acof*ae1(i,j+1) ajm(i,j+1) = (acof+rueta(i,j+1))*ae1(i,j+1) end do end do do j = 1,m1 biu(j) = gam(1,j) blu(j) = gam(l1,j) end do do i = 1,l1 bju(i) = gam(i,1) bmu(i) = gam(i,m1) end do do i = 1,l1 do j = 1,m1 conu(i,j) = con(i,j) cofu(i,j,1) = ap(i,j) cofu(i,j,2) = aip(i,j) cofu(i,j,3) = aim(i,j) cofu(i,j,4) = ajp(i,j) cofu(i,j,5) = ajm(i,j) end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = (con(i,j)+ap(i,j)*f(i,j,1))*vol(i,j) end do end do isol = 0 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) do i = 2,l1-1 do j = 2,m1-1 hutop(i,j) = (aip(i,j)*f(i+1,j,1)+aim(i,j)*f(i-1,j,1) & +ajp(i,j)*f(i,j+1,1)+ajm(i,j)*f(i,j-1,1)+con(i,j))/vol(i,j) end do end do ! ! N = 2, V VELOCITY. ! else if ( n == 2 ) then call gamsor(ap,birad,ce1,ce2,cmu,con,ewall,f,fcsl,fksl, & fma,fo,frsl,gam,gamt,grash,hamag,heta,hksi,icrys,inturb, & izone,jcrys,l1,lsolve,m1,mode,nf,pr,r,rdtm,re,recb,rect, & rho,rhocon,rpr,sige,sigk,sigt,stel,stes,tal,tas,tf,tw,x,xc,y,yc) do j = 2,m1-1 diff = 2.0*gam(2,j)/hksi(2,j) flow = -ruksi(2,j) call diflow(acof,diff,flow) aim(2,j) = acof*ak1(2,j) diff = 2.0*gam(l1-1,j)/hksi(l1-1,j) flow = ruksi(l1,j) call diflow(acof,diff,flow) aip(l1-1,j) = acof*ak1(l1,j) do i = 2,l1-2 diff = gam(i,j)*gam(i+1,j)/(0.5*hksi(i,j)*gam(i+1,j) & +0.5*hksi(i+1,j)*gam(i,j)) flow = ruksi(i+1,j) call diflow(acof,diff,flow) aip(i,j) = acof*ak1(i+1,j) aim(i+1,j) = (acof+ruksi(i+1,j))*ak1(i+1,j) end do end do do i = 2,l1-1 diff = 2.0*gam(i,2)/heta(i,2) flow = -rueta(i,2) call diflow(acof,diff,flow) ajm(i,2) = acof*ae1(i,2) diff = 2.0*gam(i,m1-1)/heta(i,m1-1) flow = rueta(i,m1) call diflow(acof,diff,flow) ajp(i,m1-1) = acof*ae1(i,m1) do j = 2,m1-2 diff = gam(i,j)*gam(i,j+1)/(0.5*heta(i,j)*gam(i,j+1) & +0.5*heta(i,j+1)*gam(i,j)) flow = rueta(i,j+1) call diflow(acof,diff,flow) ajp(i,j) = acof*ae1(i,j+1) ajm(i,j+1) = (acof+rueta(i,j+1))*ae1(i,j+1) end do end do biv(1:m1) = gam(1,1:m1) blv(1:m1) = gam(l1,1:m1) do i = 1,l1 bjv(i) = gam(i,1) bmv(i) = gam(i,m1) end do do i = 1,l1 do j = 1,m1 conv(i,j) = con(i,j) cofv(i,j,1) = ap(i,j) cofv(i,j,2) = aip(i,j) cofv(i,j,3) = aim(i,j) cofv(i,j,4) = ajp(i,j) cofv(i,j,5) = ajm(i,j) end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = (con(i,j)+ap(i,j)*f(i,j,2))*vol(i,j) end do end do isol = 0 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) do i = 2,l1-1 do j = 2,m1-1 hvtop(i,j) = (aip(i,j)*f(i+1,j,2)+aim(i,j)*f(i-1,j,2) & +ajp(i,j)*f(i,j+1,2)+ajm(i,j)*f(i,j-1,2)+con(i,j))/vol(i,j) dxdksi = 0.5*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5*(yc(i+1,j+1)+yc(i+1,j)-yc(i,j+1)-yc(i,j)) hksip(i,j) = 0.5*(dxdksi*hutop(i,j)+dydksi*hvtop(i,j)) dxdeta = 0.5*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5*(yc(i+1,j+1)+yc(i,j+1)-yc(i+1,j)-yc(i,j)) hetap(i,j) = 0.5*(dxdeta*hutop(i,j)+dydeta*hvtop(i,j)) end do end do ! ! N = 3, P PRESSURE. ! ! HKSIP (stored in CON0), HETAP (stored in AP0) ! are based on momentum interpolation. ! T3 is based on the two-dimenisonal correction of MIS ! T4 is based on the secondary correction of 2-d of MIS ! else if ( n == 3 ) then do i = 2,l1-1 do j = 2,m1-1 con(i,j) = 0.0 ap(i,j) = (aip(i,j)+aim(i,j)+ajm(i,j)+ajp(i,j))/vol(i,j)/rho(i,j) tmp(i,j) = ap(i,j) end do end do if ( .not.lortho ) then do j = 2,m1-1 do i = 1,l1 qt(i) = 0.5 * (rueta(i,j)+rueta(i,j+1)) end do do i = 2,l1-1 tem1 = hksi(i+1,j)+hksi(i,j) tem2 = hksi(i-1,j)+hksi(i,j) t1 = hksi(i,j)/tem1 t2 = hksi(i+1,j)/tem1 t3 = hksi(i,j)/tem2 t4 = hksi(i-1,j)/tem2 con(i,j) = con(i,j)+(t1*qt(i+1)+t2*qt(i))*ak2(i+1,j) & -(t4*qt(i)+t3*qt(i-1))*ak2(i,j) end do end do do i = 2, l1-1 do j = 1, m1 qt(j) = 0.5 * (ruksi(i,j)+ruksi(i+1,j)) end do do j = 2,m1-1 tem1 = heta(i,j+1)+heta(i,j) tem2 = heta(i,j-1)+heta(i,j) t1 = heta(i,j)/tem1 t2 = heta(i,j+1)/tem1 t3 = heta(i,j)/tem2 t4 = heta(i,j-1)/tem2 con(i,j) = con(i,j)+(t1*qt(j+1)+t2*qt(j))*ae2(i,j+1) & -(t4*qt(j)+t3*qt(j-1))*ae2(i,j) end do end do end if do j = 2,m1-1 temp = 0.0 do i = 2,l1-2 a1 = 0.5 * (ak1(i,j)+ak1(i+1,j)) a2 = 0.5 * (ak1(i+1,j)+ak1(i+2,j)) denom = 0.5 * (ap(i,j)*hksi(i,j)/a1+ap(i+1,j)*hksi(i+1,j)/a2) apr = 2.0 * denom/(hksi(i,j)+hksi(i+1,j)) dxdksi = 0.5 * (xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5 * (yc(i+1,j+1)+yc(i+1,j)-yc(i,j+1)-yc(i,j)) rukij = (dxdksi*f(i,j,1)+dydksi*f(i,j,2))*rho(i,j)/hksi(i,j)*a1 dxdksi = 0.5 * (xc(i+2,j+1)+xc(i+2,j)-xc(i+1,j+1)-xc(i+1,j)) dydksi = 0.5 * (yc(i+2,j+1)+yc(i+2,j)-yc(i+1,j+1)-yc(i+1,j)) ruki1j = (dxdksi*f(i+1,j,1)+dydksi*f(i+1,j,2))*rho(i+1,j) & /hksi(i+1,j)*a2 t1 = rukij*(apr*hksi(i+1,j)-ap(i,j)*hksi(i,j)/a1) t2 = ruki1j*(apr*hksi(i,j)-ap(i+1,j)*hksi(i+1,j)/a2) t3 = 0.5 * (t1+t2) if ( i == 2 ) then rukt = ruksi(2,j) else rukt = temp end if temp = ruksi(i+1,j) frc = hksi(i+1,j)/(hksi(i,j)+hksi(i+1,j)) t4 = 0.5*(frc*(ruksi(i+1,j)*ak1(i+1,j)-rukt*ak1(i,j)) & +(1.0-frc)*(ruksi(i+1,j)*ak1(i+1,j)-ruksi(i+2,j)*ak1(i+2,j))) ruksi(i+1,j) = ((hksip(i,j)+hksip(i+1,j)+t3)/denom+t4)/ak1(i+1,j) aip(i,j) = 1.0/denom aim(i+1,j) = aip(i,j) end do aip(l1-1,j) = 0.0 aim(2,j) = 0.0 end do do i = 2,l1-1 do j = 2,m1-2 a1 = 0.5*(ae1(i,j)+ae1(i,j+1)) a2 = 0.5*(ae1(i,j+1)+ae1(i,j+2)) denom = 0.5*(ap(i,j)*heta(i,j)/a1+ap(i,j+1)*heta(i,j+1)/a2) apr = 2.0*denom/(heta(i,j)+heta(i,j+1)) dxdeta = 0.5*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5*(yc(i+1,j+1)+yc(i,j+1)-yc(i+1,j)-yc(i,j)) rueij = (dxdeta*f(i,j,1)+dydeta*f(i,j,2))*rho(i,j)/heta(i,j)*a1 dxdeta = 0.5*(xc(i+1,j+2)+xc(i,j+2)-xc(i+1,j+1)-xc(i,j+1)) dydeta = 0.5*(yc(i+1,j+2)+yc(i,j+2)-yc(i+1,j+1)-yc(i,j+1)) rueij1 = (dxdeta*f(i,j+1,1)+dydeta*f(i,j+1,2))*rho(i,j+1) & /heta(i,j+1)*a2 t1 = rueij*(apr*heta(i,j+1)-ap(i,j)*heta(i,j)/a1) t2 = rueij1*(apr*heta(i,j)-ap(i,j+1)*heta(i,j+1)/a2) t3 = 0.5*(t1+t2) ruet = rueta(i,2) if ( j /= 2) ruet = temp temp = rueta(i,j+1) frc = heta(i,j+1)/(heta(i,j)+heta(i,j+1)) t4 = 0.5*(frc*(rueta(i,j+1)*ae1(i,j+1)-ruet*ae1(i,j)) & +(1.0-frc)*(rueta(i,j+1)*ae1(i,j+1)-rueta(i,j+2)*ae1(i,j+2))) rueta(i,j+1) = ((hetap(i,j)+hetap(i,j+1)+t3)/denom+t4)/ae1(i,j+1) ajp(i,j) = 1.0/denom ajm(i,j+1) = ajp(i,j) end do ajp(i,m1-1) = 0.0 ajm(i,2) = 0.0 end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)+ak1(i,j)*ruksi(i,j)-ak1(i+1,j)*ruksi(i+1,j) & +ae1(i,j)*rueta(i,j)-ae1(i,j+1)*rueta(i,j+1) ap(i,j) = aip(i,j)+aim(i,j)+ajp(i,j)+ajm(i,j) end do end do smax = 0.0 ssum = 0.0 do i = 2,l1-1 do j = 2,m1-1 soor = ap(i,j)*f(i,j,3)-aip(i,j)*f(i+1,j,3) & -aim(i,j)*f(i-1,j,3)-ajp(i,j)*f(i,j+1,3) & -ajm(i,j)*f(i,j-1,3)-con(i,j) ssum = ssum+soor/rho(i,j) smax = max(abs(soor)/1.0d4,smax) end do end do do i = 2,l1-1 do j = 2,m1-1 ap(i,j) = ap(i,j)/relax(np) con(i,j) = con(i,j)+(1.0-relax(np))*ap(i,j)*f(i,j,3) end do end do call solve1(aim,aip,ajm,ajp,ap,cappa,cd,cmu,con,f,fo, & heta,hksi,icrys,izone,jcrys,l1,m1,nf) call solve2(aim,aip,ajm,ajp,ap,con,f,l1,lblk,m1,nf,nsolve) do j = 2,m1-1 do i = 2,l1-2 ruksi(i+1,j) = ruksi(i+1,j)+aip(i,j)/ak1(i+1,j) & *(f(i,j,nf)-f(i+1,j,nf)) end do end do do i = 2,l1-1 do j = 2,m1-2 rueta(i,j+1) = rueta(i,j+1)+ajp(i,j)/ae1(i,j+1) & *(f(i,j,nf)-f(i,j+1,nf)) end do end do do j = 2,m1-1 do i = 2,l1-2 bpi = f(i,j,3) bpi1 = f(i+1,j,3)-hksip(i,j)-hksip(i+1,j) em = hksi(i,j)*tmp(i,j)/(ak1(i,j)+ak1(i+1,j)) ep = hksi(i+1,j)*tmp(i+1,j)/(ak1(i+1,j)+ak1(i+2,j)) t1 = 0.5*em*ep*(ruksi(i+2,j)*ak1(i+2,j) & -ruksi(i,j)*ak1(i,j))/(em+ep) bpm = (bpi*ep+bpi1*em)/(em+ep)+t1 pksi(i+1,j) = bpm+hksip(i,j) end do pksi(l1,j) = 2.0*f(l1-1,j,nf)-pksi(l1-1,j) pksi(2,j) = 2.0*f(2,j,nf)-pksi(3,j) f(1,j,nf) = pksi(2,j) f(l1,j,nf) = pksi(l1,j) end do do i = 2,l1-1 do j = 2,m1-2 bpj = f(i,j,3) bpj1 = f(i,j+1,3)-hetap(i,j)-hetap(i,j+1) em = heta(i,j)*tmp(i,j)/(ae1(i,j)+ae1(i,j+1)) ep = heta(i,j+1)*tmp(i,j+1)/(ae1(i,j+1)+ae1(i,j+2)) t1 = 0.5*em*ep*(rueta(i,j+2)*ae1(i,j+2) & -rueta(i,j)*ae1(i,j))/(em+ep) bpm = (bpj*ep+bpj1*em)/(em+ep)+t1 peta(i,j+1) = bpm+hetap(i,j) end do peta(i,m1) = 2.0*f(i,m1-1,nf)-peta(i,m1-1) peta(i,2) = 2.0*f(i,2,nf)-peta(i,3) f(i,1,nf) = peta(i,2) f(i,m1,nf) = peta(i,m1) end do do j = 1,m1 gam(1,j) = biu(j) gam(l1,j) = blu(j) end do do i = 1,l1 gam(i,1) = bju(i) gam(i,m1) = bmu(i) end do do i = 1,l1 do j = 1,m1 con(i,j) = conu(i,j) ap(i,j) = cofu(i,j,1) aip(i,j) = cofu(i,j,2) aim(i,j) = cofu(i,j,3) ajp(i,j) = cofu(i,j,4) ajm(i,j) = cofu(i,j,5) end do end do do i = 2,l1-1 if ( gam(i,1) == 0.0 ) then ajm(i,2) = 0.0 end if if ( gam(i,m1) == 0.0 ) then ajp(i,m1-1) = 0.0 end if end do do j = 2,m1-1 if ( gam(1,j) == 0.0 ) then aim(2,j) = 0.0 end if if ( gam(l1,j) == 0.0 ) then aip(l1-1,j) = 0.0 end if end do do i = 2,l1-1 do j = 2,m1-1 xr = 0.5*(xc(i+1,j+1)+xc(i+1,j)) yr = 0.5*(yc(i+1,j+1)+yc(i+1,j)) xl = 0.5*(xc(i,j+1)+xc(i,j)) yl = 0.5*(yc(i,j+1)+yc(i,j)) xu = 0.5*(xc(i+1,j+1)+xc(i,j+1)) yu = 0.5*(yc(i+1,j+1)+yc(i,j+1)) xd = 0.5*(xc(i+1,j)+xc(i,j)) yd = 0.5*(yc(i+1,j)+yc(i,j)) dxdksi = xr-xl dydksi = yr-yl dxdeta = xu-xd dydeta = yu-yd dpksi = pksi(i+1,j)-pksi(i,j) dpeta = peta(i,j+1)-peta(i,j) t1 = -dydeta*dpksi+dydksi*dpeta t2 = dxdeta*dpksi-dxdksi*dpeta if ( mode == 1 ) then t1 = t1*r(i,j) t2 = t2*r(i,j) end if tmp(i,j) = t2 con(i,j) = con(i,j)*vol(i,j)+t1 ap(i,j) = -ap(i,j)*vol(i,j)+aip(i,j)+aim(i,j)+ajp(i,j)+ajm(i,j) end do end do nf = 1 isol = 1 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl,cappa, & cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) do j = 1,m1 gam(1,j) = biv(j) gam(l1,j) = blv(j) end do do i = 1,l1 gam(i,1) = bjv(i) gam(i,m1) = bmv(i) end do do i = 1,l1 do j = 1,m1 con(i,j) = conv(i,j) ap(i,j) = cofv(i,j,1) aip(i,j) = cofv(i,j,2) aim(i,j) = cofv(i,j,3) ajp(i,j) = cofv(i,j,4) ajm(i,j) = cofv(i,j,5) end do end do nf = 2 do i = 2,l1-1 if ( gam(i,1) == 0.0 ) then ajm(i,2) = 0.0 end if if ( gam(i,m1) == 0.0 ) then ajp(i,m1-1) = 0.0 end if end do do j = 2,m1-1 if ( gam(1,j) == 0.0 ) then aim(2,j) = 0.0 end if if ( gam(l1,j) == 0.0 ) then aip(l1-1,j) = 0.0 end if end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)*vol(i,j)+tmp(i,j) ap(i,j) = -ap(i,j)*vol(i,j)+aip(i,j)+aim(i,j)+ajp(i,j)+ajm(i,j) end do end do isol = 1 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) ! ! N = 4, PC CORRECTED PRESSURE. ! else if ( n == 4 ) then do j = 1,m1 gam(1,j) = biu(j) gam(l1,j) = blu(j) end do do i = 1,l1 gam(i,1) = bju(i) gam(i,m1) = bmu(i) end do do i = 1,l1 do j = 1,m1 con(i,j) = conu(i,j) ap(i,j) = cofu(i,j,1) aip(i,j) = cofu(i,j,2) aim(i,j) = cofu(i,j,3) ajp(i,j) = cofu(i,j,4) ajm(i,j) = cofu(i,j,5) end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = (con(i,j)+ap(i,j)*f(i,j,1))*vol(i,j) end do end do nf = 1 isol = 0 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) nf = npc do i = 2,l1-1 do j = 2,m1-1 hutop(i,j) = (aip(i,j)*f(i+1,j,1)+aim(i,j)*f(i-1,j,1) & +ajp(i,j)*f(i,j+1,1)+ajm(i,j)*f(i,j-1,1)+con(i,j))/vol(i,j) end do end do do j = 2,m1-1 do i = 2,l1-1 con(i,j) = 0.0 ap(i,j) = 0.0 end do end do do j = 1,m1 gam(1,j) = biv(j) gam(l1,j) = blv(j) end do do i = 1,l1 gam(i,1) = bjv(i) gam(i,m1) = bmv(i) end do do i = 1,l1 do j = 1,m1 con(i,j) = conv(i,j) ap(i,j) = cofv(i,j,1) aip(i,j) = cofv(i,j,2) aim(i,j) = cofv(i,j,3) ajp(i,j) = cofv(i,j,4) ajm(i,j) = cofv(i,j,5) end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = (con(i,j)+ap(i,j)*f(i,j,2))*vol(i,j) end do end do nf = 2 isol = 0 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) nf = npc do i = 2,l1-1 do j = 2,m1-1 hvtop(i,j) = (aip(i,j)*f(i+1,j,2)+aim(i,j)*f(i-1,j,2) & +ajp(i,j)*f(i,j+1,2)+ajm(i,j)*f(i,j-1,2)+con(i,j))/vol(i,j) dxdksi = 0.5*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5*(yc(i+1,j+1)+yc(i+1,j)-yc(i,j+1)-yc(i,j)) hksip(i,j) = 0.5*(dxdksi*hutop(i,j)+dydksi*hvtop(i,j)) dxdeta = 0.5*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5*(yc(i+1,j+1)+yc(i,j+1)-yc(i+1,j)-yc(i,j)) hetap(i,j) = 0.5*(dxdeta*hutop(i,j)+dydeta*hvtop(i,j)) end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = 0.0 ap(i,j) = (aip(i,j)+aim(i,j)+ajm(i,j)+ajp(i,j))/vol(i,j)/rho(i,j) tmp(i,j) = ap(i,j) end do end do if ( .not. lortho ) then do j = 2,m1-1 do i = 1,l1 qt(i) = 0.5*(rueta(i,j)+rueta(i,j+1)) end do do i = 2,l1-1 tem1 = hksi(i+1,j)+hksi(i,j) tem2 = hksi(i-1,j)+hksi(i,j) t1 = hksi(i,j)/tem1 t2 = hksi(i+1,j)/tem1 t3 = hksi(i,j)/tem2 t4 = hksi(i-1,j)/tem2 con(i,j) = con(i,j)+(t1*qt(i+1)+t2*qt(i))*ak2(i+1,j) & -(t4*qt(i)+t3*qt(i-1))*ak2(i,j) end do end do do i = 2,l1-1 do j = 1,m1 qt(j) = 0.5*(ruksi(i,j)+ruksi(i+1,j)) end do do j = 2,m1-1 tem1 = heta(i,j+1)+heta(i,j) tem2 = heta(i,j-1)+heta(i,j) t1 = heta(i,j)/tem1 t2 = heta(i,j+1)/tem1 t3 = heta(i,j)/tem2 t4 = heta(i,j-1)/tem2 con(i,j) = con(i,j)+(t1*qt(j+1)+t2*qt(j))*ae2(i,j+1) & -(t4*qt(j)+t3*qt(j-1))*ae2(i,j) end do end do end if do j = 2,m1-1 temp = 0.0 do i = 2,l1-2 a1 = 0.5*(ak1(i,j)+ak1(i+1,j)) a2 = 0.5*(ak1(i+1,j)+ak1(i+2,j)) denom = 0.5*(ap(i,j)*hksi(i,j)/a1+ap(i+1,j)*hksi(i+1,j)/a2) apr = 2.0*denom/(hksi(i,j)+hksi(i+1,j)) dxdksi = 0.5*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5*(yc(i+1,j+1)+yc(i+1,j)-yc(i,j+1)-yc(i,j)) rukij = (dxdksi*f(i,j,1)+dydksi*f(i,j,2))*rho(i,j)/hksi(i,j)*a1 dxdksi = 0.5*(xc(i+2,j+1)+xc(i+2,j)-xc(i+1,j+1)-xc(i+1,j)) dydksi = 0.5*(yc(i+2,j+1)+yc(i+2,j)-yc(i+1,j+1)-yc(i+1,j)) ruki1j = (dxdksi*f(i+1,j,1)+dydksi*f(i+1,j,2))*rho(i+1,j) & /hksi(i+1,j)*a2 t1 = rukij*(apr*hksi(i+1,j)-ap(i,j)*hksi(i,j)/a1) t2 = ruki1j*(apr*hksi(i,j)-ap(i+1,j)*hksi(i+1,j)/a2) t3 = 0.5*(t1+t2) if ( i == 2 ) then rukt = ruksi(2,j) else rukt = temp end if temp = ruksi(i+1,j) frc = hksi(i+1,j)/(hksi(i,j)+hksi(i+1,j)) t4 = 0.5*(frc*(ruksi(i+1,j)*ak1(i+1,j)-rukt*ak1(i,j)) & +(1.0-frc)*(ruksi(i+1,j)*ak1(i+1,j) & -ruksi(i+2,j)*ak1(i+2,j))) ruksi(i+1,j) = ((hksip(i,j)+hksip(i+1,j)+t3)/denom+t4)/ak1(i+1,j) ruksi(i+1,j) = ruksi(i+1,j)+(f(i,j,3)-f(i+1,j,3))/denom/ak1(i+1,j) aip(i,j) = 1.0/denom aim(i+1,j) = aip(i,j) end do aip(l1-1,j) = 0.0 aim(2,j) = 0.0 end do do i = 2,l1-1 do j = 2,m1-2 a1 = 0.5*(ae1(i,j)+ae1(i,j+1)) a2 = 0.5*(ae1(i,j+1)+ae1(i,j+2)) denom = 0.5*(ap(i,j)*heta(i,j)/a1+ap(i,j+1)*heta(i,j+1)/a2) apr = 2.0*denom/(heta(i,j)+heta(i,j+1)) dxdeta = 0.5*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5*(yc(i+1,j+1)+yc(i,j+1)-yc(i+1,j)-yc(i,j)) rueij = (dxdeta*f(i,j,1)+dydeta*f(i,j,2))*rho(i,j)/heta(i,j)*a1 dxdeta = 0.5*(xc(i+1,j+2)+xc(i,j+2)-xc(i+1,j+1)-xc(i,j+1)) dydeta = 0.5*(yc(i+1,j+2)+yc(i,j+2)-yc(i+1,j+1)-yc(i,j+1)) rueij1 = (dxdeta*f(i,j+1,1)+dydeta*f(i,j+1,2))*rho(i,j+1) & /heta(i,j+1)*a2 t1 = rueij*(apr*heta(i,j+1)-ap(i,j)*heta(i,j)/a1) t2 = rueij1*(apr*heta(i,j)-ap(i,j+1)*heta(i,j+1)/a2) t3 = 0.5*(t1+t2) ruet = rueta(i,2) if ( j /= 2) ruet = temp temp = rueta(i,j+1) frc = heta(i,j+1)/(heta(i,j)+heta(i,j+1)) t4 = 0.5*(frc*(rueta(i,j+1)*ae1(i,j+1)-ruet*ae1(i,j)) & +(1.0-frc)*(rueta(i,j+1)*ae1(i,j+1)-rueta(i,j+2)*ae1(i,j+2))) rueta(i,j+1) = ((hetap(i,j)+hetap(i,j+1)+t3)/denom+t4)/ae1(i,j+1) rueta(i,j+1) = rueta(i,j+1) & +(f(i,j,3)-f(i,j+1,3))/denom/ae1(i,j+1) ajp(i,j) = 1.0/denom ajm(i,j+1) = ajp(i,j) end do ajp(i,m1-1) = 0.0 ajm(i,2) = 0.0 end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)+ak1(i,j)*ruksi(i,j) & -ak1(i+1,j)*ruksi(i+1,j) & +ae1(i,j)*rueta(i,j)-ae1(i,j+1)*rueta(i,j+1) ap(i,j) = aip(i,j)+aim(i,j)+ajp(i,j)+ajm(i,j) end do end do call solve1(aim,aip,ajm,ajp,ap,cappa,cd,cmu,con,f,fo, & heta,hksi,icrys,izone,jcrys,l1,m1,nf) call solve2(aim,aip,ajm,ajp,ap,con,f,l1,lblk,m1,nf,nsolve) do j = 2,m1-1 do i = 2,l1-2 ruksi(i+1,j) = ruksi(i+1,j)+aip(i,j)/ak1(i+1,j) & *(f(i,j,nf)-f(i+1,j,nf)) end do end do do i = 2,l1-1 do j = 2,m1-2 rueta(i,j+1) = rueta(i,j+1)+ajp(i,j)/ae1(i,j+1) & *(f(i,j,nf)-f(i,j+1,nf)) end do end do do j = 2,m1-1 do i = 2,l1-2 bpi = f(i,j,3) bpi1 = f(i+1,j,3)-hksip(i,j)-hksip(i+1,j) em = hksi(i,j)*tmp(i,j)/(ak1(i,j)+ak1(i+1,j)) ep = hksi(i+1,j)*tmp(i+1,j)/(ak1(i+1,j)+ak1(i+2,j)) pksi(i+1,j) = (f(i,j,nf)*ep+f(i+1,j,nf)*em)/(em+ep) end do pksi(l1,j) = 2.0*f(l1-1,j,nf)-pksi(l1-1,j) pksi(2,j) = 2.0*f(2,j,nf)-pksi(3,j) f(1,j,nf) = pksi(2,j) f(l1,j,nf) = pksi(l1,j) end do do i = 2,l1-1 do j = 2,m1-2 bpj = f(i,j,3) bpj1 = f(i,j+1,3)-hetap(i,j)-hetap(i,j+1) em = heta(i,j)*tmp(i,j)/(ae1(i,j)+ae1(i,j+1)) ep = heta(i,j+1)*tmp(i,j+1)/(ae1(i,j+1)+ae1(i,j+2)) peta(i,j+1) = (f(i,j,nf)*ep+f(i,j+1,nf)*em)/(em+ep) end do peta(i,m1) = 2.0*f(i,m1-1,nf)-peta(i,m1-1) peta(i,2) = 2.0*f(i,2,nf)-peta(i,3) f(i,1,nf) = peta(i,2) f(i,m1,nf) = peta(i,m1) end do do i = 2,l1-1 do j = 2,m1-1 xr = 0.5*(xc(i+1,j+1)+xc(i+1,j)) yr = 0.5*(yc(i+1,j+1)+yc(i+1,j)) xl = 0.5*(xc(i,j+1)+xc(i,j)) yl = 0.5*(yc(i,j+1)+yc(i,j)) xu = 0.5*(xc(i+1,j+1)+xc(i,j+1)) yu = 0.5*(yc(i+1,j+1)+yc(i,j+1)) xd = 0.5*(xc(i+1,j)+xc(i,j)) yd = 0.5*(yc(i+1,j)+yc(i,j)) dxdksi = xr-xl dydksi = yr-yl dxdeta = xu-xd dydeta = yu-yd dpksi = pksi(i+1,j)-pksi(i,j) dpeta = peta(i,j+1)-peta(i,j) t1 = -dydeta*dpksi+dydksi*dpeta t2 = dxdeta*dpksi-dxdksi*dpeta if ( mode == 1 ) then t1 = t1*r(i,j) t2 = t2*r(i,j) end if f(i,j,1) = f(i,j,1)+t1/tmp(i,j)/rho(i,j)/vol(i,j) f(i,j,2) = f(i,j,2)+t2/tmp(i,j)/rho(i,j)/vol(i,j) end do end do ! ! N > 4, T, W, TK, TE, E, PSI ! else call gamsor(ap,birad,ce1,ce2,cmu,con,ewall,f,fcsl,fksl, & fma,fo,frsl,gam,gamt,grash,hamag,heta,hksi,icrys,inturb, & izone,jcrys,l1,lsolve,m1,mode,nf,pr,r,rdtm,re,recb,rect, & rho,rhocon,rpr,sige,sigk,sigt,stel,stes,tal,tas,tf,tw,x, & xc,y,yc) do j = 2,m1-1 diff = 2.0*gam(1,j)/hksi(2,j) flow = -ruksi(2,j) call diflow(acof,diff,flow) aim(2,j) = acof*ak1(2,j) diff = 2.0*gam(l1,j)/hksi(l1-1,j) flow = ruksi(l1,j) call diflow(acof,diff,flow) aip(l1-1,j) = acof*ak1(l1,j) do i = 2,l1-2 diff = gam(i,j)*gam(i+1,j)/(0.5*hksi(i,j)*gam(i+1,j) & +0.5*hksi(i+1,j)*gam(i,j)) flow = ruksi(i+1,j) call diflow(acof,diff,flow) aip(i,j) = acof*ak1(i+1,j) aim(i+1,j) = (acof+ruksi(i+1,j))*ak1(i+1,j) end do end do do i = 2,l1-1 diff = 2.0*gam(i,1)/heta(i,2) flow = -rueta(i,2) call diflow(acof,diff,flow) ajm(i,2) = acof*ae1(i,2) diff = 2.0*gam(i,m1)/heta(i,m1-1) flow = rueta(i,m1) call diflow(acof,diff,flow) ajp(i,m1-1) = acof*ae1(i,m1) do j = 2,m1-2 diff = gam(i,j)*gam(i,j+1)/(0.5*heta(i,j)*gam(i,j+1) & +0.5*heta(i,j+1)*gam(i,j)) flow = rueta(i,j+1) call diflow(acof,diff,flow) ajp(i,j) = acof*ae1(i,j+1) ajm(i,j+1) = (acof+rueta(i,j+1))*ae1(i,j+1) end do end do do i = 2,l1-1 do j = 2,m1-1 con(i,j) = con(i,j)*vol(i,j) ap(i,j) = -ap(i,j)*vol(i,j)+aip(i,j)+aim(i,j)+ajp(i,j)+ajm(i,j) end do end do isol = 1 call flux(ae1,ae2,aim,aip,ajm,ajp,ak1,ak2,ap,b1jbl,b2jbl, & cappa,cd,cmu,con,epsil,f,fjeta,fjksi,fmax,fn,fo,gam,heta,hksi, & icrys,isol,iter,izone,jcrys,l1,lblk,lconv,lortho, & m1,mode,nf,nsolve,ntimes,r,relax,res,rueta,ruksi) end if end if end do return end subroutine setx ( l, m, ni, nj, x, xc, y, yc ) ! !*********************************************************************** ! !! SETX locates the primary nodes from the corner nodes. ! ! ! Discussion: ! ! SETX is given (XC,YC), the locations of the "corners" of the ! control volumes, and calculates the locations of the primary nodes ! (X,Y). ! ! (XC,YC) sits in the center of the control volume. On the other ! hand, the point (X,Y) is NOT necessarily the center of the four ! nearby values of XC and YC. ! ! SETX may be called to do a portion of the region, or the ! full region (in which case L = L0, M=M0). ! implicit none ! integer ni integer nj ! integer i integer j integer l integer m real x(ni,nj) real xc(ni,nj) real y(ni,nj) real yc(ni,nj) ! ! The central set of values. ! do i = 2,l-1 do j = 2,m-1 x(i,j) = 0.25*(xc(i,j)+xc(i+1,j)+xc(i,j+1)+xc(i+1,j+1)) y(i,j) = 0.25*(yc(i,j)+yc(i+1,j)+yc(i,j+1)+yc(i+1,j+1)) end do end do ! ! The first and last columns, I = 2 to L-1, J=1 and J=M. ! do i = 2,l-1 x(i,1) = 0.5*(xc(i,2)+xc(i+1,2)) y(i,1) = 0.5*(yc(i,2)+yc(i+1,2)) x(i,m) = 0.5*(xc(i,m)+xc(i+1,m)) y(i,m) = 0.5*(yc(i,m)+yc(i+1,m)) end do ! ! The first and last rows, I = 1 and I=L, J=2 to M-1. ! do j = 2,m-1 x(1,j) = 0.5*(xc(2,j)+xc(2,j+1)) y(1,j) = 0.5*(yc(2,j)+yc(2,j+1)) x(l,j) = 0.5*(xc(l,j)+xc(l,j+1)) y(l,j) = 0.5*(yc(l,j)+yc(l,j+1)) end do ! ! The corners. ! x(1,1) = xc(2,2) y(1,1) = yc(2,2) x(1,m) = xc(2,m) y(1,m) = yc(2,m) x(l,1) = xc(l,2) y(l,1) = yc(l,2) x(l,m) = xc(l,m) y(l,m) = yc(l,m) return end subroutine solve1 ( aim, aip, ajm, ajp, ap, cappa, cd, cmu, con, f, fo, & heta, hksi, icrys, izone, jcrys, l1, m1, nf ) ! !*********************************************************************** ! !! SOLVE1 sets or modifies the linear system to be handled by SOLVE2. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! real aim(ni,nj) real aip(ni,nj) real ajm(ni,nj) real ajp(ni,nj) real ap(ni,nj) real cappa real cd real cm4 real cmu real con(ni,nj) real f(ni,nj,ns) real fo(ni,nj,ns) real heta(ni,nj) real hksi(ni,nj) integer i integer icrys integer izone integer j integer jcrys integer l1 integer m1 integer nf ! ! Temperature calculation in solid zone. ! if ( izone == 1 .and. nf == 5 ) then do i = 1,icrys do j = 1,m1 ap(i,j) = 1.0 aip(i,j) = 0.0 aim(i,j) = 0.0 ajp(i,j) = 0.0 ajm(i,j) = 0.0 con(i,j) = fo(i,j,5) end do end do end if ! ! Magnetic stream function calculation in solid zone. ! if ( izone == 1 .and. nf == 9 ) then do i = 1,icrys-1 do j = 1,m1 ap(i,j) = 1.0 aip(i,j) = 0.0 aim(i,j) = 0.0 ajp(i,j) = 0.0 ajm(i,j) = 0.0 con(i,j) = fo(i,j,9) end do end do end if ! ! U, V, P, PC, or T calculation in liquid zone. ! if ( izone == 2 .and. nf <= 5 ) then do j = 2,jcrys ap(icrys,j) = 1.0 aip(icrys,j) = 0.0 aim(icrys,j) = 0.0 ajp(icrys,j) = 0.0 ajm(icrys,j) = 0.0 con(icrys,j) = 0.0 end do end if ! ! Turbulent dissipation (TE) calculation. ! if ( nf == 8 ) then cm4 = cmu**0.25 do i = 2,l1-1 con(i,2) = cd*f(i,2,7)**1.5/(cappa*cm4*0.5*heta(i,2)) ap(i,2) = 1.0 con(i,m1-1) = cd*f(i,m1-1,7)**1.5/(cappa*cm4*0.5*heta(i,m1-1)) ap(i,m1-1) = 1.0 end do do j = 2,m1-1 con(2,j) = cd*f(2,j,7)**1.5/(cappa*cm4*0.5*hksi(2,j)) ap(2,j) = 1.0 con(l1-1,j) = cd*f(l1-1,j,7)**1.5/(cappa*cm4*0.5*hksi(l1-1,j)) ap(l1-1,j) = 1.0 end do end if return end subroutine solve2 ( aim, aip, ajm, ajp, ap, con, f, l1, lblk, m1, & nf, nsolve ) ! !*********************************************************************** ! !! SOLVE2 is the tridiagonal matrix solver. ! ! ! Discussion: ! ! SOLVE2 must solve a set of linear equations defined on a five point ! stencil, of the form: ! ! A(I,J)* U(I,J) ! - A(I-1,J)*U(I-1,J) - A(I+1,J)*U(I+1,J) ! - A(I,J-1)*U(I,J-1) - A(I,J+1)*U(I,J+1) = CON. ! ! The equations must be solved for indices 2 < = I < L1-1 ! and 2 <= J <= M1-1. ! ! If I = 1 or I=L1 or J=1 or J=M1, then U(I,J) already contains its ! correct value, and these correct values must be included in the ! above equations. ! ! The solution of the equations is stored in F(*,*,NF). ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! real aim(ni,nj) real aip(ni,nj) real ajm(ni,nj) real ajp(ni,nj) real ap(ni,nj) real bl real blc real blm real blp real con(ni,nj) real denom real f(ni,nj,ns) integer i integer j integer l1 logical lblk(ns) integer m1 integer nf integer nsolve(ns) integer nt real pt(nmaxij) real qt(nmaxij) ! do nt = 1,nsolve(nf) if ( lblk(nf) ) then pt(1) = 0.0E+00 qt(1) = 0.0E+00 do i = 2,l1-1 bl = 0.0E+00 blp = 0.0E+00 blm = 0.0E+00 blc = 0.0E+00 do j = 2,m1-1 bl = bl+ap(i,j) if ( j /= m1-1) bl = bl-ajp(i,j) if ( j /= 2) bl = bl-ajm(i,j) blp = blp+aip(i,j) blm = blm+aim(i,j) blc = blc+con(i,j)+aip(i,j)*f(i+1,j,nf)+aim(i,j)*f(i-1,j,nf) & +ajp(i,j)*f(i,j+1,nf)+ajm(i,j)*f(i,j-1,nf)-ap(i,j)*f(i,j,nf) end do denom = bl-pt(i-1)*blm if ( abs(denom/bl) < 1.0E-10 ) then pt(i) = 0.0E+00 qt(i) = 0.0E+00 else pt(i) = blp/denom qt(i) = (blc+blm*qt(i-1))/denom end if end do bl = 0.0E+00 do i = l1-1,2,-1 bl = bl*pt(i)+qt(i) do j = 2,m1-1 f(i,j,nf) = f(i,j,nf)+bl end do end do pt(1) = 0.0E+00 qt(1) = 0.0E+00 do j = 2,m1-1 bl = 0.0E+00 blp = 0.0E+00 blm = 0.0E+00 blc = 0.0E+00 do i = 2,l1-1 bl = bl+ap(i,j) if ( i /= l1-1) bl = bl-aip(i,j) if ( i /= 2) bl = bl-aim(i,j) blp = blp+ajp(i,j) blm = blm+ajm(i,j) blc = blc+con(i,j)+aip(i,j)*f(i+1,j,nf)+aim(i,j)*f(i-1,j,nf) & +ajp(i,j)*f(i,j+1,nf)+ajm(i,j)*f(i,j-1,nf)-ap(i,j)*f(i,j,nf) end do denom = bl-pt(j-1)*blm if ( abs(denom/bl) < 1.0E-10 ) then pt(j) = 0.0E+00 qt(j) = 0.0E+00 else pt(j) = blp/denom qt(j) = (blc+blm*qt(j-1))/denom end if end do bl = 0.0E+00 do j = m1-1,2,-1 bl = bl*pt(j)+qt(j) do i = 2,l1-1 f(i,j,nf) = f(i,j,nf)+bl end do end do end if ! ! Sweep 1: For each column J increasing. ! do j = 2,m1-1 pt(1) = 0.0E+00 qt(1) = f(1,j,nf) do i = 2,l1-1 pt(i) = aip(i,j)/(ap(i,j)-pt(i-1)*aim(i,j)) qt(i) = (con(i,j)+ajp(i,j)*f(i,j+1,nf)+ajm(i,j)*f(i,j-1,nf) & +aim(i,j)*qt(i-1))/(ap(i,j)-pt(i-1)*aim(i,j)) end do do i = l1-1,2,-1 f(i,j,nf) = f(i+1,j,nf)*pt(i)+qt(i) end do end do ! ! Sweep 2: For each column J decreasing. ! do j = m1-2,2,-1 pt(1) = 0.0E+00 qt(1) = f(1,j,nf) do i = 2,l1-1 pt(i) = aip(i,j)/(ap(i,j)-pt(i-1)*aim(i,j)) qt(i) = (con(i,j)+ajp(i,j)*f(i,j+1,nf)+ajm(i,j)*f(i,j-1,nf) & +aim(i,j)*qt(i-1))/(ap(i,j)-pt(i-1)*aim(i,j)) end do do i = l1-1,2,-1 f(i,j,nf) = f(i+1,j,nf)*pt(i)+qt(i) end do end do ! ! Sweep 3: For each row I increasing. ! do i = 2,l1-1 pt(1) = 0.0E+00 qt(1) = f(i,1,nf) do j = 2,m1-1 pt(j) = ajp(i,j)/(ap(i,j)-pt(j-1)*ajm(i,j)) qt(j) = (con(i,j)+aip(i,j)*f(i+1,j,nf)+aim(i,j)*f(i-1,j,nf) & +ajm(i,j)*qt(j-1))/(ap(i,j)-pt(j-1)*ajm(i,j)) end do do j = m1-1,2,-1 f(i,j,nf) = f(i,j+1,nf)*pt(j)+qt(j) end do end do ! ! Sweep 4: For each row I decreasing. ! do i = l1-1,2,-1 pt(1) = 0.0E+00 qt(1) = f(i,1,nf) do j = 2,m1-1 pt(j) = ajp(i,j)/(ap(i,j)-pt(j-1)*ajm(i,j)) qt(j) = (con(i,j)+aip(i,j)*f(i+1,j,nf)+aim(i,j)*f(i-1,j,nf) & +ajm(i,j)*qt(j-1))/(ap(i,j)-pt(j-1)*ajm(i,j)) end do do j = m1-1,2,-1 f(i,j,nf) = f(i,j+1,nf)*pt(j)+qt(j) end do end do end do return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end