“…If (time.eq.1.or.mod(time,int((ttotal/60/delt))).eq.0) Then Do k=z2+1,z3 qz(k-z2)=(3.0*T(imax/2,y2+1,k)-4.0* T(imax/2,y2+2,k) & +T(imax/2,y2+3,k))*kcon(imax/2,y2+1,k)/2.0/dely/10000.0 End Do write(4,40)time*delt,qz(1), qz(2), qz(3), qz(4), qz (5), & qz (6), qz (7), qz (8), qz(9), qz(10), & qz (11), qz (12), qz (13), qz (14),qz (15), & qz (16), qz (17), qz (18), qz (19),qz (20), & qz (21), qz (22), qz (23), qz(24),qz (25), & qz (26), qz (27), qz (28), qz (29),qz (30), & qz (31), qz (32), qz (33), qz (34),qz (35), & qz (36), qz (37), qz (38), qz 39 Do k=5,kmax-4,2 rhoustar(imax-1,j,k)=rhoustar(imax-3,j,k) End Do End Do Do j=3,jmax-2,2 Do k=3,kmax-2,2 rhovstar(imax,j,k)=rhovstar(imax-2,j,k) End Do End Do Do j=4,jmax-3,2 Do k=4,kmax-3,2 rhowstar(imax-1,j,k)=rhowstar(imax-3,j,k) End Do End Do !using rhoustar, rhovstar, rhowstar, solve for the corrected pressure ! (denoted as pp)using the pressure correction formula which is solved !using the successive overrelaxation subroutine (sor) !prepare e matrix for sor subroutine !e can be thought of as a mass source term and represents how far the !continuity equation deviates from the converged answer.…”