SUBROUTINE JKMHD(DT,RHO,VX,VY,VZ,C,BX,BY,BZ,BI,BJ,BK, $ RHO1,VX1,VY1,VZ1,C1,BX1,BY1,BZ1,BI1,BJ1,BK1) #include "param.inc" #include "global_dims.inc" #include "mhd_P++_pointers.inc" #include "fortran_P++_arrays.inc" #include "var2.inc" #include "help.inc" #include "mdims.inc" #include "bzero.inc" #include "boundx.inc" dimension bound_mix(lip1) C DIMENSION RHO(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ VX(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ VY(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ VZ(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) DIMENSION C(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BX(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BY(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BZ(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) DIMENSION BI(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) $ ,BJ(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BK(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) DIMENSION RHO1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ VX1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ VY1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ VZ1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) DIMENSION C1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BX1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BY1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BZ1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) DIMENSION BI1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) $ ,BJ1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH), $ BK1(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH) #include "twod.inc" C C NOW THE ARRAYS FOR THE JHYDRO SWEEP C #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in jkmhd.F::JKMHD(...)" #endif nl_ip1 = nl_i+1 nl_jp1 = nl_j+1 nl_kp1 = nl_k+1 isweep = 2 CALL MIJSET(nl_j,nl_k) C DO 499 I=1,nl_i C level3 = i C DO 200 K=llow,nl_k+no2 DO 200 J=llow,NL_J+no2 VOLd(J,K) = VOLUME(I,J,K) BXd(J,K) = bx(i,j,k) BYd(J,K) = by(i,j,k) BZd(J,K) = bz(i,j,k) BX1d(J,K) = bx1(i,j,k) BY1d(J,K) = by1(i,j,k) BZ1d(J,K) = bz1(i,j,k) RHOd(J,K) = RHO(I,J,K) VXd(J,K) = VX(I,J,K) VYd(J,K) = VY(I,J,K) VZd(J,K) = VZ(I,J,K) Cd(J,K) = C(I,J,K) RHO1d(J,K) = RHO1(I,J,K) VX1d(J,K) = VX1(I,J,K) VY1d(J,K) = VY1(I,J,K) VZ1d(J,K) = VZ1(I,J,K) C1d(J,K) = C1(I,J,K) 200 continue do k=1,nl_k do j=1,nl_j RHO2d(J,K) = RHO2(I,J,K) VX2d(J,K) = VX2(I,J,K) VY2d(J,K) = VY2(I,J,K) VZ2d(J,K) = VZ2(I,J,K) C2d(J,K) = C2(I,J,K) BDVX2d(J,K) = BDVX2(I,J,K) BDVY2d(J,K) = BDVY2(I,J,K) BDVZ2d(J,K) = BDVZ2(I,J,K) enddo enddo C DO 220 K=1,NL_K DO 220 J=1,NL_JP1 FACEd(J,K) = FACE_j(I,J,K) do ltran = 1,6 transd(j,k,ltran) = tran_j(i,j,k,ltran) enddo dtx = transd(j,k,2)*transd(j,k,6)-transd(j,k,3)*transd(j,k,5) dty = transd(j,k,3)*transd(j,k,4)-transd(j,k,1)*transd(j,k,6) dtz = transd(j,k,1)*transd(j,k,5)-transd(j,k,2)*transd(j,k,4) dinv = 1./sqrt(dtx**2+dty**2+dtz**2) transd(j,k,7) = dtx*dinv transd(j,k,8) = dty*dinv transd(j,k,9) = dtz*dinv bnormd(j,k) = bj1(i,j,k) 220 CONTINUE C do 228 k=1,nl_k do 228 j=1,nl_jp1 bnqd(j,k) = bnqfj(i,j,k) bsqqd(j,k) = bsqqj(i,j,k) bxqd(j,k) = bxqfj(i,j,k) byqd(j,k) = byqfj(i,j,k) bzqd(j,k) = bzqfj(i,j,k) bxqnd(j,k) = bxqnj(i,j,k) byqnd(j,k) = byqnj(i,j,k) bzqnd(j,k) = bzqnj(i,j,k) 228 continue C CALL HYDRO(DT,BNORMD,VOLD,FACED,TRANSD, $ BNqD,bsqqd,BXqD,BYqD,BZqD,bxqnd,byqnd,bzqnd, $ RHOD,VXD,VYD,VZD,CD,BXD,BYD,BZD, $ RHO1D,VX1D,VY1D,VZ1D,C1D,BX1D,BY1D,BZ1D, $ RHO2D,VX2D,VY2D,VZ2D,C2D,BDVX2D,BDVY2D,BDVZ2D) C C C now the bfield push, I electric fields C C DO 300 K=1,NL_KP1 DO 300 J=-no2+1,nl_j+no2 BJd(J,K) = BK(I,J,K) BJ1d(J,K) = BK1(I,J,K) 300 CONTINUE do k=1,nl_k+1 do j=1,nl_j BJ2d(J,K) = BK2(I,J,K) enddo enddo C DO 310 K=1,NL_K DO 310 J=1,NL_JP1 BI2d(J,K) = BJ2(I,J,K) 310 CONTINUE do k=-no2+1,nl_k+no2 do j=1,nl_jp1 BId(J,K) = Bj(I,J,K) BI1d(J,K) = BJ1(I,J,K) enddo enddo c DO 320 K=1,NL_KP1 DO 320 J=1,NL_J 320 FACEJd(J,K) = FACE_K(I,J,1) C DO 332 K=1,NL_KP1 DO 330 J=1,NL_JP1 330 EDGEd(J,K) = EDGE_I(I,J,K) 332 CONTINUE C DO 340 K=0,NL_K+2 DO 340 J=0,NL_J+2 Xqd(J,K) = 0.5*(X(I,J,K)+X(I+1,J,K)) Yqd(J,K) = 0.5*(Y(I,J,K)+Y(I+1,J,K)) Zqd(J,K) = 0.5*(Z(I,J,K)+Z(I+1,J,K)) 340 CONTINUE C DO 352 K=1,NL_KP1 DO 352 J=1,NL_JP1 tran2d(j,k,1) = vec_i(i,j,k,1) tran2d(j,k,2) = vec_i(i,j,k,2) tran2d(j,k,3) = vec_i(i,j,k,3) tran2d(j,k,7) = vec_i(i,j,k,7) tran2d(j,k,8) = vec_i(i,j,k,8) tran2d(j,k,9) = vec_i(i,j,k,9) tran2d(j,k,4) = vec_i(i,j,k,8)*vec_i(i,j,k,3) - $ vec_i(i,j,k,9)*vec_i(i,j,k,2) tran2d(j,k,5) = vec_i(i,j,k,9)*vec_i(i,j,k,1) - $ vec_i(i,j,k,7)*vec_i(i,j,k,3) tran2d(j,k,6) = vec_i(i,j,k,7)*vec_i(i,j,k,2) - $ vec_i(i,j,k,8)*vec_i(i,j,k,1) dinv = 1./sqrt(tran2d(j,k,4)**2+tran2d(j,k,5)**2+ $ tran2d(j,k,6)**2 ) tran2d(j,k,4) = dinv*tran2d(j,k,4) tran2d(j,k,5) = dinv*tran2d(j,k,5) tran2d(j,k,6) = dinv*tran2d(j,k,6) e2(j,k) = ei_g(i,j,k) bound_mix(k) = 0. ! no-op for jk sweep 352 CONTINUE C do k=llow,nl_k+no2 do j=1,nl_jp1 xnormdi(j,k) =tran_j(i,j,k,1) ynormdi(j,k) =tran_j(i,j,k,2) znormdi(j,k) =tran_j(i,j,k,3) facedi(j,k) = face_j(i,j,k) enddo enddo C C do k=1,nl_kp1 do j=llow,nl_j+no2 xnormdj(j,k) =tran_k(i,j,k,1) ynormdj(j,k) =tran_k(i,j,k,2) znormdj(j,k) =tran_k(i,j,k,3) facedj(j,k) = face_k(i,j,k) enddo enddo C C do 356 k=1,nl_kp1 do 356 j=1,nl_jp1 bx0d(j,k) = bxqi(i,j,k) by0d(j,k) = byqi(i,j,k) bz0d(j,k) = bzqi(i,j,k) 356 continue C C C * WRITE (9,*) '************* I=', I,' **************************' * CALL WROTE('t(1)',TRAN2D(1,1,1),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(2)',TRAN2D(1,1,2),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(3)',TRAN2D(1,1,3),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(4)',TRAN2D(1,1,4),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(5)',TRAN2D(1,1,5),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(6)',TRAN2D(1,1,6),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(7)',TRAN2D(1,1,7),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(8)',TRAN2D(1,1,8),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) * CALL WROTE('t(9)',TRAN2D(1,1,9),1,nl_jp1,nl_jp1,1,nl_kp1,nl_kp1) C * write(9,*) '**************************************************' * write(9,*) 'JKMHD, i= ',i * write(9,*) '**************************************************' * write(9,*) '**************************************************' * call print2d(vxd,'vxd',mi,mj,li,lj) * call print2d(vyd,'vyd',mi,mj,li,lj) * call print2d(vzd,'vzd',mi,mj,li,lj) * call print2d(bi1d,'bi1d',mip1,mj,lip1,lj) * call print2d(bj1d,'bj1d',mi,mjp1,li,ljp1) *c CALL BFIELD(DT,EDGEd, $ TRAN2D,facedi,facedj, $ Xqd,Yqd,Zqd, $ RHOD,VXD,VYD,VZD, $ BID,BJD, $ RHO1D,VX1D,VY1D,VZ1D, $ BI1D,BJ1D, $ BI2D,BJ2d,e2,bound_mix, $ xnormdj,ynormdj,znormdj, $ xnormdi,ynormdi,znormdi, $ bx0d,by0d,bz0d) C DO 450 K=1,NL_K DO 450 J=1,NL_J RHO2(I,J,K) = RHO2d(J,K) VX2(I,J,K) = VX2d(J,K) VY2(I,J,K) = VY2d(J,K) VZ2(I,J,K) = VZ2d(J,K) C2(I,J,K) = C2d(J,K) BDVX2(I,J,K) = BDVX2d(J,K) BDVY2(I,J,K) = BDVY2d(J,K) 450 BDVZ2(I,J,K) = BDVZ2d(J,K) C DO 470 K=1,NL_KP1 DO 470 J=1,NL_J 470 BK2(I,J,K) = BJ2d(J,K) C DO 480 K=1,NL_K DO 480 J=1,NL_JP1 480 BJ2(I,J,K) = BI2d(J,K) C do k=1,nl_k+1 do j=1,nl_j+1 ei_g(i,j,k) = e2(j,k) enddo enddo C 499 CONTINUE C C C NOW BACK TO MHD C C RETURN END