SUBROUTINE IJMHD(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" 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" dimension bound_mix(lip1) #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in ijmhd.F::IJMHD(...)" #endif C C NOW THE ARRAYS FOR THE IHYDRO SWEEP C #ifdef PSCTEST * * test code myunit = mype+40 k=25 level3=25 if (level3 .eq. 25) then write (myunit,*) '>>>>>>>>>>>>>>>>> level3 =',level3,k do j=llow,nl_j+no2 do i=1,nl_i+1 write (myunit,978) i,j,k,tran_i(i,j,k,1), $ tran_i(i,j,k,2), $ tran_i(i,j,k,3) call flush(myunit) 978 format(' s ',3i6,3(1pe13.3)) enddo * call sleepit(1) enddo endif * * #endif nl_ip1 = nl_i+1 nl_jp1 = nl_j+1 nl_kp1 = nl_k+1 isweep = 1 CALL MIJSET(nl_i,nl_j) C DO 499 k=1,nl_k C level3 = k C DO 200 J=llow,nl_j+no2 DO 200 i=llow,NL_I+no2 VOLd(i,j) = VOLUME(I,J,K) BXd(i,j) = bx(i,j,k) BYd(i,j) = by(i,j,k) BZd(i,j) = bz(i,j,k) BX1d(i,j) = bx1(i,j,k) BY1d(i,j) = by1(i,j,k) BZ1d(i,j) = bz1(i,j,k) RHOd(i,j) = RHO(I,J,K) VXd(i,j) = VX(I,J,K) VYd(i,j) = VY(I,J,K) VZd(i,j) = VZ(I,J,K) Cd(i,j) = C(I,J,K) RHO1d(i,j) = RHO1(I,J,K) VX1d(i,j) = VX1(I,J,K) VY1d(i,j) = VY1(I,J,K) VZ1d(i,j) = VZ1(I,J,K) C1d(i,j) = C1(I,J,K) 200 continue do j=1,nl_j do i=1,nl_i RHO2d(i,j) = RHO2(I,J,K) VX2d(i,j) = VX2(I,J,K) VY2d(i,j) = VY2(I,J,K) VZ2d(i,j) = VZ2(I,J,K) C2d(i,j) = C2(I,J,K) BDVX2d(i,j) = BDVX2(I,J,K) BDVY2d(i,j) = BDVY2(I,J,K) BDVZ2d(i,j) = BDVZ2(I,J,K) enddo enddo C DO 220 J=1,NL_J DO 220 I=1,NL_IP1 FACEd(i,j) = FACE_i(I,J,K) do ltran = 1,6 transd(i,j,ltran) = tran_i(i,j,k,ltran) enddo dtx = transd(i,j,2)*transd(i,j,6)-transd(i,j,3)*transd(i,j,5) dty = transd(i,j,3)*transd(i,j,4)-transd(i,j,1)*transd(i,j,6) dtz = transd(i,j,1)*transd(i,j,5)-transd(i,j,2)*transd(i,j,4) dinv = 1./sqrt(dtx**2+dty**2+dtz**2) transd(i,j,7) = dtx*dinv transd(i,j,8) = dty*dinv transd(i,j,9) = dtz*dinv bnormd(i,j) = bi1(i,j,k) 220 CONTINUE C do 228 j=1,nl_j do 228 i=1,nl_ip1 bnqd(i,j) = bnqfi(i,j,k) bsqqd(i,j) = bsqqi(i,j,k) bxqd(i,j) = bxqfi(i,j,k) byqd(i,j) = byqfi(i,j,k) bzqd(i,j) = bzqfi(i,j,k) bxqnd(i,j) = bxqni(i,j,k) byqnd(i,j) = byqni(i,j,k) bzqnd(i,j) = bzqni(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 J=1,NL_JP1 DO 300 I=-no2+1,NL_I+no2 BJd(i,j) = Bj(I,J,K) BJ1d(i,j) = Bj1(I,J,K) 300 CONTINUE DO j=1,nl_j+1 do i =1,nl_i bj2d(i,j) = bj2(i,j,k) enddo enddo C DO 310 J=1,NL_J DO 310 I=1,NL_IP1 BI2d(i,j) = Bi2(I,J,K) 310 CONTINUE do j=-no2+1,nl_j+no2 do i=1,nl_i+1 BId(i,j) = Bi(I,J,K) BI1d(i,j) = Bi1(I,J,K) enddo enddo DO 320 J=1,NL_JP1 DO 320 I=1,NL_I 320 FACEjd(i,j) = FACE_J(I,J,k) C DO 332 J=1,NL_JP1 DO 330 I=1,NL_IP1 330 EDGEd(i,j) = EDGE_k(I,J,k) 332 CONTINUE C DO 340 J=0,NL_J+2 DO 340 I=0,NL_I+2 Xqd(i,j) = 0.5*(X(I,J,K)+X(I,J,K+1)) Yqd(i,j) = 0.5*(Y(I,J,K)+Y(I,J,K+1)) Zqd(i,j) = 0.5*(Z(I,J,K)+Z(I,J,K+1)) 340 CONTINUE C DO 352 J=1,NL_JP1 DO 352 I=1,NL_IP1 tran2d(i,j,1) = vec_k(i,j,k,1) tran2d(i,j,2) = vec_k(i,j,k,2) tran2d(i,j,3) = vec_k(i,j,k,3) tran2d(i,j,7) = vec_k(i,j,k,7) tran2d(i,j,8) = vec_k(i,j,k,8) tran2d(i,j,9) = vec_k(i,j,k,9) tran2d(i,j,4) = vec_k(i,j,k,8)*vec_k(i,j,k,3) - $ vec_k(i,j,k,9)*vec_k(i,j,k,2) tran2d(i,j,5) = vec_k(i,j,k,9)*vec_k(i,j,k,1) - $ vec_k(i,j,k,7)*vec_k(i,j,k,3) tran2d(i,j,6) = vec_k(i,j,k,7)*vec_k(i,j,k,2) - $ vec_k(i,j,k,8)*vec_k(i,j,k,1) dinv = 1./sqrt(tran2d(i,j,4)**2+tran2d(i,j,5)**2+ $ tran2d(i,j,6)**2 ) tran2d(i,j,4) = dinv*tran2d(i,j,4) tran2d(i,j,5) = dinv*tran2d(i,j,5) tran2d(i,j,6) = dinv*tran2d(i,j,6) e2(i,j) = ek_g(i,j,k) bound_mix(j) = bmix_ij(j,k) 352 CONTINUE C do j=llow,nl_j+no2 do i=1,nl_ip1 xnormdi(i,j) =tran_i(i,j,k,1) ynormdi(i,j) =tran_i(i,j,k,2) znormdi(i,j) =tran_i(i,j,k,3) facedi(i,j) = face_i(i,j,k) enddo enddo C C do j=1,nl_jp1 do i=llow,nl_i+no2 xnormdj(i,j) =tran_j(i,j,k,1) ynormdj(i,j) =tran_j(i,j,k,2) znormdj(i,j) =tran_j(i,j,k,3) facedj(i,j) = face_j(i,j,k) enddo enddo C do 356 j=1,nl_jp1 do 356 i=1,nl_ip1 bx0d(i,j) = bxqk(i,j,k) by0d(i,j) = byqk(i,j,k) bz0d(i,j) = bzqk(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 #ifdef PSCTEST * * test code myunit = mype+40 kq =25 write (myunit,*) '>>>>>>>>>>>>>>>>> level3 =',level3,kq do j=llow,nl_j+no2 do i=1,nl_i+1 write (myunit,979) i,j,kq,tran_i(i,j,kq,1), $ tran_i(i,j,kq,2), $ tran_i(i,j,kq,3) call flush(myunit) 979 format(' r ',3i6,3(1pe13.3)) enddo * call sleepit(1) enddo * * myunit = mype +40 if ( level3 .ge. 24 ) then write (myunit,*) '>>>>>>>>>>>>>>>>> level3 =',level3,k do j=llow,nl_j+no2 do i=1,nl_ip1 write (myunit,977) i,j,k,tran_i(i,j,k,1),tran_i(i,j,k,2), $ tran_i(i,j,k,3),xnormdi(i,j),ynormdi(i,j), $ znormdi(i,j) call flush(myunit) 977 format(3i6,6(1pe13.3)) enddo enddo endif #endif 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 J=1,NL_J DO 450 I=1,NL_I RHO2(I,J,K) = RHO2d(i,j) VX2(I,J,K) = VX2d(i,j) VY2(I,J,K) = VY2d(i,j) VZ2(I,J,K) = VZ2d(i,j) C2(I,J,K) = C2d(i,j) BDVX2(I,J,K) = BDVX2d(i,j) BDVY2(I,J,K) = BDVY2d(i,j) 450 BDVZ2(I,J,K) = BDVZ2d(i,j) C DO 470 J=1,NL_JP1 DO 470 I=1,NL_I 470 Bj2(I,J,K) = BJ2d(i,j) C DO 480 J=1,NL_J DO 480 I=1,NL_IP1 480 BI2(I,J,K) = BI2d(i,j) C if ( ilow_bound) then lower_i = 2 else lower_i =1 endif do j=1,nl_j+1 do i=lower_i,nl_i+1 ek_g(i,j,k) = e2(i,j) enddo enddo C 499 CONTINUE C C C NOW BACK TO MHD C C RETURN END