C SUBROUTINE KIMHD(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 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) C C DIMENSIONS AND EQUIVALENCES FOR THE (K,I) SWEEP C #include "twod.inc" dimension bound_mix(lip1) #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in kimhd.F::KIMHD(...)" #endif C C nl_ip1 = nl_i+1 nl_jp1 = nl_j+1 nl_kp1 = nl_k+1 C isweep = 3 CALL MIJSET(Nl_k,nl_i) C C NOW FOR THE KHYDRO SWEEP C DO 999 J=1,Nl_j C level3 = j C C DO 200 I=llow,nl_i+no2 DO 200 K=llow,NL_K+no2 VOLd(K,I) = VOLUME(I,J,K) BXd(K,I) = BX(I,J,K) BYd(K,I) = BY(I,J,K) BZd(K,I) = BZ(I,J,K) BX1d(K,I) = BX1(I,J,K) BY1d(K,I) = BY1(I,J,K) BZ1d(K,I) = BZ1(I,J,K) RHOd(K,I) = RHO(I,J,K) VXd(K,I) = VX(I,J,K) VYd(K,I) = VY(I,J,K) VZd(K,I) = VZ(I,J,K) Cd(K,I) = C(I,J,K) RHO1d(K,I) = RHO1(I,J,K) VX1d(K,I) = VX1(I,J,K) VY1d(K,I) = VY1(I,J,K) VZ1d(K,I) = VZ1(I,J,K) C1d(K,I) = C1(I,J,K) 200 continue do i=1,nl_i do k=1,nl_k RHO2d(K,I) = RHO2(I,J,K) VX2d(K,I) = VX2(I,J,K) VY2d(K,I) = VY2(I,J,K) VZ2d(K,I) = VZ2(I,J,K) C2d(K,I) = C2(I,J,K) BDVX2d(K,I) = BDVX2(I,J,K) BDVY2d(K,I) = BDVY2(I,J,K) BDVZ2d(K,I) = BDVZ2(I,J,K) enddo enddo C DO 220 I=1,NL_I DO 220 K=1,Nl_KP1 FACEd(K,I) = FACE_k(I,J,K) do ltran = 1,6 transd(k,i,ltran) = tran_k(i,j,k,ltran) enddo dtx = transd(k,i,2)*transd(k,i,6)-transd(k,i,3)*transd(k,i,5) dty = transd(k,i,3)*transd(k,i,4)-transd(k,i,1)*transd(k,i,6) dtz = transd(k,i,1)*transd(k,i,5)-transd(k,i,2)*transd(k,i,4) dinv = 1./sqrt(dtx**2+dty**2+dtz**2) transd(k,i,7) = dtx*dinv transd(k,i,8) = dty*dinv transd(k,i,9) = dtz*dinv BNORMd(K,I) = BK1(I,J,K) 220 CONTINUE C do 230 i=1,nl_i do 230 k=1,nl_kp1 Bnqd(k,i) = bnqfk(i,j,k) bsqqd(k,i) = bsqqk(i,j,k) bxqd(k,i) = bxqfk(i,j,k) byqd(k,i) = byqfk(i,j,k) bzqd(k,i) = bzqfk(i,j,k) bxqnd(k,i) = bxqnk(i,j,k) byqnd(k,i) = byqnk(i,j,k) bzqnd(k,i) = bzqnk(i,j,k) 230 continue C 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 DO 352 I=1,NL_IP1 DO 352 k=1,NL_kP1 tran2d(k,i,1) = vec_j(i,j,k,1) tran2d(k,i,2) = vec_j(i,j,k,2) tran2d(k,i,3) = vec_j(i,j,k,3) tran2d(k,i,7) = vec_j(i,j,k,7) tran2d(k,i,8) = vec_j(i,j,k,8) tran2d(k,i,9) = vec_j(i,j,k,9) tran2d(k,i,4) = vec_j(i,j,k,8)*vec_j(i,j,k,3) - $ vec_j(i,j,k,9)*vec_j(i,j,k,2) tran2d(k,i,5) = vec_j(i,j,k,9)*vec_j(i,j,k,1) - $ vec_j(i,j,k,7)*vec_j(i,j,k,3) tran2d(k,i,6) = vec_j(i,j,k,7)*vec_j(i,j,k,2) - $ vec_j(i,j,k,8)*vec_j(i,j,k,1) dinv = 1./sqrt(tran2d(k,i,4)**2+tran2d(k,i,5)**2+ $ tran2d(k,i,6)**2 ) tran2d(k,i,4) = dinv*tran2d(k,i,4) tran2d(k,i,5) = dinv*tran2d(k,i,5) tran2d(k,i,6) = dinv*tran2d(k,i,6) e2(k,i) = ej_g(i,j,k) bound_mix(k) = bmix_ki(j,k) edged(k,i) = edge_j(i,j,k) 352 CONTINUE C do i=llow,nl_i+no2 do k=1,nl_kp1 xnormdi(k,i) =tran_k(i,j,k,1) ynormdi(k,i) =tran_k(i,j,k,2) znormdi(k,i) =tran_k(i,j,k,3) facedi(k,i) = face_k(i,j,k) enddo enddo C C do i=1,nl_ip1 do k=llow,nl_k+no2 xnormdj(k,i) =tran_i(i,j,k,1) ynormdj(k,i) =tran_i(i,j,k,2) znormdj(k,i) =tran_i(i,j,k,3) facedj(k,i) = face_i(i,j,k) enddo enddo C DO 320 I=0,NL_I+2 DO 320 K=0,Nl_K+2 XqD(K,I) = 0.5*(X(I,J,K)+X(I,J+1,K)) YqD(K,I) = 0.5*(Y(I,J,K)+Y(I,J+1,K)) ZqD(K,I) = 0.5*(Z(I,J,K)+Z(I,J+1,K)) 320 CONTINUE C DO 340 I=1,NL_I DO 340 K=1,NL_KP1 BI2d(K,I) = BK2(I,J,K) 340 CONTINUE do i=-no2+1,nl_i+no2 do k=1,nl_k+1 BID(K,I) = Bk(I,J,K) BI1d(K,I) = BK1(I,J,K) enddo enddo C DO 350 I=1,NL_IP1 DO 350 K=1,NL_K BJ2d(K,I) = BI2(I,J,K) 350 CONTINUE do i=1,nl_i+1 do k=-no2+1,nl_k+no2 bjd(k,i) = bi(i,j,k) bj1d(k,i) = bi1(i,j,k) enddo enddo C do 356 i=1,nl_ip1 do 356 k=1,nl_kp1 bx0d(k,i) = bxqj(i,j,k) by0d(k,i) = byqj(i,j,k) bz0d(k,i) = bzqj(i,j,k) 356 continue C C C * WRITE (9,*) '************** J=', J,' **************************' * CALL WROTE('TRANS(1)',TRAN2d(1,1,1),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(2)',TRAN2d(1,1,2),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(3)',TRAN2d(1,1,3),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(4)',TRAN2d(1,1,4),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(5)',TRAN2d(1,1,5),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(6)',TRAN2d(1,1,6),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(7)',TRAN2d(1,1,7),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(8)',TRAN2d(1,1,8),1,NIP1,NIP1,1,nkp1,nkp1) * CALL WROTE('TRANS(9)',TRAN2d(1,1,9),1,NIP1,NIP1,1,nkp1,nkp1) C * write(9,*) '**************************************************' * write(9,*) 'KIMHD, j= ',j * 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 C DO 500 I=1,NL_I DO 500 K=1,NL_K RHO2(I,J,K) = RHO2d(K,I) VX2(I,J,K) = VX2d(K,I) VY2(I,J,K) = VY2d(K,I) VZ2(I,J,K) = VZ2d(K,I) C2(I,J,K) = C2d(K,I) BDVX2(I,J,K) = BDVX2d(K,I) BDVY2(I,J,K) = BDVY2d(K,I) 500 BDVZ2(I,J,K) = BDVZ2d(K,I) C C DO 520 I=1,NL_I DO 520 K=1,NL_KP1 520 Bk2(I,J,K) = Bi2d(K,I) DO 530 I=1,NL_IP1 DO 530 K=1,Nl_K 530 BI2(I,J,K) = Bj2d(K,I) C if ( ilow_bound) then lower_i = 2 else lower_i =1 endif do i=lower_i,nl_i+1 do k=1,nl_k+1 ej_g(i,j,k) = e2(k,i) enddo enddo C C 999 CONTINUE C C C C NOW RETURN TO MHD C C RETURN END