00001 subroutine current_MHD_on_SFC
00002 use phys_constant, only : long, pi
00003 use grid_parameter
00004 use coordinate_grav_r, only : rg
00005 use trigonometry_grav_phi, only : sinphig, cosphig
00006 use def_matter, only : emd, hhf, utf, uxdf, uydf, uzdf
00007 use def_matter_parameter, only : ome, ber, radi
00008 use def_metric_on_SFC_CF, only : alphf, psif, bvxuf, bvzuf
00009 use def_emfield, only : vaxd, vayd, vazd, jtuf, jxuf, jyuf, jzuf
00010 use def_faraday_tensor, only : fiju_grid, fxu_grid, fzu_grid
00011 use def_vector_x, only : vec_xf
00012 use def_vector_phi, only : vec_phif
00013 use integrability_fnc_MHD
00014 use make_array_3d
00015 use interface_interpo_gr2fl
00016 use interface_flgrad_4th_gridpoint
00017 implicit none
00018 real(long), pointer :: vaphidf(:,:,:)
00019 real(long), pointer :: vaxdf(:,:,:), vaydf(:,:,:), vazdf(:,:,:)
00020 real(long), pointer :: huxdf(:,:,:), huzdf(:,:,:), huphidf(:,:,:)
00021 real(long), pointer :: fxzug(:,:,:), fxzuf(:,:,:), bmf(:,:,:)
00022 real(long), pointer :: fxuf(:,:,:), fzuf(:,:,:)
00023 real(long) :: vecphif(3)
00024 real(long) :: Aphi, dxAphi, dyAphi, dzAphi
00025 real(long) :: Bphi, dxAx, dyAx, dzAx, dxAz, dyAz, dzAz, BB
00026 real(long) :: vorphi, dxhux, dyhux, dzhux, dxhuz, dyhuz, dzhuz
00027 real(long) :: huphi, dxhuphi, dyhuphi, dzhuphi
00028 real(long) :: dxbmf, dybmf, dzbmf, pi4inv = 0.25/pi
00029 real(long) :: hh, ut, ux, uy, uz, vx, vy, vz, pre, rho, ene, qq
00030 real(long) :: alphff, psiff, alps6f, jphiuf, jt, f11inv, radi2inv
00031 integer :: irf, itf, ipf
00032
00033 call alloc_array3d(vaxdf, 0, nrf, 0, ntf, 0, npf)
00034 call alloc_array3d(vaydf, 0, nrf, 0, ntf, 0, npf)
00035 call alloc_array3d(vazdf, 0, nrf, 0, ntf, 0, npf)
00036 call alloc_array3d(vaphidf, 0, nrf, 0, ntf, 0, npf)
00037 call alloc_array3d(huxdf, 0, nrf, 0, ntf, 0, npf)
00038 call alloc_array3d(huzdf, 0, nrf, 0, ntf, 0, npf)
00039 call alloc_array3d(huphidf, 0, nrf, 0, ntf, 0, npf)
00040 call alloc_array3d(bmf, 0, nrf, 0, ntf, 0, npf)
00041 call alloc_array3d(fxzuf, 0, nrf, 0, ntf, 0, npf)
00042 call alloc_array3d(fxzug, 0, nrg, 0, ntg, 0, npg)
00043 call alloc_array3d(fxuf, 0, nrf, 0, ntf, 0, npf)
00044 call alloc_array3d(fzuf, 0, nrg, 0, ntg, 0, npg)
00045
00046 call interpo_gr2fl(vaxd, vaxdf)
00047 call interpo_gr2fl(vayd, vaydf)
00048 call interpo_gr2fl(vazd, vazdf)
00049 vaphidf(0:nrf,0:ntf,0:npf)=vaydf(0:nrf,0:ntf,0:npf) &
00050 & *vec_phif(0:nrf,0:ntf,0:npf,2)
00051 huxdf(0:nrf,0:ntf,0:npf) = hhf(0:nrf,0:ntf,0:npf)*uxdf(0:nrf,0:ntf,0:npf)
00052 huzdf(0:nrf,0:ntf,0:npf) = hhf(0:nrf,0:ntf,0:npf)*uzdf(0:nrf,0:ntf,0:npf)
00053 huphidf(0:nrf,0:ntf,0:npf)=hhf(0:nrf,0:ntf,0:npf)*uydf(0:nrf,0:ntf,0:npf) &
00054 & *vec_phif(0:nrf,0:ntf,0:npf,2)
00055
00056 fxzug(0:nrg,0:ntg,0:npg) = fiju_grid(0:nrg,0:ntg,0:npg,2)
00057 call interpo_gr2fl(fxzug, fxzuf)
00058 call interpo_gr2fl(fxu_grid, fxuf)
00059 call interpo_gr2fl(fzu_grid, fzuf)
00060 bmf(0:nrf,0:ntf,0:npf)=(-fxzuf(0:nrf,0:ntf,0:npf) &
00061 & +(bvxuf(0:nrf,0:ntf,0:npf)*fzuf(0:nrf,0:ntf,0:npf) &
00062 & - bvzuf(0:nrf,0:ntf,0:npf)*fxuf(0:nrf,0:ntf,0:npf)) &
00063 & /alphf(0:nrf,0:ntf,0:npf)) &
00064 & *alphf(0:nrf,0:ntf,0:npf)*psif(0:nrf,0:ntf,0:npf)**6 &
00065 & *vec_xf(0:nrf,0:ntf,0:npf,1)
00066
00067
00068 ipf = 0
00069 do itf = 0, ntf
00070 do irf = 0, nrf
00071
00072 psiff = psif(irf,itf,ipf)
00073 alphff = alphf(irf,itf,ipf)
00074 alps6f = alphff*psiff**6
00075 Aphi = vaphidf(irf,itf,ipf)
00076 huphi = huphidf(irf,itf,ipf)
00077 ut = utf(irf,itf,ipf)
00078 jt = jtuf(irf,itf,ipf)
00079 vecphif(1)= vec_phif(irf,itf,ipf,1)
00080 vecphif(2)= vec_phif(irf,itf,ipf,2)
00081 vecphif(3)= vec_phif(irf,itf,ipf,3)
00082 call calc_integrability_fnc_MHD(Aphi)
00083 call flgrad_4th_gridpoint(vaphidf,dxAphi,dyAphi,dzAphi,irf,itf,ipf)
00084 call flgrad_4th_gridpoint(huphidf,dxhuphi,dyhuphi,dzhuphi,irf,itf,ipf)
00085 call flgrad_4th_gridpoint(vaxdf,dxAx,dyAx,dzAx,irf,itf,ipf)
00086 call flgrad_4th_gridpoint(vazdf,dxAz,dyAz,dzAz,irf,itf,ipf)
00087 Bphi = -(dxAz - dzAx)
00088 call flgrad_4th_gridpoint(huxdf,dxhux,dyhux,dzhux,irf,itf,ipf)
00089 call flgrad_4th_gridpoint(huzdf,dxhuz,dyhuz,dzhuz,irf,itf,ipf)
00090 vorphi = dxhuz - dzhux
00091 call flgrad_4th_gridpoint(bmf,dxbmf,dybmf,dzbmf,irf,itf,ipf)
00092 qq = emd(irf,itf,ipf)
00093 call peos_q2hprho(qq, hh, pre, rho, ene)
00094
00095
00096
00097 f11inv = 0.0d0
00098 if (irf.ne.0.and.itf.ne.0.and.itf.ne.ntf) &
00099 & f11inv = 1.0d0/vec_xf(irf,itf,ipf,1)
00100
00101
00102 radi2inv = 1.0d0/radi**2
00103
00104
00105
00106
00107
00108
00109
00110
00111 jxuf(irf,itf,ipf) = radi2inv*f11inv/alps6f &
00112 & * MHDfnc_dLambda_phi*(-dzAphi)
00113
00114 jzuf(irf,itf,ipf) = radi2inv*f11inv/alps6f &
00115 & * MHDfnc_dLambda_phi*(+dxAphi)
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 jphiuf = (radi2inv*f11inv/alps6f)*MHDfnc_dLambda_phi*Bphi &
00129 & - rho*hh*MHDfnc_dLambda_GS &
00130 & - MHDfnc_dAt*jt
00131 jyuf(irf,itf,ipf) = jphiuf*vecphif(2)
00132
00133
00134
00135
00136
00137 end do
00138 end do
00139
00140
00141
00142
00143
00144 do ipf = 1, npf
00145 do itf = 0, ntf
00146 do irf = 0, nrf
00147 jtuf(irf,itf,ipf) = jtuf(irf,itf,0)
00148 jxuf(irf,itf,ipf) = cosphig(ipf)*jxuf(irf,itf,0) &
00149 & - sinphig(ipf)*jyuf(irf,itf,0)
00150 jyuf(irf,itf,ipf) = sinphig(ipf)*jxuf(irf,itf,0) &
00151 & + cosphig(ipf)*jyuf(irf,itf,0)
00152 jzuf(irf,itf,ipf) = jzuf(irf,itf,0)
00153 end do
00154 end do
00155 end do
00156
00157 itf = ntgeq; ipf = 0
00158
00159 open(15,file='test_vec_cur',status='unknown')
00160 do irf = 0, nrf
00161 write(15,'(1p,9e20.12)') rg(irf), jtuf(irf,itf,ipf) &
00162 & , jxuf(irf,itf,ipf) &
00163 & , jyuf(irf,itf,ipf) &
00164 & , jzuf(irf,itf,ipf)
00165 end do
00166 close(15)
00167
00168 deallocate(vaxdf)
00169 deallocate(vaydf)
00170 deallocate(vazdf)
00171 deallocate(vaphidf)
00172 deallocate(huxdf)
00173 deallocate(huzdf)
00174 deallocate(huphidf)
00175 deallocate(bmf)
00176 deallocate(fxzuf)
00177 deallocate(fxzug)
00178 deallocate(fxuf)
00179 deallocate(fzuf)
00180
00181 end subroutine current_MHD_on_SFC
