program crystal_qed ! !*********************************************************************** ! !! CRYSTAL_QED seeks parameters that minimize a crystallization cost functional. ! ! ! Discussion: ! ! The DQED package is used to solve the bounded and constrained least ! squares problem. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! implicit none ! ! Parameters that don't depend on anything else. ! integer, parameter :: liopt = 17 integer, parameter :: lropt = 1 integer, parameter :: mcon = 6 integer, parameter :: mequa = 1 integer, parameter :: mvars = 7 integer, parameter :: nt = 5 ! ! First order parameters. ! integer, parameter :: ldfj = mcon+mequa integer, parameter :: nall = mcon+2*mvars+nt+1 ! ! Second order parameters. ! integer, parameter :: liwork = 3*mcon+9*mvars+4*nt+nall+11 integer, parameter :: lwork = nall*nall+4*nall+mvars*nt+33*mvars & +mequa*nt+mvars*nt+14*nt+9*mcon+26+3*nall+2 ! double precision bl(mvars+mcon) double precision bu(mvars+mcon) double precision dnrm2 external dqedev double precision fj(ldfj,mvars+1) double precision fnorm external func double precision fx(mcon+mequa) integer i integer igo integer ind(mvars+mcon) integer iopt(liopt) integer iopti integer ivars integer iwork(liwork) integer nvars double precision ropt(lropt) real tarray(2) double precision temp real time1 real time2 double precision value double precision work(lwork) double precision xpar(mvars) ! call timestamp ( ) ! ! Read the system CPU clock at the starting time. ! call etime ( tarray ) time1 = tarray(1)+tarray(2) ! ! Say hello. ! call hello(liwork,lwork) ! ! Tell DQED the size of the work arrays. ! iwork(1) = lwork iwork(2) = liwork ! ! Set the initial parameter values. ! ivars = 0 ! ! These are Hui Zhang's original crucible shape parameters. ! if ( ivars == 0 ) then nvars = 7 xpar(1) = 0.015D+00 xpar(2) = 0.03D+00 xpar(3) = 0.046D+00 xpar(4) = 0.07D+00 xpar(5) = 0.10D+00 xpar(6) = 0.14D+00 xpar(7) = 0.18D+00 ! ! These crucible shape parameters represent the optimum solution ! when the nodes are monotonically constrained. ! else if ( ivars == 1 ) then nvars = 7 xpar(1) = 0.064D+00 xpar(2) = 0.064D+00 xpar(3) = 0.064D+00 xpar(4) = 0.25D+00 xpar(5) = 0.25D+00 xpar(6) = 0.25D+00 xpar(7) = 0.25D+00 ! ! These crucible shape parameters cause the grid routine to fail, ! by generating degenerate control volumes. ! else if ( ivars == 2 ) then nvars = 7 value = 0.30D+00 xpar(1) = 3.0D+00*value/5.0D+00 xpar(2) = value xpar(3) = value-(value-0.25D+00)/5.0D+00 xpar(4) = value-2.0D+00*(value-0.25D+00)/5.0D+00 xpar(5) = value-3.0D+00*(value-0.25D+00)/5.0D+00 xpar(6) = value-4.0D+00*(value-0.25D+00)/5.0D+00 xpar(7) = value-4.5D+00*(value-0.25D+00)/5.0D+00 end if write ( *, '(a)' ) ' ' write ( *, * ) 'Initial parameter values:' write ( *, '(a)' ) ' ' do i = 1,nvars write ( *, * ) i,xpar(i) end do ! ! Set the constraints on the variables. ! My first guess is to just restrain all the variables to be ! between 0 and 0.25. ! if ( ivars == 0 .or. ivars == 1 .or. ivars == 2 ) then ind(1) = 3 bl(1) = 0.0D+00 bu(1) = 0.25D+00 ind(2) = 3 bl(2) = 0.0D+00 bu(2) = 0.25D+00 ind(3) = 3 bl(3) = 0.0D+00 bu(3) = 0.25D+00 ind(4) = 3 bl(4) = 0.0D+00 bu(4) = 0.25D+00 ind(5) = 3 bl(5) = 0.0D+00 bu(5) = 0.25D+00 ind(6) = 3 bl(6) = 0.0D+00 bu(6) = 0.25D+00 ind(7) = 3 bl(7) = 0.0D+00 bu(7) = 0.25D+00 ! ! The following 6 constraints force monotonicity. ! ind(8) = 1 bl(8) = 0.0D+00 bu(8) = 0.0D+00 ind(9) = 1 bl(9) = 0.0D+00 bu(9) = 0.0D+00 ind(10) = 1 bl(10) = 0.0D+00 bu(10) = 0.0D+00 ind(11) = 1 bl(11) = 0.0D+00 bu(11) = 0.0D+00 ind(12) = 1 bl(12) = 0.0D+00 bu(12) = 0.0D+00 ind(13) = 1 bl(13) = 0.0D+00 bu(13) = 0.0D+00 end if ! ! Set IOPTI, the switch which chooses ! just one function call (for checking), ! just one jacobian call (for checking), ! or optimization. ! iopti = 0 write ( *, '(a)' ) ' ' if ( iopti == 0 ) then write ( *, '(a)' ) 'Just call the function evaluator once.' else if ( iopti == 1 ) then write ( *, '(a)' ) 'Just call the jacobian approximator once.' else write ( *, '(a)' ) 'Optimize the problem.' end if ! ! If IOPTI = 0, just call the function evaluator once. ! if ( iopti == 0 ) then call func(fx,iopt,mcon,mequa,nvars,ropt,xpar) fnorm = dnrm2(mequa,fx(mcon+1),1) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) 'Two-norm of residual is ', fnorm ! ! If IOPTI = 1, just call the jacobian approximator once. ! else if ( iopti == 1 ) then call diffor(fj,func,fx,iopt,ldfj,mcon,mequa,nvars,ropt,xpar) ! ! If IOPTI = 2, then carry out an optimization. ! else igo = 0 ! ! IOPT(1) = 2 means change the value of ITMAX to IOPT(2). ! iopt(1) = 2 iopt(2) = 10 write ( *, '(a)' ) ' ' write ( *, * ) 'Maximum number of minimization steps is ',iopt(2) ! ! IOPT(3) = 99 means no more options. ! iopt(3) = 99 call dqed(dqedev,mequa,nvars,mcon,ind,bl,bu,xpar,fj,ldfj, & fnorm,igo,iopt,ropt,iwork,work) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'DQED return flag IGO = ', igo write ( *, '(a,i6)' ) 'DQED real work array requirement = ', iwork(1) write ( *, '(a,i6)' ) 'DQED integer work array requirement = ', iwork(2) write ( *, '(a,g14.6)' ) 'Two-norm of residual is ', fnorm write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQED returned the parameter values:' write ( *, '(a)' ) ' ' do i = 1,nvars write ( *, '(i6,g14.6)' ) i,xpar(i) end do write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) 'with the function value ', fnorm call func(fx,iopt,mcon,mequa,nvars,ropt,xpar) fnorm = dnrm2(mequa,fx(mcon+1),1) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I check the output function value!' write ( *, '(a,g14.6)' ) 'I get a function value ', fnorm end if ! ! Having completed the desired task, read the system clock and stop. ! call etime(tarray) time2 = tarray(1)+tarray(2) write ( *, '(a)' ) 'CRYSTAL_QED:' temp = time2-time1 write ( *, '(a,g14.6,a)' ) ' Total elapsed CPU time = ',temp,' seconds,' temp = temp/60.0 write ( *, '(a,g14.6,a)' ) ' = ',temp,' minutes.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL_QED:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine dqedev ( xpar, fj, ldfj, igo, iopt, ropt ) ! !*********************************************************************** ! !! DQEDEV evaluates certain functions being treated by DQED. ! ! ! Discussion: ! ! DQEDEV also evaluates partial derivatives. ! ! The user has NVARS variables XPAR(I), and is trying to minimize ! the square root of the sum of the squares of MEQUA functions ! F(I)(XPAR), subject to MCON constraints which have the form ! ! BL(I) < = G(I)(XPAR) <= BU(I) ! ! where either the left or right bounding inequality may be dropped. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision XPAR(*). ! XPAR is an array of length NVARS, containing the values of the ! independent variables at which the functions and partial ! derivatives should be evaluated. ! ! Output, double precision FJ(LDFJ,NVARS+1). ! ! If IGO is nonzero, then partial derivatives must ! be placed in the first NVARS columns of FJ, as ! follows: ! ! Rows I = 1 to MCON, and columns J = 1 to NVARS ! should contain the partials of the I-th constraint ! function G(I) with respect to the J-th variable. ! ! Rows I = MCON+1 to MCON+MEQUA, and columns J = 1 to NVARS ! should contain the partials of the (I-MCON)-th nonlinear ! function F(I-MCON) with respect to the J-th variable. ! ! Regardless of the value of IGO, column NVARS+1 must be ! assigned the values of the constraints and functions, ! as follows: ! ! Rows I = 1 to MCON, column J = NVARS+1, should contain ! the value of G(I). ! ! Rows I = MCON+1 to MCON+MEQUA, column J = NVARS+1, should ! contain the value of F(I-MCON). ! ! Input, integer LDFJ. ! LDFJ is the leading dimension of FJ, which must be at least ! MCON+MEQUA. ! ! Input/output, integer IGO. ! On input, IGO tells the user whether the partial derivatives ! are needed. ! ! 0, the partial derivatives are not needed. ! nonzero, the partial derivatives are needed. ! ! On output, the user may reset the input value of IGO if one ! of two situations is encountered: ! ! 99, the functions, constraints, or partial derivatives ! could not be evaluated at the given input point X. Request ! that DQED reject that point, and try a different one. ! ! Any other value, abort the run. ! ! Input, integer IOPT(*), the integer option array. ! ! Input, double precision ROPT(*), the real option array. ! implicit none ! integer ldfj integer, parameter :: mcon = 6 integer, parameter :: mequa = 1 integer, parameter :: nvars = 7 ! double precision fj(ldfj,nvars+1) external func double precision fx(mcon+mequa) integer igo integer iopt(*) double precision ropt(*) double precision xpar(nvars) ! if ( igo /= 0 ) then call diffor(fj,func,fx,iopt,ldfj,mcon,mequa,nvars,ropt,xpar) end if call func(fj(1,nvars+1),iopt,mcon,mequa,nvars,ropt,xpar) return end subroutine func ( fx, iopt, mcon, mequa, nvars, ropt, xpar ) ! !*********************************************************************** ! !! FUNC communicates between the optimizer and the constitutive routines. ! ! ! Discussion: ! ! FUNC receives a set of parameter values as input, solves the ! resulting constitutive equations, and evaluates the derived ! cost functional, whose value is returned to the calling routine. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IOPT(*), the integer option array. ! ! Input, integer MCON. ! MCON is the number of constraints imposed on the independent ! variables that compose a feasible solution. ! ! Input, integer MEQUA. ! MEQUA is the number of components of the nonlinear function. ! ! Input, integer NVARS. ! NVARS is the number of parameters or independent variables ! upon which the function depends. ! ! Input, double precision ROPT(*). ! ROPT is the real option array. ! ! Input, double precision XPAR(*). ! XPAR is an array of length NVARS, containing the values of the ! parameters. ! implicit none ! integer, parameter :: maxbot = 20 integer mcon integer mequa integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: ns = 10 integer nvars ! double precision ae1(ni,nj) double precision ae2(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision area(ni,nj) double precision areal double precision areas double precision areat double precision b double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision birad double precision bo double precision cappa double precision cd double precision ce1 double precision ce2 double precision cfo double precision cinc double precision cinco double precision cmu double precision cost double precision cvn double precision damax double precision delt double precision dtm double precision epsad double precision epsil double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fjeta(ni,nj) double precision fjksi(ni,nj) double precision fks double precision fksl double precision fma double precision fmax(ns) double precision fn(ni,nj,ns) double precision fnu double precision fo(ni,nj,ns) double precision fr double precision frsl double precision fx(mcon+mequa) double precision gam(ni,nj) double precision gamt(ni,nj) double precision grash double precision hamag double precision heta(ni,nj) double precision hf double precision hksi(ni,nj) integer i integer icost integer icrys integer inturb integer iopt(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 np integer npc integer nsolve(ns) integer ntimes(ns) double precision orth double precision pr double precision r(ni,nj) double precision ra double precision rdtm double precision re double precision recb double precision rect double precision relax(nk) double precision res(ns) double precision rho(ni,nj) double precision rhocon double precision ropt(*) double precision rpr double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision sige double precision sigk double precision sigma double precision sigt double precision smax double precision smooth double precision ssum double precision stel double precision stes double precision t1 double precision t2 double precision t3 double precision tal double precision tanca double precision tanca2 double precision tas double precision tend double precision tf double precision tinit character ( len = 25 ) title(ns) double precision tnow double precision tw double precision vave double precision vol(ni,nj) double precision vort(ni,nj) double precision x(ni,nj) double precision xbot(maxbot) double precision xc(ni,nj) double precision xlen double precision xpar(nvars) double precision y(ni,nj) double precision ybot(maxbot) double precision yc(ni,nj) double precision ylen ! ncall = ncall+1 ! ! Initialize the data. ! call inidat(ae1,ae2,ak1,ak2,area,b,b1jbl,b2jbl,b3jbl, & 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,gam,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,npc,ns,nsolve,ntimes,orth, & 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,vort,x,xbot,xc,xlen,y,ybot,yc,ylen) ! ! TEMPORARY ! do i = 1,7 xbot(i+1) = xpar(i) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNC:' write ( *, '(a)' ) ' Input parameters:' do i = 1,7 write ( *, '(i6,g14.6)' ) i,xpar(i) end do ! ! 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 setp(fjeta,fjksi,heta,hksi,l1,m1,ni,nj,f(1,1,3)) call setru(heta,hksi,l1,m1,ni,nj,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,b3jbl,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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNC - Warning!' write ( *, '(a)' ) ' The solid iteration has not converged' write ( *, '(a,i6,a)' ) ' after ',iter,' iterations.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Solid norm and max relative change:' write ( *, '(a)' ) ' ' do i = 1,ns if ( lsolve(i) ) then fmax(i) = damax ( l0*m0, f(1,1,i), 1 ) 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 do i = 1,icrys ruksi(i,1) = 0.0D+00 ruksi(i,m0) = 0.0D+00 rueta(i,2) = 0.0D+00 rueta(i,m0) = 0.0D+00 end do do j = 1,m0 rueta(1,j) = 0.0D+00 rueta(icrys,j) = 0.0D+00 ruksi(2,j) = 0.0D+00 ruksi(icrys,j) = 0.0D+00 end do 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.0D+00 end do end do else if ( inturb == 1 ) then do i = 2,l1-1 do j = 2,m1-1 gamt(i,j) = (1.0D+00-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,b3jbl,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.0D+00 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.5D+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.5D+00*(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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNC - Warning!' write ( *, '(a)' ) ' The liquid iteration has not converged' write ( *, '(a,i6,a)' ) ' after ',iter,' iterations.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Liquid norms and max relative change:' write ( *, '(a)' ) ' ' do i = 1,ns if ( lsolve(i) ) then fmax(i) = damax(l0*m0,f(1,1,i),1) 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 ! ! Compute the vorticity VORT. ! call vortic ( heta, hksi, icrys, l0, m0, f(1,1,1), f(1,1,2), vort, xc, yc ) ! ! Evaluate the cost function integrand at the current time. ! cinco = cinc call setcst(area,cinc,gam,icost,icrys,jcrys,l0,m0,ni,nj,f(1,1,5), & tf,f(1,1,1),f(1,1,2),vave,vort,xc,yc) ! ! Update the total cost, by adding the estimated contribution ! from the current time interval. ! ! TEMPORARY ! if ( ndt > 0 ) then ! cost = cost+0.5D+00*delt*(cinco+cinc) cost = cost+0.5D+00*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) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HEY, ABOUT TO CALL RSWRIT.' call rswrit(b1jbl,b2jbl,b3jbl,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),vort,f(1,1,6),x,xbot,xc,y,ybot,yc) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'HEY, ABOUT TO CALL WRTEC.' call wrtec(gamt,l1,m1,f(1,1,3),f(1,1,10),f(1,1,5),f(1,1,8), & f(1,1,7),f(1,1,1),f(1,1,2),f(1,1,6),x,y) end if write ( *, '(a,g14.6)' ) 'FUNC: COST = ',cost do i = 1,mcon fx(i) = xpar(i+1)-xpar(i) end do fx(mcon+1) = 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. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision a1 double precision a11 double precision a12 double precision a2 double precision a21 double precision a22 double precision a3 double precision b1 double precision b2 double precision b3 double precision c1 double precision c2 double precision c3 double precision cvn double precision det double precision dot double precision dx double precision dxc(ni,nj) double precision dxmax double precision dy double precision dyc(ni,nj) double precision epsad double precision gn double precision gt double precision gx double precision gy integer i integer icrys integer iprint integer j integer jcrys integer l1 integer m1 integer nin double precision orth double precision pb(ni,nj) double precision pc1 double precision pc2 double precision pd(ni,nj) double precision ratio double precision res1 double precision res2 double precision rmax double precision s1 double precision s2 double precision scale double precision smooth double precision ssmot double precision ssweg double precision vmag2 double precision wn double precision wt double precision wx double precision wy double precision xbot(nbot) double precision xc(ni,nj) double precision xdisp double precision xij double precision xjac double precision xn double precision xnn double precision xpp double precision xt double precision xtn double precision xtt double precision ybot(nbot) double precision yc(ni,nj) double precision ydisp double precision yij double precision yn double precision ynn double precision ypp double precision yt double precision ytn double precision ytt ! dxmax = 0.0D+00 do i = 1,l1 do j = 1,m1 dxc(i,j) = 0.0D+00 dyc(i,j) = 0.0D+00 end do end do do i = 2,l1 do j = 2,m1 if ( i < icrys ) then pc1 = 0.95D+00*((i-(2+icrys)/2)/10.0D+00)**2+0.05D+00 else pc1 = 0.8D+00*((i-l1)/12.0)**2+0.2D+00 end if if ( j <= jcrys ) then pc2 = 0.9D+00*((j-(2+jcrys)/2)/10.0D+00)**2+0.05D+00 else pc2 = 0.9D+00*((j-(m1+jcrys)/2)/10.0D+00)**2+0.05D+00 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.0D+00 ssweg = 0.0D+00 rmax = 0.0D+00 if ( iprint > 0 ) then if ( nin == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) '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.5D+00*(xc(i+1,j)-xc(i-1,j)) xn = 0.5D+00*(xc(i,j+1)-xc(i,j-1)) yt = 0.5D+00*(yc(i+1,j)-yc(i-1,j)) yn = 0.5D+00*(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.0D+00*smooth*(xt*yt+xn*yn)*(xn*xn+yn*yn)/xjac**3 a2 = 4.0D+00*smooth*(xt*yt+xn*yn)*(xt*xn+yt*yn)/xjac**3 a3 = -2.0D+00*smooth*(xt*yt+xn*yn)*(xt*xt+yt*yt)/xjac**3 b1 = 2.0D+00*smooth*(yt*yt+yn*yn)*(xn*xn+yn*yn)/xjac**3 b2 = -4.0D+00*smooth*(yt*yt+yn*yn)*(xt*xn+yt*yn)/xjac**3 b3 = 2.0D+00*smooth*(yt*yt+yn*yn)*(xt*xt+yt*yt)/xjac**3 c1 = 2.0D+00*smooth*(xt*xt+xn*xn)*(xn*xn+yn*yn)/xjac**3 c2 = -4.0D+00*smooth*(xt*xt+xn*xn)*(xt*xn+yt*yn)/xjac**3 c3 = 2.0D+00*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.0D+00*orth*pd(i,j)*(2.0D+00*xt*xn+yt*yn) b3 = b3+orth*pd(i,j)*xt*xt c1 = c1+orth*pd(i,j)*yn*yn c2 = c2+2.0D+00*orth*pd(i,j)*(xt*xn+2.0D+00*yt*yn) c3 = c3+orth*pd(i,j)*yt*yt a1 = a1-2.0D+00*cvn*pb(i,j)*xn*yn a2 = a2+2.0D+00*cvn*pb(i,j)*(xt*yn+xn*yt) a3 = a3-2.0D+00*cvn*pb(i,j)*xt*yt b1 = b1+2.0D+00*cvn*pb(i,j)*yn*yn b2 = b2-4.0D+00*cvn*pb(i,j)*yt*yn b3 = b3+2.0D+00*cvn*pb(i,j)*yt*yt c1 = c1+2.0D+00*cvn*pb(i,j)*xn*xn c2 = c2-4.0D+00*cvn*pb(i,j)*xt*xn c3 = c3+2.0D+00*cvn*pb(i,j)*xt*xt ! ! Now compute terms from the derivatives of the weight functions. ! gt = 0.5D+00*(pd(i+1,j)-pd(i-1,j)) gn = 0.5D+00*(pd(i,j+1)-pd(i,j-1)) gx = orth*0.5D+00*(gt*yn-gn*yt)*(xt*xn+yt*yn)**2/xjac gy = orth*0.5D+00*(gn*xt-gt*xn)*(xt*xn+yt*yn)**2/xjac wt = 0.5D+00*(pb(i+1,j)-pb(i-1,j)) wn = 0.5D+00*(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.0D+00*xc(i,j)+xc(i-1,j) xnn = xc(i,j+1)-2.0D+00*xc(i,j)+xc(i,j-1) ytt = yc(i+1,j)-2.0D+00*yc(i,j)+yc(i-1,j) ynn = yc(i,j+1)-2.0D+00*yc(i,j)+yc(i,j-1) xtn = 0.25D+00*(xc(i+1,j+1)+xc(i-1,j-1)-xc(i-1,j+1)-xc(i+1,j-1)) ytn = 0.25D+00*(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.0D+00*(b1+b3) a12 = -2.0D+00*(a1+a3) a22 = -2.0D+00*(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.25D+00 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.5D+00,0.25D+00/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.4D+00+0.02D+00*sin(xc(i,jcrys)*20.0D+00*3.14159D+00) 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.25D+00 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.3D+00,0.25D+00/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.002D+00*dble(i-2)) ) then xc(i,m1) = xc(2,m1)+0.002D+00*dble(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: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NBOT, the number of points to be used to ! define the bottom of the crucible. ! implicit none ! integer nbot integer, parameter :: ni = 64 ! double precision dome integer i integer ist double precision p2 double precision p3 double precision shift double precision t1 double precision t2 double precision tm(ni) double precision tt(ni) double precision x1 double precision x2 double precision xbot(nbot) double precision xl(ni) double precision xnew double precision y1 double precision y2 double precision ybot(nbot) double precision yl(ni) double precision ymin double precision ynew ! if ( ynew < ybot(1) .or. ybot(nbot).lt.ynew ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CUBIC - Fatal error!' write ( *, '(a,g14.6)' ) ' YNEW = ',ynew write ( *, '(a,g14.6)' ) ' outside of range YBOT(1) = ',ybot(1) write ( *, '(a,g14.6)' ) ' through YBOT(NBOT) = ',ybot(nbot) stop end if do i = 1,ni xl(i) = 0.0D+00 yl(i) = 0.0D+00 tt(i) = 0.0D+00 tm(i) = 0.0D+00 end do 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.0D-10) then tt(i) = 0.0D+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.0D+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.0D+00*(x2-x1)/(y2-y1)-2.0D+00*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 function damax ( n, dx, incx ) ! !*********************************************************************** ! !! DAMAX returns the maximum absolute value of the entries of a vector. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries to examine. ! ! Input, double precision DX(*), the vector whose entries ! are to be examined. ! ! Input, integer INCX, the increment between successive entries ! of DX to be searched. ! ! Output, double precision DAMAX, the value of the maximum ! absolute value of the examined entries of DX. ! implicit none ! integer n ! double precision damax double precision dx(*) integer i integer incx integer j ! damax = 0.0D+00 i = 1 j = 1 10 continue if ( i <= n ) then if ( abs(dx(j)) > damax ) then damax = abs(dx(j)) end if i = i+1 j = j+incx go to 10 end if return end subroutine diflow ( acof, diff, flow ) ! !*********************************************************************** ! !! DIFLOW computes the convection-diffusion coefficient ACOF. ! ! ! Discussion: ! ! DIFLOW is 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. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! double precision acof double precision diff double precision flow ! acof = diff if ( diff == 0.0D+00) return ! ! 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.0D+00-0.1D+00*abs(flow/diff))**5) if ( flow < 0.0D+00) acof = acof-flow return end subroutine doarea(area,areal,areas,areat,icrys,jcrys,l0,m0,ni,nj,xc,yc) ! !*********************************************************************** ! !! DOAREA calculates control volume areas. ! ! ! 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 which is bounded ! by the N points (X(I),Y(I)) 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. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision AREA(NI,NJ). ! AREA contains the area of the control volume (I,J) which is ! bounded by the corner nodes [I,J], [I+1,J], [I,J+1], and [I+1,J+1]. ! ! Output, double precision AREAL, the area of the liquid region. ! ! Output, double precision AREAS. ! AREAS is the area of the solid or crystal region. ! ! Output, double precision AREAT, the total area. ! ! Input, integer ICRYS. ! ICRYS specifies the end of the crystal in the I array direction, ! and in the vertical coordinate direction. ! ! Input, integer JCRYS. ! JCRYS specifies the end of the crystal in the J array direction, ! and in the horizontal coordinate direction. ! ! Input, integer L0. ! L0 is the extent of the grid in the vertical or "I" coordinate. ! ! Input, integer M0. ! M0 is the extent of the grid in the horizontal or "J" coordinate. ! ! Input, integer NI. ! NI is the maximum number of grid points in the "I" or vertical ! direction. ! ! Input, integer NJ. ! NJ is the maximum number of grid points in the "J" or horizontal ! direction. ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer ni integer nj ! double precision area(ni,nj) double precision areal double precision areas double precision areat integer i integer icrys integer j integer jcrys integer l0 integer m0 integer nbad double precision xc(ni,nj) double precision 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.0D+00 else area(i,j) = 0.5D+00*( & 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & '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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & '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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & '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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,2i6)' ) ' Zero length side for cell I,J:',i,j write ( *, '(a,2g14.6)' ) & '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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,i6,a)' ) ' 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.0D+00 ) then nbad = nbad+1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,i6,a,i6)' ) ' Zero area for cell I = ',i,' J=',j end if end do end do if ( nbad > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SETARE - Fatal error!' write ( *, '(a,i6,a)' ) ' 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 = 0.0 do i = 1,l0 do j = 1,m0 areat = areat+area(i,j) end do end do return end subroutine findp ( j, xold, yold, xnew, ynew, xc, yc ) ! !*********************************************************************** ! !! FINDP finds the boundary nodes associated with a Neumann condition. ! ! ! Discussion: ! ! FINDP is given a point (XOLD,YOLD). ! ! FINDP returns a pair of values (XNEW,YNEW), which represent ??? ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision a1 double precision ai double precision aj double precision alp double precision b1 double precision bet double precision c1 double precision dels double precision denm double precision dxb1 double precision dxb2 double precision dyb1 double precision dyb2 integer j double precision r0 double precision r1 double precision r11 double precision r2 double precision r22 double precision r33 double precision slpi double precision xc(ni,nj) double precision xnew double precision xnew1 double precision xnew2 double precision xold double precision yc(ni,nj) double precision yn1 double precision yn2 double precision ynew double precision 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.0001D+00 ) then if ( abs(dxb1) >= 0.001D+00 ) 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.5D+00*(r22-r11) r2 = 0.5D+00*(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.0D+00+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, & b3jbl,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 flux of various quantities. ! ! ! It looks like if ISOL is 0, you just go through the big loop once. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! double precision acof double precision ae1(ni,nj) double precision ae2(ni,nj) double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision ap(ni,nj) double precision ap0(ni,nj) double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision cappa double precision cd double precision cmu double precision con(ni,nj) double precision con0(ni,nj) double precision diff double precision epsil double precision error double precision f(ni,nj,ns) double precision fjbi1(nj,ns) double precision fjbj1(ni,ns) double precision fjbl1(nj,ns) double precision fjbm1(ni,ns) double precision fjeta(ni,nj) double precision fjksi(ni,nj) double precision flow double precision fmax(ns) double precision fn(ni,nj,ns) double precision fo(ni,nj,ns) double precision gam(ni,nj) double precision heta(ni,nj) double precision 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) double precision pt(nmaxij) double precision qt(nmaxij) double precision r(ni,nj) double precision relax(nk) double precision res(ns) double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision t1 double precision t2 double precision t3 double precision t4 double precision tem1 double precision tem2 double precision xje double precision xjn double precision xjw double precision xjs double precision xme double precision xmn double precision xmw double precision xms ! ! Save copies of CON and AP. ! do i = 1,l1 do j = 1,m1 con0(i,j) = con(i,j) ap0(i,j) = ap(i,j) end do end do ! ! 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.5D+00*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.5D+00*hksi(i,jj)*gam(i-1,jj) & +0.5D+00*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.5D+00*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.5D+00*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.5D+00*heta(ii,j)*gam(ii,j-1) & +0.5D+00*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.5D+00*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.0D+00 ) then fjksi(2,j) = (con(1,j)*heta(1,j)+0.5D+00*ak2(2,j) & *(fjeta(1,j)+fjeta(1,j+1)))/ak1(2,j) end if if ( gam(l1,j) == 0.0D+00 ) then fjksi(l1,j) = (con(l1,j)*heta(l1,j)+0.5D+00*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.0D+00 ) then fjeta(i,2) = (con(i,1)*hksi(i,1)+0.5D+00*ae2(i,2) & *(fjksi(i,1)+fjksi(i+1,1)))/ae1(i,2) end if if ( gam(i,m1) == 0.0D+00 ) then fjeta(i,m1) = (con(i,m1)*hksi(i,m1)+0.5D+00*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.5D+00*(fjksi(i,1)+fjksi(i+1,1)) t2 = 0.5D+00*(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.5D+00*(fjeta(1,j)+fjeta(1,j+1)) t2 = 0.5D+00*(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. ! We compute ! B1JBL, the heat flux between rows (?) ICRYS-1 and ICRYS, ! for columns 2 through JCRYS and column M1. ! B2JBL, the heat flux between rows (?) ICRYS and ICRYS+1, ! for columns 2 through JCRYS and column M1. ! B3JBL, the heat flux between rows (?) and (?), ! for columns JCRYS+1 to M1-1. ! if ( nf == 5 ) then do j = 2,jcrys-1 t1 = gam(icrys-1,j)*(f(icrys-1,j,5)-f(icrys,j,5))/hksi(icrys-1,j) 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)*(f(icrys,j,5)-f(icrys+1,j,5))/hksi(icrys+1,j) 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)*(f(icrys-1,j,5)-f(icrys,j,5))/hksi(icrys-1,j) t2 = gam(icrys-1,j)*(f(icrys,j+1,5)-f(icrys,j,5))/heta(icrys,j) 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.0D+00 ) 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.0D+00 .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.0D+00 ) 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.0D+00 .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.5D+00*(fjeta(i,j)+fjeta(i,j+1)) qt(i) = 0.5D+00*(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.5D+00*(fjksi(i,j)+fjksi(i+1,j)) qt(j) = 0.5D+00*(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.0D+00 if ( gam(1,j) == 0.0D+00 ) then t1 = fjksi(2,j)-ruksi(2,j)*f(2,j,nf) diff = gam(2,j)/(0.5D+00*hksi(2,j)) flow = -ruksi(2,j) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) 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.0D+00 if ( gam(l1,j) == 0.0D+00 ) then t1 = fjksi(l1,j)-ruksi(l1,j)*f(l1-1,j,nf) diff = gam(l1-1,j)/(0.5D+00*hksi(l1-1,j)) flow = ruksi(l1,j) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) 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.0D+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.0D+00 if ( gam(i,1) == 0.0D+00 ) then t1 = fjeta(i,2)-rueta(i,2)*f(i,2,nf) diff = gam(i,2)/(0.5D+00*heta(i,2)) flow = -rueta(i,2) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) 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.0D+00 if ( gam(i,m1) == 0.0D+00 ) then t1 = fjeta(i,m1)-rueta(i,m1)*f(i,m1-1,nf) diff = gam(i,m1-1)/(0.5D+00*heta(i,m1-1)) flow = rueta(i,m1) call diflow(acof,diff,flow) if ( acof /= 0.0D+00 ) 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.0D+00 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.0D+00-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.0D+00 do i = 1,l1 do j = 1,m1 error = abs(f(i,j,nf)-fn(i,j,nf)) if ( fmax(nf) > 0.0D+00 ) 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 ( *, '(a,i6)' ) 'FLUX - Convergence for variable ',nf ! write ( *, '(a,i6)' ) ' on step ',nt ! write ( *, '(a,g14.6)' ) ' RES(NF) = ',res(nf) ! write ( *, '(a,g14.6)' ) ' FMAX(NF)=',fmax(nf) ! return end if end do ! ! End of big iteration ! ! Compute the largest solution component. ! ! fmax(nf) = 0.0D+00 ! 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.0D+00 ! do i = 1,l1 ! do j = 1,m1 ! ! error = abs(f(i,j,nf)-fn(i,j,nf)) ! ! if ( fmax(nf) > 0.0D+00 ) 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 ( *, '(a,i6)' ) 'FLUX - No convergence for variable ',nf ! write ( *, '(a,i6)' ) ' on step ',nt ! write ( *, '(a,g14.6)' ) ' RES(NF) = ',res(nf) ! write ( *, '(a,g14.6)' ) ' 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. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! double precision ap(ni,nj) double precision ar double precision birad double precision ce1 double precision ce2 double precision cm4 double precision cmu double precision con(ni,nj) double precision dcdr double precision depdr double precision depdx double precision dudx(ni,nj) double precision dudy(ni,nj) double precision dvdx(ni,nj) double precision dvdy(ni,nj) double precision dwdx(ni,nj) double precision dwdy(ni,nj) double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fksl double precision fma double precision fo(ni,nj,ns) double precision frsl double precision gam(ni,nj) double precision gamt(ni,nj) double precision gdudx(ni,nj) double precision gdudy(ni,nj) double precision gdvdx(ni,nj) double precision gdvdy(ni,nj) double precision grash double precision hamag double precision heta(ni,nj) double precision 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 double precision p11 double precision p12 double precision p21 double precision p22 double precision p31 double precision p32 double precision p41 double precision p42 double precision pbig double precision pr double precision prod(ni,nj) double precision r(ni,nj) double precision rdtm double precision re double precision recb double precision rect double precision rho(ni,nj) double precision rhocon double precision rpr double precision sige double precision sigk double precision sigt double precision stel double precision stes double precision tal double precision tas double precision tf double precision tw double precision x(ni,nj) double precision xc(ni,nj) double precision xplus(nmaxij) double precision y(ni,nj) double precision yc(ni,nj) double precision yplus(nmaxij) ! do i = 1,l1 do j = 1,m1 ap(i,j) = -rdtm*rho(i,j) con(i,j) = 0.0D+00 gam(i,j) = rhocon/re+gamt(i,j) end do end do ! ! 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.0D+00 gdudy(i,j) = 0.0D+00 gdvdx(i,j) = 0.0D+00 gdvdy(i,j) = 0.0D+00 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.63D+00 ) then gam(i,1) = rhocon/re else gam(i,1) = yplus(i)/(2.5D+00*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.5D+00*heta(i,m1-1)/rhocon if ( yplus(i) <= 11.63D+00 ) then gam(i,m1) = rhocon/re else gam(i,m1) = yplus(i)/(2.5D+00*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.5D+00*hksi(2,j)/rhocon if ( xplus(j) <= 11.63D+00 ) then gam(1,j) = rhocon/re else gam(1,j) = xplus(j)/(2.5D+00*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.5D+00*hksi(l1-1,j)/rhocon if ( xplus(j) <= 11.63D+00 ) then gam(l1,j) = rhocon/re else gam(l1,j) = xplus(j)/(2.5D+00*log(ewall*xplus(j))) end if end do end if do i = 2,l1-1 gam(i,1) = 0.0D+00 end do do j = 2,m1-1 f(1,j,1) = 0.0D+00 f(l1,j,1) = 0.0D+00 end do f(l1,1,1) = 0.0D+00 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.0D+00 ! ! 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.0D+00 gdudy(i,j) = 0.0D+00 gdvdx(i,j) = 0.0D+00 gdvdy(i,j) = 0.0D+00 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.25 if ( mode == 0 ) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5D+00*heta(i,2)/rhocon gam(i,1) = rhocon/re if ( yplus(i) > 11.63D+00) then gam(i,1) = yplus(i)/(2.5D+00*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.63D+00 ) then gam(i,m1) = yplus(i)/(2.5D+00*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.5D+00*hksi(2,j)/rhocon gam(1,j) = rhocon/re if ( xplus(j) > 11.63D+00 ) then gam(1,j) = xplus(j)/(2.5D+00*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.5D+00*hksi(l1-1,j)/rhocon gam(l1,j) = rhocon/re if ( xplus(j) > 11.63D+00 ) then gam(l1,j) = xplus(j)/(2.5D+00*log(ewall*xplus(j))) end if end do end if do i = 2,l1-1 f(i,1,2) = 0.0D+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.0D+00 do j = 2,jcrys f(1,j,2) = 0.0D+00 f(l1,j,2) = 0.0D+00 end do f(l1,m1-1,2) = 0.0D+00 f(l1,m1,2) = 0.0D+00 ! ! 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.0D+00 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.0D+00 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.24D+00*((pr/sigt)**0.75-1.0D+00) & *(1.0D+00+0.28D+00*exp(-0.007D+00*pr/sigt)) if ( mode == 0) then do i = 2,l1-1 yplus(i) = re*rho(i,2)*sqrt(f(i,2,7))*cm4*0.5D+00*heta(i,2)/rhocon gam(i,1) = rpr if ( yplus(i) > 11.63D+00 ) then gam(i,1) = max(rhocon/re, & yplus(i)/(sigt*(2.5D+00*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.5D+00*heta(i,m1-1)/rhocon gam(i,m1) = rpr if ( yplus(i) > 11.63D+00 ) then gam(i,m1) = max(rhocon/re, & yplus(i)/(sigt*(2.5D+00*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.5D+00*hksi(2,j)/rhocon gam(1,j) = rpr if ( xplus(j) > 11.63D+00 ) 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.63D+00 ) then gam(l1,j) = max(rhocon/re, & xplus(j)/(sigt*(2.5D+00*log(ewall*xplus(j))+pbig))) end if end do end if do j = 2,jcrys f(l1,j,5) = 0.0D+00 end do do j = jcrys+1,m1 f(l1,j,5) = f(l1-1,j,5)-0.1D+00*birad*((f(l1,j,5)*(tw-tf)+tf)**4-tal**4) & *(x(l1,j)-x(l1-1,j))/100000000.0D+00 end do do j = 1,m1 f(1,j,5) = 0.5D+00 + 0.5D+00*yc(2,j) end do do i = 1,l1 f(i,m1,5) = 1.0D+00 end do do i = 1,l1 gam(i,1) = 0.0D+00 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.0D+00 ) then modelm = 0.0D+00 else if ( hamag > 0.0D+00 .and. hamag <= 20.0D+00 ) then modelm = 2 else if ( hamag > 20.0D+00 .and. hamag <= 40.0D+00 ) then modelm = 5 else if ( hamag > 40.0D+00 ) 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.0D+00/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.0D+00 end do do j = jcrys+1,m1 gam(l1,j) = 0.0D+00 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.0D+00*(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.0D+00*(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.0D+00*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.0D-05) 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.0D-05)) ap(i,j) = ap(i,j)-ce2*cmu*rho(i,j)**2*abs(f(i,j,7))/(gamt(i,j)+1.0D-05) end do end do end if return end subroutine gradnt(heta,hksi,i,j,phi,dphidx,dphidy,xc,yc) ! !*********************************************************************** ! !! GRADNT estimats the gradient of PHI at a primary node index (I,J). ! ! ! Discussion: ! ! The routine is given: ! ! PHI, the value of a state quantity at the primary nodes; ! HETA, HKSI, the two dimensions of the cell, as measured ! in the X, Y coordinate system; ! I, J, the indices of the primary node at which a gradient ! is desired; ! XC, YC, the X, Y coordinates of the corner nodes that define the ! shape of the cells, and the location of the midside nodes, and ! primary nodes. ! ! GRADNT must estimate ! ! DPHIDX, DPHIDY, an estimate for the gradient (d PHI/d X, d PHI/d Y) ! of the physical quantity PHI at the primary node (I,J). ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision HETA(NI,NJ). ! HETA is the length in the (X,Y) coordinate system, of the cell ! in the ETA direction. This is the actual length of the line ! connecting the mid-side nodes that separate primary node (I,J) ! from primary nodes (I,J-1) and (I,J+1). ! ! Input, double precision HKSI(NI,NJ). ! HKSI is the length in the (X,Y) coordinate system of the cell ! in the KSI direction. This is the actual length of the line ! connecting the mid-side nodes that separate primary node (I,J) ! from primary nodes (I-1,J) and (I+1,J). ! ! Input, integer I, J. ! I and J are the row and column indices of the ! primary node at which the gradient of PHI is desired. ! ! Input, double precision PHI(NI,NJ). ! PHI is the physical quantity whose value is given at each of ! the primary nodes. ! ! Output, double precision DPHIDX, DPHIDY. ! DPHIDX and DPHIDY are the computed approximations to the values ! of the partial derivatives d PHI/d X and d PHI/d Y. ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision detadx double precision detady double precision dpdeta double precision dpdksi double precision dphidx double precision dphidy double precision dxdeta double precision dxdksi double precision dksidx double precision dksidy double precision dydeta double precision dydksi double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer j double precision phi(ni,nj) double precision phie double precision phin double precision phis double precision phiw double precision t1 double precision t2 double precision volume double precision xc(ni,nj) double precision xe double precision xn double precision xs double precision xw double precision yc(ni,nj) double precision ye double precision yn double precision ys double precision yw ! ! We are given the value of some quantity PHI at the primary nodes. ! ! We want to estimate the gradients dPHI/dKSI and dPHI/dETA at ! those primary nodes. ! ! To do so, we must first estimate the value of PHI at the midpoint nodes, ! which lie between the given primary node and each of its four neighbors, ! above, below, to the right and to the left. ! ! Our task is made a little easier because of the fact that, in the ! KSI and ETA coordinate system, primary nodes are separated by ! exactly one unit, (and so are corner nodes, and midside nodes). ! ! This will make it easy to estimate PHI as a function of KSI and ! ETA, and hence to get dPHI/dKSI and dPHI/dETA by using a difference ! of two values. ! ! However, since we really want to compute dPHI/dX and dPHI/dY, ! we must make some adjustments to this problem: ! ! 1) the interpolated value of PHI at the midside node must not be ! the average of the two nearby primary node values, but the WEIGHTED ! average. The weights come from the actual X, Y distances between ! the primary nodes and the midside node in the "real world". ! ! 2) in the X, Y coordinate system, some nodes are separated by ! zero distance. This is just to help us handle boundary conditions, ! but we have to be careful to avoid getting a zero divisor in ! our weights. ! ! 3) we start out by calculating dPHI/dKSI and dPHI/dETA. To get ! dPHI/dX and dPHI/dY, we must compute the jacobian of the ! transformation that maps KSI and ETA to X and Y. Then we ! have to compute the INVERSE of this 2 by 2 matrix, to get ! the terms dKSI/dX, dKSI/dY, dETA/dX and dETA/dY. Then we ! can use the chain rule to calculate terms like: ! ! dPHI/dX = dPHI/dKSI * dKSI/dX + dPHI/dETA * dETA/dX ! ! ! Now let us look at our primary node in the simple (KSI,ETA) ! coordinate system, where: ! ! the (I,J) primary node has (KSI,ETA) coordinates (I,J). ! ! the midside nodes that neighbor primary node (I,J) have ! (KSI,ETA) coordinates as follows: ! ! North: (I, J+1/2) ! South: (I, J-1/2) ! East: (I+1/2, J) ! West: (I-1/2, J) ! ! the corner nodes that form the corners of the cell containing primary ! node (I,J) have (KSI,ETA) coordinates as follows: ! ! Northwest: (I-1/2, J+1/2) ! Northeast: (I+1/2, J+1/2) ! Southwest: (I-1/2, J-1/2) ! Southeast: (I+1/2, J-1/2) ! ! However, if we want to get the X, Y coordinates of a corner node, ! we can't use these (KSI,ETA) subscripts as indices, because FORTRAN ! only accepts whole number indices. So the appropriate indices to ! use to access the XC, YC arrays are: ! ! Northwest: (I, J+1) ! Northeast: (I+1, J+1) ! Southwest: (I, J) ! Southeast: (I+1, J) ! ! In the following diagram, we are interested primary node (I,J). ! ! We show the cell containing (I,J), the neighbor four primary nodes, ! the corner nodes, and the midside nodes. ! ! ! ^ | (I, J+1) | ! | | | ! | | | ! | ----------[I-1/2, J+1/2]--{I, J+1/2}--[I+1/2, J+1/2]---- ! | | | ! E | | ! T (I-1, J) {I-1/2, J] (I, J) {I+1/2, J} (I+1/2, J) ! A | | ! | | | ! | ----------[I-1/2, J-1/2]--{I, J-1/2}--[I+1/2, J-1/2]------ ! | | | ! | | | ! | | (I, J-1) | ! | ! | ! +-------------------KSI ( = I for primary nodes) -------------> ! ! ! Make a weighted average of PHI at the east and at the center, ! to estimate PHIE, the value of PHI at the east midside node: ! if ( hksi(i,j)+hksi(i+1,j) /= 0.0D+00 ) then t1 = hksi(i,j)/(hksi(i,j)+hksi(i+1,j)) t2 = hksi(i+1,j)/(hksi(i,j)+hksi(i+1,j)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phie = t1*phi(i+1,j)+t2*phi(i,j) ! ! Make a weighted average of PHI at the north and at the center, ! to estimate PHIN, the value of PHI at the north midside node. ! if ( heta(i,j)+heta(i,j+1) /= 0.0D+00 ) then t1 = heta(i,j)/(heta(i,j)+heta(i,j+1)) t2 = heta(i,j+1)/(heta(i,j)+heta(i,j+1)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phin = t1*phi(i,j+1)+t2*phi(i,j) ! ! Make a weighted average of PHI at the west and at the center, ! to estimate PHIW, the value of PHI at the west midside node. ! if ( hksi(i,j)+hksi(i-1,j) /= 0.0D+00 ) then t1 = hksi(i,j)/(hksi(i,j)+hksi(i-1,j)) t2 = hksi(i-1,j)/(hksi(i,j)+hksi(i-1,j)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phiw = t1*phi(i-1,j)+t2*phi(i,j) ! ! Make a weighted average of PHI at the south and at the center, ! to estimate PHIS, the value of PHI at the south midside node. ! if ( heta(i,j)+heta(i,j-1) /= 0.0D+00 ) then t1 = heta(i,j)/(heta(i,j)+heta(i,j-1)) t2 = heta(i,j-1)/(heta(i,j)+heta(i,j-1)) else t1 = 0.5D+00 t2 = 0.5D+00 end if phis = t1*phi(i,j-1)+t2*phi(i,j) ! ! Now subtract opposing values, to estimate dPHI/dKSI and dPHI/dETA ! at the primary node. We are assuming that DELTA KSI = DELTA ETA = 1. ! dpdksi = phie-phiw dpdeta = phin-phis ! ! Now we have to convert our results, which are in terms of KSI and ! ETA, into the X and Y coordinate system. ! ! Using the X, Y coordinates of the "corner nodes", compute ! the X, Y coordinates of the midside nodes of the control volume ! interfaces that surround the primary node with coordinates (I,J). ! ! ^ ! | [I-1/2, J+1/2]-----{I, J+1/2}----[I+1/2, J+1/2] ! E | | ! T {I-1/2, J} (I, J) {I+1/2, J} ! A | | ! | [I-1/2, J-1/2]-----{I, J-1/2}----[I+1/2, J] ! | ! +----------------------XSI--------------------------> ! xe = 0.5D+00*(xc(i+1,j+1)+xc(i+1,j)) ye = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)) xw = 0.5D+00*(xc(i,j+1)+xc(i,j)) yw = 0.5D+00*(yc(i,j+1)+yc(i,j)) xn = 0.5D+00*(xc(i,j+1)+xc(i+1,j+1)) yn = 0.5D+00*(yc(i,j+1)+yc(i+1,j+1)) xs = 0.5D+00*(xc(i,j)+xc(i+1,j)) ys = 0.5D+00*(yc(i,j)+yc(i+1,j)) ! ! From the X and Y coordinates of the midside nodes, ! estimate dX/dKSI, dX/dETA, dY/dKSI and dY/dETA at the primary node, ! again assuming a spacing delta KSI or delta ETA = 1 ! for opposing midside nodes. ! dxdksi = xe-xw dxdeta = xn-xs dydksi = ye-yw dydeta = yn-ys ! ! 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 ! ! If the determinant is zero, then the control volume has zero ! volume, which should never happen. ! if ( volume == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRADNT - Fatal error!' write ( *, '(a)' ) ' The control volume has zero volume!' stop end if ! ! Now compute the elements of the inverse Jacobian. ! ! J_Inverse(XSI,ETA) = ( dKSI/dX dKSI/dY) ! ( dETA/dX dETA/dY) ! ! using the simple fact that, for a 2 by 2 matrix, ! ! A_Inverse = (a22 -a21) / determinant(A) ! (-a12 a11) ! dksidx = dydeta/volume dksidy = -dxdeta/volume detadx = -dydksi/volume detady = dxdksi/volume ! ! Now we can finally compute dPHI/dX and dPHI/dY using the chain rule. ! dphidx = dpdksi*dksidx + dpdeta*detadx dphidy = dpdeta*detady + dpdksi*dksidy return end subroutine hello ( liwork, lwork ) ! !*********************************************************************** ! !! HELLO prints a brief message identifying the program. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer liwork integer lwork ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CRYSTAL_QED:' write ( *, '(a)' ) ' Version of 29 October 1996' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' CRYSTAL_QED seeks a set of parameters that minimize a ' write ( *, '(a)' ) ' cost functional associated with a crystallization ' write ( *, '(a)' ) ' problem.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' LIWORK, integer workspace size = ',liwork write ( *, '(a,i6)' ) ' LWORK, real workspace size = ',lwork return end subroutine inidat(ae1,ae2,ak1,ak2,area,b,b1jbl,b2jbl,b3jbl, & 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, & gam,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,npc,ns, & nsolve,ntimes,orth,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,vort,x,xbot,xc,xlen,y,ybot,yc,ylen) ! !*********************************************************************** ! !! INIDAT sets the initial values of certain data. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Output, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer maxbot integer ni integer nj integer nk integer ns ! double precision ae1(ni,nj) double precision ae2(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision area(ni,nj) double precision b double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision birad double precision bo double precision cappa double precision cd double precision ce1 double precision ce2 double precision cfo double precision cinc double precision cl double precision cmu double precision cost double precision cs double precision cvn double precision delt double precision dtm double precision epsad double precision epsil double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fkl double precision fks double precision fksl double precision fma double precision fn(ni,nj,ns) double precision fnu double precision fo(ni,nj,ns) double precision fr double precision frsl double precision gam(ni,nj) double precision gamt(ni,nj) double precision grash double precision grav double precision hamag double precision heta(ni,nj) double precision hf double precision 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) double precision orth double precision pr double precision ra double precision rdtm double precision re double precision recb double precision rect double precision relax(nk) double precision rho(ni,nj) double precision rhocon double precision rhol double precision rhos double precision rpr double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision sige double precision sigk double precision sigma double precision sigt double precision smooth double precision stel double precision stes double precision tal double precision tanca double precision tanca2 double precision tas double precision tend double precision tf double precision tin double precision tinit character ( len = 25 ) title(ns) double precision tnow double precision tw double precision vave double precision vol(ni,nj) double precision vort(ni,nj) double precision x(ni,nj) double precision xbot(maxbot) double precision xc(ni,nj) double precision xlen double precision y(ni,nj) double precision ybot(maxbot) double precision yc(ni,nj) double precision ylen ! ! Set quantities whose values, initial values, or dummy values ! do not depend on other quantities. ! do i = 1,ni do j = 1,nj ae1(i,j) = 0.0D+00 ae2(i,j) = 0.0D+00 ak1(i,j) = 0.0D+00 ak2(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj area(i,j) = 0.0D+00 end do end do b = 0.1778D+00 do i = 1,nj b1jbl(i) = 0.0D+00 b2jbl(i) = 0.0D+00 end do do i = 1,ni b3jbl(i) = 0.0D+00 end do cappa = 0.4187D+00 cd = 1.0D+00 ce1 = 1.44D+00 ce2 = 1.92D+00 cinc = 0.0D+00 cl = 1000.0D+00 cmu = 0.09D+00 cost = 0.0D+00 cs = 1000.0D+00 cvn = 200000.0D+00 dtm = 100.0D+00 epsad = 0.0001D+00 epsil = 0.0001D+00 ewall = 9.793D+00 do i = 1,ni do j = 1,nj do k = 1,ns f(i,j,k) = 0.0D+00 end do end do end do fkl = 64.0D+00 fks = 22.0D+00 fma = -1000.0D+00 do i = 1,ni do j = 1,nj do k = 1,ns fn(i,j,k) = 0.0D+00 end do end do end do fnu = 3.0D-07 do i = 1,ni do j = 1,nj do k = 1,ns fo(i,j,k) = 0.0D+00 end do end do end do fr = 1.0D-10 do i = 1,ni do j = 1,nj gam(i,j) = 0.0D+00 end do end do grash = 10000000.0D+00 grav = 9.81D+00 hamag = 0.0D+00 heta(1:ni,1:nj) = 0.0D+00 hf = 1800000.0D+00 do i = 1,ni do j = 1,nj hksi(i,j) = 0.0D+00 end do end do ! ! ICOST = 1, U**2+V**2 ! 2, (T-TF)**2 ! 3, (VAVE-SQRT(U**2+V**2)) ! icost = 1 icrys = 42 inturb = 0 ! ! IPLOT = 0, do not create a plot file. ! 1, create plot file. ! 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 do i = 1,4 lblk(i) = .false. end do do i = 5,10 lblk(i) = .true. end do lortho = .false. do i = 1,ns lsolve(i) = .false. end do 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.0D+00 pr = 0.015D+00 rdtm = 0.0D+00 re = 1.0D+00 recb = 0.0D+00 rect = 0.0D+00 relax(1) = 0.3D+00 relax(2) = 0.3D+00 relax(3) = 0.8D+00 relax(4) = 1.0D+00 relax(5) = 0.7D+00 relax(6) = 0.2D+00 relax(7) = 0.5D+00 relax(8) = 0.5D+00 relax(9) = 1.0D+00 relax(10) = 1.0D+00 relax(11) = 1.0D+00 relax(12) = 0.6D+00 relax(13) = 1.0D+00 relax(14) = 1.0D+00 rhocon = 1.0D+00 rhol = 2490.0D+00 rhos = 2490.0D+00 do i = 1,ni do j = 1,nj rueta(i,j) = 0.0D+00 ruksi(i,j) = 0.0D+00 end do end do sige = 1.3D+00 sigk = 1.0D+00 sigma = 0.72D+00 sigt = 0.9D+00 smooth = 0.0001D+00 tal = 1523.0D+00 tanca = 1.0D+00 tanca2 = 1.0D+00 tas = 1523.0D+00 tf = 1683.0D+00 tin = 1523.0D+00 tinit = 0.0D+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' tw = 1713.0D+00 vave = 1850.0D+00 do i = 1,ni do j = 1,nj vol(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj vort(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj x(i,j) = 0.0D+00 end do end do do i = 1,ni do j = 1,nj xc(i,j) = 0.0D+00 end do end do xlen = 1.2 do i = 1,ni do j = 1,nj y(i,j) = 0.0D+00 yc(i,j) = 0.0D+00 end do end do ylen = 1.0 ! ! Set things that depend on things. ! birad = 0.25D+00 * 5.67D+00 * 1.5D+00 * 0.0254D+00 /fks/(tw-tf)/re bo = (rhol*grav*b**2)/sigma cfo = fnu/(b*b) if ( inturb == 1 ) then do i = 1,ni do j = 1,nj f(i,j,7) = 0.005D+00 end do end do 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 do i = 1,ni do j = 1,nj f(i,j,8) = f(i,j,7)**1.5 / 0.006D+00 gamt(i,j) = 0.006D+00 *cmu*rho(i,j)*f(i,j,7)**0.5 end do end do end if ! ! Set the crucible shape. ! xbot(1) = 0.0D+00 ybot(1) = 0.0D+00 xbot(2) = 0.015D+00 ybot(2) = 0.3D+00 xbot(3) = 0.03D+00 ybot(3) = 0.5D+00 xbot(4) = 0.046D+00 ybot(4) = 0.6D+00 xbot(5) = 0.07D+00 ybot(5) = 0.7D+00 xbot(6) = 0.10D+00 ybot(6) = 0.8D+00 xbot(7) = 0.14D+00 ybot(7) = 0.9D+00 xbot(8) = 0.18D+00 ybot(8) = 0.95D+00 xbot(9) = 0.25D+00 ybot(9) = 1.0D+00 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. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Output, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision free integer i integer icrys integer j integer jcrys integer l0 integer licrys integer ll1 integer llst1 integer m0 double precision t1 double precision t2 double precision wave double precision xbot(nbot) double precision xc(ni,nj) double precision xlen double precision ybot(nbot) double precision yc(ni,nj) double precision ylen ! ll1 = icrys/2+1 llst1 = jcrys/2+1 licrys = (jcrys+m0)/2+1 free = 0.5D+00*xlen ! ! Set the Y coordinates ! do i = 2, l0 if ( i < ll1 ) then t2 = free*(0.5D+00*(dble(i-2)/dble(ll1-2))**1.5) else if ( i < icrys ) then t2 = free*(1.0D+00-0.5D+00*(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.4D+00 + 0.02D+00 * sin ( t2 * 20.0D+00 * 3.14159D+00 ) if ( j < llst1 ) then t1 = wave/2.0D+00*(dble(j-2)/dble(llst1-2))**1.5 else if ( j < jcrys ) then t1 = wave-wave/2.0D+00*(dble(jcrys-j)/dble(jcrys-llst1))**1.5 else if ( j < licrys ) then t1 = 0.4D+00 + 0.3D+00*(dble(j-jcrys)/dble(licrys-jcrys))**1.5 else t1 = 1.0D+00 - 0.3D+00*(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.5D+00 * (dble(i-2)/dble(ll1-2))**1.5 xc(i,j) = (1.0D+00 - t1)*xc(2,j)+t1*free else if ( i <= icrys ) then t1 = 1.0D+00 - 0.5D+00 * (dble(icrys-i)/dble(icrys-ll1))**1.5 xc(i,j) = ( 1.0D+00 - t1)*xc(2,j)+t1*free else t1 = (dble(i-icrys)/dble(l0-icrys))**1.2 if ( j < jcrys ) then xc(i,j) = xc(icrys,j)+t1*((xlen-free)+0.5D+00*yc(l0,jcrys)**2 & -0.5D+00 * 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 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. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input/output, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision ac double precision as double precision b1jbl(nj) double precision b2jbl(nj) double precision bo double precision coeff double precision coff double precision d2hdy2(nj) double precision delt double precision dhdy(nj) double precision dhdym(nj) double precision dhdyp(nj) double precision dlen double precision dx1 double precision dxm1 double precision dy1 double precision dym1 double precision fksl double precision flow(ni) double precision flomax double precision flomin double precision fr double precision frsl double precision high(ni) double precision hilev double precision himax double precision himin integer i integer icrys integer interl integer iprint integer j integer jcrys integer k integer l0 integer m0 double precision omega double precision omega2 double precision p(ni,nj) double precision pcsurf double precision pr double precision pres(ni) double precision pullv double precision re double precision rr(nj) double precision stel double precision tanca double precision tanca2 double precision xb(ni) double precision xc(ni,nj) double precision xlen double precision yc(ni,nj) double precision 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.0D+00 + (dhdy(j)*dhdy(j)))**1.5 & +dhdy(j)/yc(icrys,j)/(1.0D+00 + (dhdy(j)*dhdy(j)))**0.5 end do pres(jcrys+1) = 0.0D+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.0D+00 himin = 0.0D+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.0D+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.0 ) 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.0D+00 flomin = 0.0D+00 if ( iprint > 0 ) then if ( k == 5 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a,g14.6)' ) ' Himin = ',himin write ( *, '(a,g14.6)' ) ' Himax = ',himax write ( *, '(a,g14.6)' ) ' Omega = ',omega write ( *, '(a)' ) ' ' 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a)' ) ' Free surface movement' end if end if end do omega2 = 0.0D+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.003D+00/(flomax+1.0D-06)) omega2 = min(omega2, 0.003D+00/(-flomin+1.0D-06)) 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a,g14.6)' ) ' Flomin = ',flomin write ( *, '(a,g14.6)' ) ' Flomax = ',flomax write ( *, '(a,g14.6)' ) ' Omega = ',omega2 write ( *, '(a)' ) ' ' 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MOVGRD:' write ( *, '(a)' ) ' Solid-liquid interface movement' write ( *, '(a)' ) ' ' 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.5D+00*yc(l0,jcrys)**2 - 0.5D+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 some sample data from the ongoing solution. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! double precision cfo integer iter integer izone double precision res(ns) double precision smax double precision ssum double precision t(ni,nj) double precision tnow double precision u(ni,nj) double precision w(ni,nj) ! ! Solid zone output. ! if ( izone == 1 ) then if ( iter == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'OUTPUT:' write ( *, '(a,g14.6)' ) ' Current time = ',tnow write ( *, '(a,g14.6,a)' ) ' = ',tnow*cfo,' seconds.' write ( *, '(a)' ) ' Solid zone results.' write ( *, '(a)' ) ' Iter Res(5) T(50,21) T(50,22)' write ( *, '(a)' ) ' ' 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'OUTPUT:' write ( *, '(a)' ) ' Liquid zone results.' write ( *, '(a)' ) ' SMAX is the maximum local mass imbalance.' write ( *, '(a)' ) ' SSUM is the total mass imbalance.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Iter SMAX SSUM U(5,5) W(5,5) T(5,5)' write ( *, '(a)' ) ' ' 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 interpolates the pressure at boundary corners. ! ! ! 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. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! integer i integer ipref integer j integer jcrys integer jpref integer l1 integer m1 double precision p(ni,nj) double precision 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. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer ns ! double precision b double precision birad double precision bo double precision cappa double precision cfo double precision cvn double precision delt double precision dtm double precision fcsl double precision fksl double precision fma double precision fnu double precision fr double precision frsl double precision 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) double precision orth double precision pr double precision ra double precision rdtm double precision recb double precision rect double precision rhocon double precision smooth double precision stel double precision stes double precision tanca double precision tanca2 double precision tend double precision tf double precision tinit character ( len = 25 ) title(ns) double precision tw double precision vave double precision xlen double precision ylen ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PRDAT:' write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' B, Crucible radius = ',b write ( *, '(a,g14.6)' ) ' BIRAD, thermal coefficient = ',birad write ( *, '(a,g14.6)' ) ' BO, the Bond number = ',bo write ( *, '(a,g14.6)' ) ' CAPPA, wall function constant = ',cappa write ( *, '(a,g14.6)' ) ' CFO, time scale = ',cfo write ( *, '(a,g14.6)' ) ' CVN, ADAPT cell volume weight = ',cvn write ( *, '(a,g14.6)' ) ' DELT, time step in seconds = ',delt write ( *, '(a,g14.6)' ) ' DTM, dimensionless time step = ',dtm write ( *, '(a,g14.6)' ) ' FCSL = CS/CL = ',fcsl write ( *, '(a,g14.6)' ) ' FKSL = FKS/FKL = ',fksl write ( *, '(a,g14.6)' ) ' FMA, Marangoni number FMA = ',fma write ( *, '(a,g14.6)' ) ' FR, Froude number = ',fr write ( *, '(a,g14.6)' ) ' FRSL = RHOS/RHOL = ',frsl write ( *, '(a,g14.6)' ) ' GRASH, Grashof number = ',grash write ( *, '(a,i6)' ) ' ICOST, cost functional = ',icost write ( *, '(a,i6)' ) ' ICRYS, maximum I of crystal = ',icrys write ( *, '(a,i6)' ) ' INTURB, turbulence option = ',inturb write ( *, '(a,i6)' ) ' IPRINT, printing option = ',iprint write ( *, '(a,i6)' ) ' JCRYS, maximum J of crystal = ',jcrys write ( *, '(a,i6)' ) ' L0, number of I nodes = ',l0 write ( *, '(a,i6)' ) ' LAST, number of zone iterations on ' write ( *, '(a,i6)' ) ' each time step = ',last write ( *, '(a,i6)' ) ' LASTT, number of time steps = ',lastt write ( *, '(a,i6)' ) ' M0, number of J nodes = ',m0 write ( *, '(a,i6)' ) ' MODE, 0 cartesian, 1 axisymmetric = ',mode write ( *, '(a,i6)' ) ' NBOT, number of boundary points = ',nbot write ( *, '(a,g14.6)' ) ' ORTH, ADAPT orthogonality weight = ',orth write ( *, '(a,g14.6)' ) ' PR, Prandtl number PR = ',pr write ( *, '(a,g14.6)' ) ' RA, Rayleigh number = ',ra write ( *, '(a,g14.6)' ) ' RBD = FMA/RA = ',fma/ra write ( *, '(a,g14.6)' ) ' RDTM, = 1/DTM or 0 = ',rdtm write ( *, '(a,g14.6)' ) ' RECB, crucible Reynolds number = ',recb write ( *, '(a,g14.6)' ) ' RECT, crystal Reynolds number = ',rect write ( *, '(a,g14.6)' ) ' RHOCON, density constant = ',rhocon write ( *, '(a,g14.6)' ) ' SMOOTH, ADAPT smoothness weight = ',smooth write ( *, '(a,g14.6)' ) ' STEL, liquid Stefan number STEL = ',stel write ( *, '(a,g14.6)' ) ' STES, solid Stefan number STES = ',stes write ( *, '(a,g14.6)' ) ' TANCA = ',tanca write ( *, '(a,g14.6)' ) ' TANCA2 = ',tanca2 write ( *, '(a,g14.6)' ) ' TEND, dimensionless end time = ',tend write ( *, '(a,g14.6)' ) ' TF, crystal melting temperature = ',tf write ( *, '(a,g14.6)' ) ' TINIT, dimensionless start time = ',tinit write ( *, '(a,g14.6)' ) ' TW, the wall temperature = ',tw write ( *, '(a,g14.6)' ) ' VAVE, desired average velocity = ',vave write ( *, '(a,g14.6)' ) ' WEBER = FR*BO = ',fr*bo write ( *, '(a,g14.6)' ) ' XLEN = problem region length = ',xlen write ( *, '(a,g14.6)' ) ' YLEN = problem region height = ',ylen write ( *, '(a,g14.6)' ) ' Characteristic velocity FNU/B = ',fnu/b write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Variable Nonlinear Linear' write ( *, '(a)' ) ' Iterations Iterations' write ( *, '(a)' ) ' ' do i = 1,ns write ( *, '(a,i6,i6)' ) 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 residual of the linear equations. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ap(ni,nj) double precision con(ni,nj) double precision conmax double precision diamax double precision f(ni,nj,ns) integer i integer imax integer iprint integer j integer jmax integer l1 integer m1 integer nf double precision offmax double precision ratio double precision ratmin double precision res double precision resmax double precision sum double precision xmax ! iprint = 0 ! conmax = 0.0D+00 diamax = 0.0D+00 offmax = 0.0D+00 ratmin = 100000.0D+00 resmax = 0.0D+00 xmax = 0.0D+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.0 ) 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RESID' write ( *, '(a,i6)' ) ' Maxixum residual for variable ',nf write ( *, '(a,g14.6)' ) ' is ',resmax write ( *, '(a,2i6)' ) ' at I,J = ',imax,jmax i = imax j = jmax write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' A(I,J)*X(I,J) = ',ap(i,j),f(i,j,nf) write ( *, '(a,2g14.6)' ) ' A(I-1,J)*X(I-1,J) = ',aim(i,j),f(i-1,j,nf) write ( *, '(a,2g14.6)' ) ' A(I+1,J)*X(I+1,J) = ',aip(i,j),f(i+1,j,nf) write ( *, '(a,2g14.6)' ) ' A(I,J-1)*X(I,J-1) = ',ajm(i,j),f(i,j-1,nf) write ( *, '(a,2g14.6)' ) ' A(I,J+1)*X(I,J+1) = ',ajp(i,j),f(i,j+1,nf) write ( *, '(a,g14.6)' ) ' CON(I,J) = ',con(i,j) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Norm of variable is ',xmax write ( *, '(a,g14.6)' ) ' Norm of RHS is ',conmax write ( *, '(a,g14.6)' ) ' Norm of diagonal is ',diamax write ( *, '(a,g14.6)' ) ' Norm of off-diag is ',offmax write ( *, '(a,g14.6)' ) ' Min diag/off-diag ',ratmin end if return end subroutine rswrit(b1jbl,b2jbl,b3jbl,cost,e,gamt,icrys,jcrys, & l0,m0,nbot,p,pc,psi,rueta,ruksi,t,te,tk,tnow,u,v,vort,w,x, & xbot,xc,y,ybot,yc) ! !*********************************************************************** ! !! RSWRIT writes out information which can be used for restarts or plots. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real COST, the value of the cost functional ! associated with this solution. ! ! Input, real E(NI,NJ), the magnetic stream function. ! ! Input, real GAMT(NI,NJ), the diffusion coefficient. ! ! Input, integer ICRYS, JCRYS, the I and J coordinates of ! the lower left corner of the crystal. ! ! Input, integer L0, the number of rows of data. ! ! Input, integer M0, the number of columns of data. ! ! Input, integer NBOT, the number of crucible nodes. ! ! Input, real P(NI,NJ), the pressure. ! ! Input, real PC(NI,NJ), the corrected pressure. ! ! Input, real PSI(NI,NJ), the stream function. ! ! Input, real RUETA(NI,NJ), RUKSI(NI,NJ), the ! the momentum in the ETA and KSI directions. ! ! Input, real T(NI,NJ), the temperature. ! ! Input, real TE(NI,NJ), the turbulent epsilon. ! ! Input, real TK(NI,NJ), the turbulent K. ! ! Output, real TNOW, the current time. ! ! Input, real U(NI,NJ), the horizontal velocity. ! ! Input, real V(NI,NJ), the vertical velocity. ! ! Input, real VMAG(NI,NJ), the velocity magnitude. ! ! Input, real VORT(NI,NJ), the vorticity. ! ! Input, real W(NI,NJ), the axial velocity. ! ! Input, real X(NI,NJ), the X coordinates of the primary ! nodes. ! ! Input, real XBOT(NBOT), the X coordinates of the crucible ! bottom nodes. ! ! Input, real XC(NI,NJ), the X coordinates of the corner ! nodes. ! ! Input, real Y(NI,NJ), the Y coordinates of the primary ! nodes. ! ! Input, real YBOT(NBOT), the Y coordinates of the crucible ! bottom nodes. ! ! Input, real YC(NI,NJ), the Y coordinates of the corner ! nodes. ! implicit none ! integer nbot integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision cost logical, parameter :: debug = .false. double precision e(ni,nj) double precision gamt(ni,nj) integer i integer icrys integer j integer jcrys integer l0 integer m0 double precision p(ni,nj) double precision pc(ni,nj) double precision psi(ni,nj) double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision t(ni,nj) double precision te(ni,nj) double precision tk(ni,nj) double precision tnow double precision u(ni,nj) double precision v(ni,nj) double precision vort(ni,nj) double precision w(ni,nj) double precision x(ni,nj) double precision xbot(nbot) double precision xc(ni,nj) double precision y(ni,nj) double precision ybot(nbot) double precision yc(ni,nj) ! integer imax integer imin integer jmax integer jmin double precision vmax double precision vmin ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RSWRIT:' write ( *, '(a)' ) ' Writing restart information.' write ( *, '(a)' ) ' ' ! ! Temporary check that VORT contains legitimate data. ! if ( debug ) then vmax = vort(1,1) imax = 1 jmax = 1 vmin = vort(1,1) imin = 1 jmin = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I, J, VORT(I,J)' write ( *, '(a)' ) ' ' do i = 1,l0 do j = 1,m0 if ( vort(i,j) > vmax ) then vmax = vort(i,j) imax = i jmax = j end if if ( vort(i,j) < vmin ) then vmin = vort(i,j) imin = i jmin = j end if write(*,'(2i6,g14.5)')i,j,vort(i,j) end do end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RSWRIT - Temporary note:' write ( *, '(a,g14.6)' ) ' Minimum vorticity is ',vmin write ( *, '(a,i6,a,i6)' ) ' IMIN = ',imin,' JMIN=',jmin write ( *, '(a,g14.6)' ) ' Maximum vorticity is ',vmax write ( *, '(a,i6,a,i6)' ) ' IMAX = ',imax,' JMAX=',jmax write ( *, '(a,g14.6)' ) ' VMIN = ', vmin write ( *, '(a,g14.6)' ) ' VMAX = ', vmax end if ! ! Delete any old copy of the restart file. ! open(unit = 10,file='rswrit.txt',status='old',err=10) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RSWRIT - Note:' write ( *, '(a)' ) ' Deleting a previous copy of the restart file.' close(unit = 10,status='delete') 10 continue open(unit = 10,file='rswrit.txt',status='unknown') write(10,*)cost write(10,*)l0 write(10,*)jcrys write(10,*)icrys write(10,*)m0 write(10,*)nbot write(10,*)tnow ! ! NEW NEW NEW ! write(10,'(6g14.5)')(b1jbl(i),i = 1,nj) write(10,'(6g14.5)')(b2jbl(i),i = 1,nj) write(10,'(6g14.5)')(b3jbl(i),i = 1,ni) write(10,'(6g14.5)')((e(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((gamt(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((p(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((pc(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((psi(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((rueta(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((ruksi(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((t(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((te(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((tk(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((u(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((v(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((vort(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((w(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((x(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')(xbot(i),i = 1,nbot) write(10,'(6g14.5)')((xc(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')((y(i,j),i = 1,l0),j=1,m0) write(10,'(6g14.5)')(ybot(i),i = 1,nbot) write(10,'(6g14.5)')((yc(i,j),i = 1,l0),j=1,m0) close(unit = 10) return end subroutine setcst ( area, cinc, gam, icost, icrys, jcrys, l0, m0, ni, & nj, t, tf, u, v, vave, vort, xc, yc ) ! !*********************************************************************** ! !! SETCST computes a portion on the cost functional. ! ! ! Discussion: ! ! The cost functional is currently the integral over time and ! space of the norm of the fluid velocity. ! ! 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: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision AREA(NI,NJ). ! AREA(I,J) is the area of the control volume (I,J), ! centered on primary node (I,J), with corner nodes ! (I,J), (I+1,J), (I,J+1) and (I+1,J+1). ! ! Output, double precision CINC. ! CINC is the value of F(T) at the current time, that is, the ! integral of the cost function over the current region. ! ! Input, double precision GAM(NI,NJ). ! GAM(I,J) contains the transport coefficient for the current equation ! at primary node (I,J). GAM(I,J) is general the fluid viscosity ! (equations for variables 1, 2, 3, or 4), or the thermal diffusion ! coefficient (equations for variable 5). ! GAM is only needed for cost function ICOST = 5. If SETCST is ! called with ICOST = 5, then routine GAMSOR should have just ! been called, to set the value of GAM for the heat equation. ! ! Input, integer ICOST. ! ICOST chooses the cost functional to be minimized. The cost ! functional is an integral over time (from TINIT to TEND) and space ! (the region Omega). The integration in space is normally done ! first, yielding a function F(T) to be integrated over time. ! ! Depending on the value of ICOST, F(T) is: ! ! ICOST = 1: Velocity magnitude: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! (U(T,X,Y)**2 + V(T,X,Y)**2) dX dY ) ! ! ICOST = 2: Temperature deviation from TF: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! (T(T,X,Y)-TF)**2 dX dY ) ! ! ICOST = 3: Velocity magnitude deviation from Vave: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! ( Vave - SQRT(U(T,X,Y)**2 + V(T,X,Y)**2)) dX dY ) ! ! ICOST = 4: Flow velocity under the crystal: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T) and below CRYSTAL(T) ) ! U(T,X,Y)**2 + V(T,X,Y)**2 dX dY) ) ! ! For now, we take "below the crystal" to mean ALL control volumes ! under the crystal, for I = 1 to ICRYS-1. ! ! ICOST = 5: Heat flux from liquid into crystal and free surface: ! ! F(T) = Integral ( (X,Y) in LIQUID(T) and (CRYSTAL(T) or VOID(T)) ! GAM(X,Y) * dT(X,Y)/dNormal dX dY) ! ! ICOST = 6: Vorticity of the flow field: ! ! F(T) = SQRT( Integral ( (X,Y) in LIQUID(T)) ! VORT(T,X,Y)**2 dX dY ) ! ! Input, integer ICRYS. ! ICRYS specifies the end of the crystal in the I array direction, ! and in the vertical coordinate direction. ! ! Input, integer JCRYS. ! JCRYS specifies the end of the crystal in the J array direction, ! and in the horizontal coordinate direction. ! ! Input, integer L0. ! L0 is the extent of the grid in the vertical or "I" coordinate. ! ! Input, integer M0. ! M0 is the extent of the grid in the horizontal or "J" coordinate. ! ! Input, integer NI. ! The maximum number of grid points in the "I" or vertical direction. ! ! Input, integer NJ. ! The maximum number of grid points in the "J" or horizontal direction. ! ! Input, double precision T(NI,NJ). ! The current temperature at each primary node (I,J). ! ! Input, double precision TF. ! The fusion or melting temperature of the crystal. ! ! Input, double precision U(NI,NJ). ! The U component of velocity at each primary node (I,J). ! ! Input, double precision V(NI,NJ). ! The V component of velocity at each primary node (I,J). ! ! Input, double precision VAVE. ! For cost function ICOST = 3, the desired average velocity magnitude. ! ! Input, double precision VORT(NI,NJ). ! The fluid flow vorticity at each primary node (I,J). ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer ni integer nj ! double precision area(ni,nj) double precision cinc double precision dist double precision dtdn double precision gam(ni,nj) double precision gamma double precision hf integer i integer icost integer icrys integer j integer jcrys integer l0 integer m0 double precision t(ni,nj) double precision t1 double precision t2 double precision tf double precision u(ni,nj) double precision v(ni,nj) double precision vave double precision vort(ni,nj) double precision xc(ni,nj) double precision yc(ni,nj) ! ! #1: Velocity magnitude: ! if ( icost == 1 ) then cinc = 0.0D+00 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) ! ! #2: Temperature deviation from TF: ! else if ( icost == 2 ) then cinc = 0.0D+00 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) ! ! #3: Velocity magnitude deviation from VAVE: ! else if ( icost == 3 ) then cinc = 0.0D+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) ! ! #4: Flow velocity under the crystal: ! ! For now, we take "under the crystal" to mean ALL control volumes ! under the crystal, for I = 1 to ICRYS-1. ! else if ( icost == 4 ) then cinc = 0.0D+00 do i = 1,icrys-1 do j = 1,jcrys cinc = cinc+area(i,j)*((u(i,j)**2+v(i,j)**2)) end do end do cinc = sqrt(cinc) ! ! #5: Heat flux through crystal and free surface: ! else if ( icost == 5 ) then cinc = 0.0D+00 ! ! Add up the heat flux passing into the BOTTOM of the crystal. ! I believe this involves heat passing from Control Volume ! CV(ICRYS-1,J) to Control Volume CV(ICRYS,J), for J = 2 to JCRYS. ! do j = 2,jcrys t1 = 0.5D+00 t2 = 0.5D+00 gamma = (t1*gam(icrys,j+1)+t2*gam(icrys,j))/(t1+t2) dist = sqrt((xc(icrys,j+1)-xc(icrys,j))**2+(yc(icrys,j+1)-yc(icrys,j))**2) dtdn = (t(icrys,j+1)-t(icrys,j))/dist hf = gamma*dtdn cinc = cinc+hf*(xc(icrys,j+1)-xc(icrys,j)) end do ! ! Add up the heat flux passing into the void from the free surface. ! This involves heat passing from Control Volume CV(ICRYS-1,J) ! to Control Volume CV(ICRYS,J) for J = JCRYS+1 to M0. ! do j = jcrys+1,m0 t1 = 0.5D+00 t2 = 0.5D+00 gamma = (t1*gam(icrys,j+1)+t2*gam(icrys,j))/(t1+t2) dist = sqrt((xc(icrys,j+1)-xc(icrys,j))**2+(yc(icrys,j+1)-yc(icrys,j))**2) dtdn = (t(icrys,j+1)-t(icrys,j))/dist hf = gamma*dtdn cinc = cinc+hf*(xc(icrys,j+1)-xc(icrys,j)) end do ! ! #6: Vorticity of the flow field: ! else if ( icost == 6 ) then cinc = 0.0D+00 do i = 1,l0 do j = 1,m0 cinc = cinc+area(i,j)*vort(i,j)**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. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, double precision AE1(NI,NJ). ! d ALPHA/d ETA in thesis. ! I think AE1 is the value of d ETA/d X or d ETA/d Y. ! ! Output, double precision AE2(NI,NJ). ! d BETA/d ETA in thesis. ! I think AE2 is the value of d KSI/d X or d KSI/d Y. ! ! Output, double precision AK1(NI,NJ). ! d ALPHA/d KSI in thesis. ! I think AK1 is the value of d ETA/d X or d ETA/d Y. ! ! Output, double precision AK2(NI,NJ). ! d BETA/d KSI in thesis. ! I think AK2 is the value of d KSI/d X or d KSI/d Y. ! ! Output, double precision HETA(NI,NJ). ! HETA(I,J) is the physical (that is, "X,Y") length of the control ! volume in the ETA direction. ! ! HETA(I,J) is the length of the line connecting the mid-side ! nodes that separate primary node (I,J) from primary nodes ! (I,J-1) and (I,J+1). ! ! Warning: The code only defines HETA for the solid or liquid ! region it is working on, and sets it to zero elsewhere! ! ! HKSI Output, double precision HKSI(NI,NJ). ! HKSI(I,J) is the physical (that is, "X,Y") length of the control ! volume in the KSI direction. ! ! HKSI(I,J) is the length of the line connecting the mid-side ! nodes that separate primary node (I,J) from primary nodes ! (I-1,J) and (I+1,J). ! ! Warning: The code only defines HKSI for the solid or liquid ! region it is working on, and sets it to zero elsewhere! ! ! L1 Input, integer L1. ! L1 is some portion of the vertical grid size. ! L1 could equal L0, or ICRYS, depending on which subproblem is ! being solved. ! ! M1 Input, integer M1. ! M1 is some portion of the horizontal grid size. ! M1 could equal M0, or JCRYS, depending on which subproblem is ! being solved. ! ! MODE Input, integer MODE. ! MODE determines the kind of 2D geometry to be used. ! If MODE = 0, rectangular geometry is used; ! If MODE = 1, axisymmetric cylindrical geometry is used. ! ! NI Input, integer NI. ! NI is the maximum number of grid points in the I or row or vertical ! direction. ! ! NJ Input, integer NJ. ! NJ is the maximum number of grid points in the J or column or ! horizontal direction. ! ! R Output, double precision R(NI,NJ). ! R is used in axisymmetric problems to store the radial distance. ! For MODE = 0, R(I,J) = 1. ! For MODE = 1, R(I,J) = Y(I,J), the radial distance. ! ! VOL Output, double precision VOL(NI,NJ). ! VOL contains the volume of each control volume. ! Unfortunately, this is not completely accurate. ! ! The code only defines the volume for the control volumes that happen ! to be in the current subregion of interest. Outside of that ! region, the value of VOL cannot be relied on. The subregion of ! interest is specified by the values of L1 and M1. ! ! Also, for an axisymmetric problem (MODE = 1), VOL is automatically ! multiplied by a factor representing the radial distance. ! ! X Output, double precision X(NI,NJ). ! X contains the X coordinates of the primary nodes, which are ! the centers of the control volumes. ! ! In most cases, X(I,J) is the average of the four enclosing ! corner values, XC(I,J), XC(I+1,J), XC(I,J+1), and XC(I+1,J+1). ! ! XC Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Y Output, double precision Y(NI,NJ). ! Y contains the Y coordinates of primary nodes, which are ! the centers of the control volumes. ! ! In most cases, Y(I,J) is the average of the four enclosing ! corner values, YC(I,J), YC(I+1,J), YC(I,J+1), and YC(I+1,J+1). ! ! YC Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer ni integer nj ! double precision ae1(ni,nj) double precision ae2(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision dxdeta double precision dxdksi double precision dydeta double precision dydksi double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer j integer l1 integer m1 integer mode double precision r(ni,nj) double precision t1 double precision t2 double precision t3 double precision vol(ni,nj) double precision x(ni,nj) double precision xc(ni,nj) double precision xd double precision xjacb double precision xm double precision xp double precision xu double precision y(ni,nj) double precision yc(ni,nj) double precision yd double precision ym double precision yp double precision yu ! ! Calculate HKSI and HETA, essentially the width and height of ! the control volume, and VOL, the volume of the control volume. ! 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.5D+00*(xc(i+1,j+1)+xc(i+1,j)) yp = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)) xm = 0.5D+00*(xc(i,j+1)+xc(i,j)) ym = 0.5D+00*(yc(i,j+1)+yc(i,j)) xu = 0.5D+00*(xc(i,j+1)+xc(i+1,j+1)) yu = 0.5D+00*(yc(i,j+1)+yc(i+1,j+1)) xd = 0.5D+00*(xc(i,j)+xc(i+1,j)) yd = 0.5D+00*(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.0D+00 heta(i,m1) = 0.0D+00 vol(i,1) = 0.0D+00 vol(i,m1) = 0.0D+00 if ( i /= 1 .and. i.ne.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.0D+00 hksi(i,m1) = 0.0D+00 end if end do ! ! Take care of values along the borders I = 1 and I=L1. ! do j = 1,m1 hksi(1,j) = 0.0D+00 hksi(l1,j) = 0.0D+00 vol(1,j) = 0.0D+00 vol(l1,j) = 0.0D+00 if ( j /= 1 .and. j.ne.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.0D+00 heta(l1,j) = 0.0D+00 end if end do ! ! Calculate AK1 and AK2, whose meaning has not been revealed to me. ! ! ! ^ ! | C[I,J+1] ! | | ! E | ! T P(I-1,J) $ P(I,J) ! A | ! | | ! | C[I,J] ! | ! | ! +---------XSI---------> ! 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.0D+00 dydksi = dydksi*2.0D+00 end if t1 = sqrt(dxdeta**2+dydeta**2) t2 = sqrt(dxdksi**2+dydksi**2) t3 = dxdksi*dxdeta+dydeta*dydksi xjacb = dxdksi*dydeta-dxdeta*dydksi if ( xjacb /= 0.0D+00 ) then ak1(i,j) = t2*t1*t1/xjacb ak2(i,j) = t1*t3/xjacb else ak1(i,j) = 0.0D+00 ak2(i,j) = 0.0D+00 end if end do end do ! ! Calculate AE1 and AE2, whose meaning has not been revealed to me. ! ! ^ ! | P(I,J) ! | ! E ! T C[I,J]-----$-------C[I+1,J] ! A ! | ! | P(I,J-1) ! | ! +-------------XSI-------------> ! 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.0D+00 dydeta = dydeta*2.0D+00 end if t1 = sqrt(dxdeta**2+dydeta**2) t2 = sqrt(dxdksi**2+dydksi**2) t3 = dxdksi*dxdeta+dydeta*dydksi xjacb = dxdksi*dydeta-dxdeta*dydksi if ( xjacb /= 0.0D+00 ) then ae1(i,j) = t1*t2*t2/xjacb ae2(i,j) = t2*t3/xjacb else ae1(i,j) = 0.0D+00 ae2(i,j) = 0.0D+00 end if end do end do ! ! Set the cylindrical geometry weight factor. ! if ( mode == 0 ) then do i = 1,l1 do j = 1,m1 r(i,j) = 1.0D+00 end do end do 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 setp ( fjeta, fjksi, heta, hksi, l1, m1, ni, nj, p ) ! !*********************************************************************** ! !! SETP estimates the values of pressure at the interfaces. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer ni integer nj ! double precision fjeta(ni,nj) double precision fjksi(ni,nj) double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer j integer l1 integer m1 double precision p(ni,nj) double precision tt ! do i = 2,l1 do j = 1,m1 if ( hksi(i,j)+hksi(i-1,j) /= 0.0D+00 ) then tt = hksi(i,j)/(hksi(i,j)+hksi(i-1,j)) else tt = 0.0D+00 end if fjksi(i,j) = tt*p(i-1,j)+(1.0D+00-tt)*p(i,j) end do end do do i = 1,l1 do j = 2,m1 if ( heta(i,j)+heta(i,j-1) /= 0.0D+00 ) then tt = heta(i,j)/(heta(i,j)+heta(i,j-1)) else tt = 0.0D+00 end if fjeta(i,j) = tt*p(i,j-1)+(1.0D+00-tt)*p(i,j) end do end do return end subroutine setru ( heta, hksi, l1, m1, ni, nj, rho, rueta, ruksi, u, & v, x, y ) ! !*********************************************************************** ! !! SETRU estimates the normal mass velocity component at interfaces. ! ! ! Discussion: ! ! SETRU estimates the normal mass velocity component at the east-west ! and north-south interfaces of the control volumes by interpolating ! values of velocity and density computed for the primary nodes. ! ! These values are needed only for cells which lie within the melt ! region. ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Double precision HETA(NI,NJ). ! HETA(I,J) is the physical (that is, "X,Y") length of the control ! volume in the ETA direction. ! ! HETA(I,J) is the length of the line connecting the mid-side ! nodes that separate primary node (I,J) from primary nodes ! (I,J-1) and (I,J+1). ! ! Warning: The code only defines HETA for the solid or liquid ! region it is working on, and sets it to zero elsewhere! ! ! Double precision HKSI(NI,NJ). ! HKSI(I,J) is the physical (that is, "X,Y") length of the control ! volume in the XI direction. ! ! HKSI(I,J) is the length of the line connecting the mid-side ! nodes that separate primary node (I,J) from primary nodes ! (I-1,J) and (I+1,J). ! Warning: The code only defines HKSI for the solid or liquid ! region it is working on, and sets it to zero elsewhere! ! ! Integer L1. ! L1 is some portion of the vertical grid size. ! L1 could equal L0, or ICRYS, depending on which subproblem is ! being solved. ! ! Integer M1. ! M1 is some portion of the horizontal grid size. ! M1 could equal M0, or JCRYS, depending on which subproblem is ! being solved. ! ! Integer NI. ! NI is the maximum number of grid points in the I or row or vertical ! or XI direction. ! ! Integer NJ. ! NJ is the maximum number of grid points in the J or column or ! horizontal or ETA direction. ! ! Double precision RHO(NI,NJ). ! RHO contains the fluid density at each primary node. ! In version 1.0 of MASTRAPP, RHO is has a constant uniform ! value, namely RHOCON. ! ! Double precision RUETA(NI,NJ). ! RUETA is the component of the mass velocity interpolated ! at the interface between primary nodes (I,J-1) and (I,J), ! and normal to the cell interface. ! ! Double precision RUKSI(NI,NJ). ! RUKSI is the component of the mass velocity interpolated ! at the interface between primary nodes (I-1,J) and (I,J), ! and normal to the cell interface. ! ! Double precision U(NI,NJ). ! U contains the horizontal (in the XY coordinate system) component of ! velocity at each primary node (I,J). ! ! Double precision V(NI,NJ). ! V contains the vertical (in the XY coordinate system) component of ! velocity at each primary node (I,J). ! ! Double precision X(NI,NJ). ! X contains the X coordinates of the primary nodes, which are ! the centers of the control volumes. ! In most cases, X(I,J) is the average of the four enclosing ! corner values, XC(I,J), XC(I+1,J), XC(I,J+1), and XC(I+1,J+1). ! ! Double precision Y(NI,NJ). ! Y contains the Y coordinates of primary nodes, which are ! the centers of the control volumes. ! In most cases, Y(I,J) is the average of the four enclosing ! corner values, YC(I,J), YC(I+1,J), YC(I,J+1), and YC(I+1,J+1). ! implicit none ! integer ni integer nj ! double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer j integer l1 integer m1 double precision rho(ni,nj) double precision rhov double precision rueta(ni,nj) double precision ruksi(ni,nj) double precision tt double precision u(ni,nj) double precision ucs double precision v(ni,nj) double precision vcs double precision vlen double precision x(ni,nj) double precision xcomp double precision y(ni,nj) double precision ycomp ! ! Compute RUKSI, the normal velocity at east-west cell interfaces. ! do i = 2,l1 do j = 1,m1 ! ! To interpolate a value at the cell interface, we interpolate various ! quantities using the relative distance of the (I-1,J) and (I,J) ! cells from their common interface. ! if ( hksi(i,j)+hksi(i-1,j) /= 0.0D+00 ) then tt = hksi(i,j)/(hksi(i,j)+hksi(i-1,j)) else tt = 0.0D+00 end if ! ! Interpolate the values of the velocity and density from the ! centers of the two cells to the intermediate cell interface. ! ucs = tt*u(i-1,j) +(1.0D+00-tt)*u(i,j) vcs = tt*v(i-1,j) +(1.0D+00-tt)*v(i,j) rhov = tt*rho(i-1,j)+(1.0D+00-tt)*rho(i,j) ! ! Compute the unit vector which is normal to the cell interface. ! xcomp = x(i,j)-x(i-1,j) ycomp = y(i,j)-y(i-1,j) vlen = sqrt(xcomp**2+ycomp**2) if ( vlen > 0.0D+00 ) then xcomp = xcomp/vlen ycomp = ycomp/vlen end if ! ! Interpolate the normal mass velocity to the cell interface. ! ruksi(i,j) = rhov*(ucs*xcomp+vcs*ycomp) end do end do ! ! Compute RUETA, the normal velocity at north-south cell interfaces. ! do i = 1,l1 do j = 2,m1 ! ! To interpolate a value at the cell interface, we interpolate various ! quantities using the relative distance of the (I,J-1) and (I,J) ! cells from their common interface. ! if ( heta(i,j)+heta(i,j-1) /= 0.0D+00 ) then tt = heta(i,j)/(heta(i,j)+heta(i,j-1)) else tt = 0.0D+00 end if ! ! Interpolate the values of the velocity and density from the ! centers of the two cells to the intermediate cell interface. ! ucs = tt*u(i,j-1)+(1.0D+00-tt)*u(i,j) vcs = tt*v(i,j-1)+(1.0D+00-tt)*v(i,j) rhov = tt*rho(i,j-1)+(1.0D+00-tt)*rho(i,j) ! ! Compute the unit vector which is normal to the cell interface. ! xcomp = x(i,j)-x(i,j-1) ycomp = y(i,j)-y(i,j-1) vlen = sqrt(xcomp**2+ycomp**2) if ( vlen > 0.0D+00 ) then xcomp = xcomp/vlen ycomp = ycomp/vlen end if ! ! Interpolate the normal mass velocity to the cell interface. ! rueta(i,j) = rhov*(ucs*xcomp+vcs*ycomp) end do end do return end subroutine setup(ae1,ae2,ak1,ak2,b1jbl,b2jbl,b3jbl,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. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision XC(NI,NJ). ! The X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! The Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nk = 14 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! double precision a1 double precision a2 double precision acof double precision ae1(ni,nj) double precision ae2(ni,nj) double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ak1(ni,nj) double precision ak2(ni,nj) double precision ap(ni,nj) double precision apr double precision b1jbl(nj) double precision b2jbl(nj) double precision b3jbl(ni) double precision birad double precision biu(nj) double precision biv(nj) double precision bju(ni) double precision bjv(ni) double precision blu(nj) double precision blv(nj) double precision bmu(ni) double precision bmv(ni) double precision bpi double precision bpi1 double precision bpj double precision bpj1 double precision bpm double precision cappa double precision cd double precision ce1 double precision ce2 double precision cmu double precision cofu(ni,nj,5) double precision cofv(ni,nj,5) double precision con(ni,nj) double precision conu(ni,nj) double precision conv(ni,nj) double precision denom double precision diff double precision dpeta double precision dpksi double precision dxdeta double precision dxdksi double precision dydeta double precision dydksi double precision em double precision ep double precision epsil double precision ewall double precision f(ni,nj,ns) double precision fcsl double precision fjeta(ni,nj) double precision fjksi(ni,nj) double precision fksl double precision flow double precision fma double precision fmax(ns) double precision fn(ni,nj,ns) double precision fo(ni,ni,ns) double precision frc double precision frsl double precision gam(ni,nj) double precision gamt(ni,nj) double precision grash double precision hamag double precision heta(ni,nj) double precision hetap(ni,nj) double precision hksi(ni,nj) double precision hksip(ni,nj) double precision hutop(ni,nj) double precision 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) double precision peta(ni,nj) double precision pksi(ni,nj) double precision pr double precision qt(nmaxij) double precision r(ni,nj) double precision rdtm double precision re double precision recb double precision rect double precision relax(nk) double precision res(ns) double precision rho(ni,nj) double precision rhocon double precision rpr double precision rueij double precision rueij1 double precision ruet double precision rueta(ni,nj) double precision ruki1j double precision rukij double precision ruksi(ni,nj) double precision rukt double precision sige double precision sigk double precision sigt double precision smax double precision soor double precision ssum double precision stel double precision stes double precision t1 double precision t2 double precision t3 double precision t4 double precision tal double precision tas double precision tem1 double precision tem2 double precision temp double precision tf double precision tmp(ni,nj) double precision tw double precision vol(ni,nj) double precision x(ni,nj) double precision xc(ni,nj) double precision xd double precision xl double precision xr double precision xu double precision y(ni,nj) double precision yc(ni,nj) double precision yd double precision yl double precision yr double precision 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.0D+00*gam(2,j)/hksi(2,j) flow = -ruksi(2,j) call diflow(acof,diff,flow) aim(2,j) = acof*ak1(2,j) diff = 2.0D+00*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.5D+00*hksi(i,j)*gam(i+1,j)+0.5D+00*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.0D+00*gam(i,2)/heta(i,2) flow = -rueta(i,2) call diflow(acof,diff,flow) ajm(i,2) = acof*ae1(i,2) diff = 2.0D+00*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.5D+00*heta(i,j)*gam(i,j+1)+0.5D+00*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, & b3jbl,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.0D+00*gam(2,j)/hksi(2,j) flow = -ruksi(2,j) call diflow(acof,diff,flow) aim(2,j) = acof*ak1(2,j) diff = 2.0D+00*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.5D+00*hksi(i,j)*gam(i+1,j) & +0.5D+00*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.0D+00*gam(i,2)/heta(i,2) flow = -rueta(i,2) call diflow(acof,diff,flow) ajm(i,2) = acof*ae1(i,2) diff = 2.0D+00*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.5D+00*heta(i,j)*gam(i,j+1) & +0.5D+00*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 biv(j) = gam(1,j) blv(j) = gam(l1,j) end do 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, & b3jbl,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.5D+00*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)-yc(i,j+1)-yc(i,j)) hksip(i,j) = 0.5D+00*(dxdksi*hutop(i,j)+dydksi*hvtop(i,j)) dxdeta = 0.5D+00*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5D+00*(yc(i+1,j+1)+yc(i,j+1)-yc(i+1,j)-yc(i,j)) hetap(i,j) = 0.5D+00*(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.0D+00 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.5D+00*(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.5D+00*(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.0D+00 do i = 2,l1-2 a1 = 0.5D+00*(ak1(i,j)+ak1(i+1,j)) a2 = 0.5D+00*(ak1(i+1,j)+ak1(i+2,j)) denom = 0.5D+00*(ap(i,j)*hksi(i,j)/a1+ap(i+1,j)*hksi(i+1,j)/a2) apr = 2.0D+00*denom/(hksi(i,j)+hksi(i+1,j)) dxdksi = 0.5D+00*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5D+00*(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.5D+00*(xc(i+2,j+1)+xc(i+2,j)-xc(i+1,j+1)-xc(i+1,j)) dydksi = 0.5D+00*(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.5D+00*(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.5D+00*(frc*(ruksi(i+1,j)*ak1(i+1,j)-rukt*ak1(i,j)) & +(1.0D+00-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.0D+00/denom aim(i+1,j) = aip(i,j) end do aip(l1-1,j) = 0.0D+00 aim(2,j) = 0.0D+00 end do do i = 2,l1-1 do j = 2,m1-2 a1 = 0.5D+00*(ae1(i,j)+ae1(i,j+1)) a2 = 0.5D+00*(ae1(i,j+1)+ae1(i,j+2)) denom = 0.5D+00*(ap(i,j)*heta(i,j)/a1+ap(i,j+1)*heta(i,j+1)/a2) apr = 2.0D+00*denom/(heta(i,j)+heta(i,j+1)) dxdeta = 0.5D+00*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5D+00*(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.5D+00*(xc(i+1,j+2)+xc(i,j+2)-xc(i+1,j+1)-xc(i,j+1)) dydeta = 0.5D+00*(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.5D+00*(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.5D+00*(frc*(rueta(i,j+1)*ae1(i,j+1)-ruet*ae1(i,j)) & +(1.0D+00-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.0D+00/denom ajm(i,j+1) = ajp(i,j) end do ajp(i,m1-1) = 0.0D+00 ajm(i,2) = 0.0D+00 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.0D+00 ssum = 0.0D+00 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.0D+04,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.0D+00-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.5D+00*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.0D+00*f(l1-1,j,nf)-pksi(l1-1,j) pksi(2,j) = 2.0D+00*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.5D+00*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.0D+00 *f(i,m1-1,nf)-peta(i,m1-1) peta(i,2) = 2.0D+00 *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.0D+00 ) then ajm(i,2) = 0.0D+00 end if if ( gam(i,m1) == 0.0D+00 ) then ajp(i,m1-1) = 0.0D+00 end if end do do j = 2,m1-1 if ( gam(1,j) == 0.0D+00 ) then aim(2,j) = 0.0D+00 end if if ( gam(l1,j) == 0.0D+00 ) then aip(l1-1,j) = 0.0D+00 end if end do do i = 2,l1-1 do j = 2,m1-1 xr = 0.5D+00*(xc(i+1,j+1)+xc(i+1,j)) yr = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)) xl = 0.5D+00*(xc(i,j+1)+xc(i,j)) yl = 0.5D+00*(yc(i,j+1)+yc(i,j)) xu = 0.5D+00*(xc(i+1,j+1)+xc(i,j+1)) yu = 0.5D+00*(yc(i+1,j+1)+yc(i,j+1)) xd = 0.5D+00*(xc(i+1,j)+xc(i,j)) yd = 0.5D+00*(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, & b3jbl,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.0D+00 ) then ajm(i,2) = 0.0D+00 end if if ( gam(i,m1) == 0.0D+00 ) then ajp(i,m1-1) = 0.0D+00 end if end do do j = 2,m1-1 if ( gam(1,j) == 0.0D+00 ) then aim(2,j) = 0.0D+00 end if if ( gam(l1,j) == 0.0D+00 ) then aip(l1-1,j) = 0.0D+00 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, & b3jbl,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, & b3jbl,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.0D+00 ap(i,j) = 0.0D+00 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, & b3jbl,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.5D+00*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)-yc(i,j+1)-yc(i,j)) hksip(i,j) = 0.5D+00*(dxdksi*hutop(i,j)+dydksi*hvtop(i,j)) dxdeta = 0.5D+00*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5D+00*(yc(i+1,j+1)+yc(i,j+1)-yc(i+1,j)-yc(i,j)) hetap(i,j) = 0.5D+00*(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.0D+00 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.5D+00*(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.5D+00*(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.0D+00 do i = 2,l1-2 a1 = 0.5D+00*(ak1(i,j)+ak1(i+1,j)) a2 = 0.5D+00*(ak1(i+1,j)+ak1(i+2,j)) denom = 0.5D+00*(ap(i,j)*hksi(i,j)/a1+ap(i+1,j)*hksi(i+1,j)/a2) apr = 2.0D+00*denom/(hksi(i,j)+hksi(i+1,j)) dxdksi = 0.5D+00*(xc(i+1,j+1)+xc(i+1,j)-xc(i,j+1)-xc(i,j)) dydksi = 0.5D+00*(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.5D+00*(xc(i+2,j+1)+xc(i+2,j)-xc(i+1,j+1)-xc(i+1,j)) dydksi = 0.5D+00*(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.5D+00*(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.5D+00*(frc*(ruksi(i+1,j)*ak1(i+1,j)-rukt*ak1(i,j)) & +(1.0D+00-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.0D+00/denom aim(i+1,j) = aip(i,j) end do aip(l1-1,j) = 0.0D+00 aim(2,j) = 0.0D+00 end do do i = 2,l1-1 do j = 2,m1-2 a1 = 0.5D+00*(ae1(i,j)+ae1(i,j+1)) a2 = 0.5D+00*(ae1(i,j+1)+ae1(i,j+2)) denom = 0.5D+00*(ap(i,j)*heta(i,j)/a1+ap(i,j+1)*heta(i,j+1)/a2) apr = 2.0D+00*denom/(heta(i,j)+heta(i,j+1)) dxdeta = 0.5D+00*(xc(i+1,j+1)+xc(i,j+1)-xc(i+1,j)-xc(i,j)) dydeta = 0.5D+00*(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.5D+00*(xc(i+1,j+2)+xc(i,j+2)-xc(i+1,j+1)-xc(i,j+1)) dydeta = 0.5D+00*(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.5D+00*(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.5D+00*(frc*(rueta(i,j+1)*ae1(i,j+1)-ruet*ae1(i,j)) & +(1.0D+00-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.0D+00/denom ajm(i,j+1) = ajp(i,j) end do ajp(i,m1-1) = 0.0D+00 ajm(i,2) = 0.0D+00 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.0D+00*f(l1-1,j,nf)-pksi(l1-1,j) pksi(2,j) = 2.0D+00*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.0D+00*f(i,m1-1,nf)-peta(i,m1-1) peta(i,2) = 2.0D+00*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.5D+00*(xc(i+1,j+1)+xc(i+1,j)) yr = 0.5D+00*(yc(i+1,j+1)+yc(i+1,j)) xl = 0.5D+00*(xc(i,j+1)+xc(i,j)) yl = 0.5D+00*(yc(i,j+1)+yc(i,j)) xu = 0.5D+00*(xc(i+1,j+1)+xc(i,j+1)) yu = 0.5D+00*(yc(i+1,j+1)+yc(i,j+1)) xd = 0.5D+00*(xc(i+1,j)+xc(i,j)) yd = 0.5D+00*(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 ! ! Get the transport coefficients, GAM. ! 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.0D+00*gam(1,j)/hksi(2,j) flow = -ruksi(2,j) call diflow(acof,diff,flow) aim(2,j) = acof*ak1(2,j) diff = 2.0D+00*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.5D+00*hksi(i,j)*gam(i+1,j)+0.5D+00*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.0D+00*gam(i,1)/heta(i,2) flow = -rueta(i,2) call diflow(acof,diff,flow) ajm(i,2) = acof*ae1(i,2) diff = 2.0D+00*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.5D+00*heta(i,j)*gam(i,j+1)+0.5D+00*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, & b3jbl,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 calculates the locations of the primary 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). ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer L. ! L specifies the portion of the vertical grid size to use. ! ! Input, integer M. ! M specifies the portion of the horizontal grid size to use. ! ! Input, integer NI. ! NI is the maximum number of grid points in the "I" or vertical ! direction. ! ! Input, integer NJ. ! NJ is the maximum number of grid points in the "J" or horizontal ! direction. ! ! Output, double precision X(NI,NJ). ! X coordinate of primary nodes, which are the centers of ! control volumes. ! In most cases, X(I,J) is the average of the four enclosing ! corner values, XC(I,J), XC(I+1,J), XC(I,J+1), and XC(I+1,J+1). ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Output, double precision Y(NI,NJ). ! Y coordinate of primary nodes, which are the centers of ! control volumes. ! In most cases, Y(I,J) is the average of the four enclosing ! corner values, YC(I,J), YC(I+1,J), YC(I,J+1), and YC(I+1,J+1). ! ! Input, double precision YC(NI,NJ). ! YC contains the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer ni integer nj ! integer i integer j integer l integer m double precision x(ni,nj) double precision xc(ni,nj) double precision y(ni,nj) double precision yc(ni,nj) ! ! The central set of values. ! do i = 2,l-1 do j = 2,m-1 x(i,j) = 0.25D+00*(xc(i,j)+xc(i+1,j)+xc(i,j+1)+xc(i+1,j+1)) y(i,j) = 0.25D+00*(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.5D+00*(xc(i,2)+xc(i+1,2)) y(i,1) = 0.5D+00*(yc(i,2)+yc(i+1,2)) x(i,m) = 0.5D+00*(xc(i,m)+xc(i+1,m)) y(i,m) = 0.5D+00*(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.5D+00*(xc(2,j)+xc(2,j+1)) y(1,j) = 0.5D+00*(yc(2,j)+yc(2,j+1)) x(l,j) = 0.5D+00*(xc(l,j)+xc(l,j+1)) y(l,j) = 0.5D+00*(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 certain entries of the linear system. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: ns = 10 ! double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ap(ni,nj) double precision cappa double precision cd double precision cm4 double precision cmu double precision con(ni,nj) double precision f(ni,nj,ns) double precision fo(ni,nj,ns) double precision heta(ni,nj) double precision 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.0D+00 aip(i,j) = 0.0D+00 aim(i,j) = 0.0D+00 ajp(i,j) = 0.0D+00 ajm(i,j) = 0.0D+00 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.0D+00 aip(i,j) = 0.0D+00 aim(i,j) = 0.0D+00 ajp(i,j) = 0.0D+00 ajm(i,j) = 0.0D+00 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.0D+00 aip(icrys,j) = 0.0D+00 aim(icrys,j) = 0.0D+00 ajp(icrys,j) = 0.0D+00 ajm(icrys,j) = 0.0D+00 con(icrys,j) = 0.0D+00 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.5D+00*heta(i,2)) ap(i,2) = 1.0D+00 con(i,m1-1) = cd*f(i,m1-1,7)**1.5/(cappa*cm4*0.5D+00*heta(i,m1-1)) ap(i,m1-1) = 1.0D+00 end do do j = 2,m1-1 con(2,j) = cd*f(2,j,7)**1.5/(cappa*cm4*0.5D+00*hksi(2,j)) ap(2,j) = 1.0D+00 con(l1-1,j) = cd*f(l1-1,j,7)**1.5/(cappa*cm4*0.5D+00*hksi(l1-1,j)) ap(l1-1,j) = 1.0D+00 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). ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 integer, parameter :: nmaxij = 64 integer, parameter :: ns = 10 ! double precision aim(ni,nj) double precision aip(ni,nj) double precision ajm(ni,nj) double precision ajp(ni,nj) double precision ap(ni,nj) double precision bl double precision blc double precision blm double precision blp double precision con(ni,nj) double precision denom double precision f(ni,nj,ns) integer i integer j integer l1 logical lblk(ns) integer m1 integer nf integer nsolve(ns) integer nt double precision pt(nmaxij) double precision qt(nmaxij) ! do nt = 1,nsolve(nf) if ( lblk(nf) ) then pt(1) = 0.0D+00 qt(1) = 0.0D+00 do i = 2,l1-1 bl = 0.0D+00 blp = 0.0D+00 blm = 0.0D+00 blc = 0.0D+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.0D-10 ) then pt(i) = 0.0D+00 qt(i) = 0.0D+00 else pt(i) = blp/denom qt(i) = (blc+blm*qt(i-1))/denom end if end do bl = 0.0D+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.0D+00 qt(1) = 0.0D+00 do j = 2,m1-1 bl = 0.0D+00 blp = 0.0D+00 blm = 0.0D+00 blc = 0.0D+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.0D-10 ) then pt(j) = 0.0D+00 qt(j) = 0.0D+00 else pt(j) = blp/denom qt(j) = (blc+blm*qt(j-1))/denom end if end do bl = 0.0D+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.0D+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.0D+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.0D+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.0D+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 subroutine vortic ( heta, hksi, icrys, l0, m0, u, v, vort, xc, yc ) ! !*********************************************************************** ! !! VORTIC computes the value of the vorticity. ! ! ! Discussion: ! ! VORTIC estimates the value of the vorticity function VORT(I,J) at ! each of the primary nodes (I,J), based on the values of the ! horizontal and vertical velocities U and V, the XSI and ETA spacings ! of the control volumes, and the X and Y coordinates of the ! corners of the control volumes. ! ! The fluid flow vorticity is defined as ! ! VORT(X,Y) = dV(X,Y)/dX - dU(X,Y)/dY ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision HETA(NI,NJ). ! HETA is the physical length of the control volume in the ETA direction. ! This is the actual length of the line connecting the mid-side ! nodes that separate primary node (I,J) from primary nodes ! (I,J-1) and (I,J+1). ! ! Input, double precision HKSI(NI,NJ). ! HKSI is the physical length of the control volume in the KSI direction. ! This is the actual length of the line connecting the mid-side ! nodes that separate primary node (I,J) from primary nodes ! (I-1,J) and (I+1,J). ! ! Input, integer ICRYS. ! ICRYS specifies the end of the crystal in the I array direction, ! and in the vertical coordinate direction. ! ! Input, integer L0. ! L0 is the extent of the grid in the vertical or "I" coordinate. ! ! Input, integer M0. ! M0 is the extent of the grid in the horizontal or "J" coordinate. ! ! Input, double precision U(NI,NJ). ! U is the horizontal component of velocity at each primary node (I,J). ! ! Input, double precision V(NI,NJ). ! V is the vertical component of velocity at each primary node (I,J). ! ! Output, double precision VORT(NI,NJ). ! VORT is the fluid flow vorticity at each primary node (I,J). ! ! Input, double precision XC(NI,NJ). ! XC contains the X coordinate of nodes which are the "corners" of ! control volumes. ! ! Input, double precision YC(NI,NJ). ! YC is the Y coordinate of nodes which are the "corners" of ! control volumes. ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! double precision dudx double precision dudy double precision dvdx double precision dvdy double precision heta(ni,nj) double precision hksi(ni,nj) integer i integer icrys integer j integer l0 integer m0 double precision u(ni,nj) double precision v(ni,nj) double precision vort(ni,nj) double precision xc(ni,nj) double precision yc(ni,nj) ! ! At all the interior nodes, get the vorticitiy. ! In the liquid region, get the spatial derivatives of the flow, ! but in the void and crystal regions, just set it to zero. ! do i = 2, l0-1 do j = 2, m0-1 if ( i < icrys ) then call gradnt(heta,hksi,i,j,u,dudx,dudy,xc,yc) call gradnt(heta,hksi,i,j,v,dvdx,dvdy,xc,yc) vort(i,j) = dvdx-dudy else vort(i,j) = 0.0D+00 end if end do end do ! ! Copy the vorticity data to the sides and corners. ! do i = 2,l0-1 vort(i,1) = vort(i,2) vort(i,m0) = vort(i,m0-1) end do do j = 2,m0 vort(1,j) = vort(2,j) vort(l0,j) = vort(l0-1,j) end do vort(1,1) = vort(2,2) vort(1,m0) = vort(2,m0-1) vort(l0,1) = vort(l0-1,2) vort(l0,m0) = vort(l0-1,m0-1) return end subroutine wrtec ( gamt, l1, m1, p, psi, t, te, tk, u, v, w, x, y ) ! !*********************************************************************** ! !! WRTEC writes solution data to a TECPLOT input file. ! ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real GAMT(NI,NJ). ! GAMT(I,J) is the diffusion coefficient at primary node (I,J). ! ! Input L1. ! L1 is the maximum value for I when indexing the data arrays. ! ! Input M1. ! M1 is the maximum value for J when indexing the data arrays. ! ! Input, real P(NI,NJ). ! P(I,J) is the pressure at primary node (I,J). ! ! Input, real PSI(NI,NJ). ! PSI(I,J) is the stream function at primary node (I,J). ! ! Input, real T(NI,NJ). ! T(I,J) is the temperature at primary node (I,J). ! ! Input, real TE(NI,NJ). ! TE(I,J) is the turbulent epsilon coefficient at ! primary node (I,J). ! ! Input, real TK(NI,NJ). ! TK(I,J) is the turbulent K coefficient at primary node (I,J). ! ! Input, real U(NI,NJ). ! U(I,J) is the horizontal velocity at primary node (I,J). ! ! Input, real V(NI,NJ). ! V(I,J) is the vertical velocity at primary node (I,J). ! ! Input, real W(NI,NJ). ! W(I,J) is the axial velocity at primary node (I,J). ! ! Input, real X(NI,NJ). ! X(I,J) is the X coordinate of primary node (I,J). ! ! Input, real Y(NI,NJ). ! Y(I,J) is the Y coordinate of primary node (I,J). ! implicit none ! integer, parameter :: ni = 64 integer, parameter :: nj = 64 ! character ( len = 11 ) filtec double precision gamt(ni,nj) integer i integer iunit integer j integer l1 integer m1 double precision p(ni,nj) double precision psi(ni,nj) double precision t(ni,nj) double precision te(ni,nj) double precision tk(ni,nj) double precision u(ni,nj) double precision v(ni,nj) double precision w(ni,nj) double precision x(ni,nj) double precision y(ni,nj) ! ! Some local variables. ! filtec = 'tecplot.txt' iunit = 10 ! ! Delete any old copy of the TECPLOT data file. ! open(unit = iunit,file=filtec,status='old',err=10) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'WRTEC - Note:' write ( *, '(a)' ) ' Deleting a previous copy of '//filtec close(unit = iunit,status='delete') 10 continue ! ! Open a new copy of the TECPLOT data file. ! open ( unit = iunit, file = filtec, status = 'new' ) ! ! Write the data. ! Note that the crystal program has X and Y reversed. ! To keep TECPLOT from going mad, we reverse X and Y, and U and V. ! write ( iunit, '(a)' )' title = ','"zhanghui"' write ( iunit, '(a)' )' variables = "x","y","u","v","w",'// & '"Press", "Temp","Psi","TK","TE","Gamt"' write ( iunit, '(a)' )' zone t = "zone 1" , i=',l1,' , j=',m1,' , f= block' write ( iunit, '(5g15.6)' ) ((y(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((x(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((v(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((u(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((w(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((p(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((t(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((psi(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((tk(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((te(i,j),i = 1,l1),j=1,m1) write ( iunit, '(5g15.6)' ) ((gamt(i,j),i = 1,l1),j=1,m1) ! ! Close the file. ! close ( unit = iunit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'WRTEC - Note:' write ( *, '(a)' ) ' TECPLOT data written to ' // trim ( filtec ) return end