Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/hydrostatic_eq_WL

00001 subroutine hydrostatic_eq_WL_MHD(emd,utf,uxf,uyf,uzf)
00002   use phys_constant, only  : long
00003   use grid_parameter
00004   use coordinate_grav_r, only : rg
00005   use trigonometry_grav_phi, only   : sinphig, cosphig
00006 !  use def_matter, only : utdf, uxdf, uydf, uzdf
00007   use def_matter_parameter, only : ome, ber
00008   use def_metric_on_SFC_CF
00009   use def_metric_on_SFC_WL
00010   use def_emfield, only : vaxd, vayd, vazd
00011   use def_vector_phi, only : vec_phif 
00012   use integrability_fnc_MHD
00013   use make_array_3d
00014   use interface_interpo_gr2fl
00015   use interface_flgrad_4th_gridpoint
00016   implicit none
00017   real(long), pointer :: emd(:,:,:)
00018   real(long), pointer :: utf(:,:,:), uxf(:,:,:), uyf(:,:,:), uzf(:,:,:)
00019   real(long), pointer :: vaxdf(:,:,:), vaydf(:,:,:), vazdf(:,:,:)
00020   real(long), pointer :: vaphidf(:,:,:)
00021   real(long) :: vecphif(3)
00022   real(long) :: ovufc(3), ovdfc(3), bvdfc(3)
00023   real(long) :: Aphi, dxAphi, dyAphi, dzAphi
00024   real(long) :: Bphi, dxAx, dyAx, dzAx, dxAz, dyAz, dzAz
00025   real(long) :: hh, ut, ux, uy, uz, vx, vy, vz, pre, rho, ene, qq
00026   real(long) :: uphid, vphi, utd, zfac
00027   real(long) :: gmxxdf, gmxydf, gmxzdf, gmyxdf, gmyydf, gmyzdf, 
00028                gmzxdf, gmzydf, gmzzdf, ovovf,  alphff, psiff
00029   integer    :: irf, itf, ipf
00030 !
00031   call alloc_array3d(vaxdf, 0, nrf, 0, ntf, 0, npf)
00032   call alloc_array3d(vaydf, 0, nrf, 0, ntf, 0, npf)
00033   call alloc_array3d(vazdf, 0, nrf, 0, ntf, 0, npf)
00034   call alloc_array3d(vaphidf, 0, nrf, 0, ntf, 0, npf)
00035 !
00036   call interpo_gr2fl(vaxd, vaxdf)
00037   call interpo_gr2fl(vayd, vaydf)
00038   call interpo_gr2fl(vazd, vazdf)
00039   vaphidf(0:nrf,0:ntf,0:npf) = vaydf(0:nrf,0:ntf,0:npf) &
00040   &                        *vec_phif(0:nrf,0:ntf,0:npf,2)
00041 !
00042 ! Solve on phi=0 plane
00043   ipf = 0
00044   do itf = 0, ntf
00045     do irf = 0, nrf
00046 !
00047       psiff  = psif(irf,itf,ipf)
00048       alphff = alphf(irf,itf,ipf)
00049       gmxxdf = 1.0d0 + hxxdf(irf,itf,ipf)
00050       gmxydf =         hxydf(irf,itf,ipf)
00051       gmxzdf =         hxzdf(irf,itf,ipf)
00052       gmyydf = 1.0d0 + hyydf(irf,itf,ipf)
00053       gmyzdf =         hyzdf(irf,itf,ipf)
00054       gmzzdf = 1.0d0 + hzzdf(irf,itf,ipf)
00055       gmyxdf = gmxydf
00056       gmzxdf = gmxzdf
00057       gmzydf = gmyzdf
00058       Aphi = vaphidf(irf,itf,ipf)
00059       vecphif(1)= vec_phif(irf,itf,ipf,1)
00060       vecphif(2)= vec_phif(irf,itf,ipf,2)
00061       vecphif(3)= vec_phif(irf,itf,ipf,3)
00062       call calc_integrability_fnc_MHD(Aphi)
00063       call flgrad_4th_gridpoint(vaphidf,dxAphi,dyAphi,dzAphi,irf,itf,ipf)
00064       call flgrad_4th_gridpoint(vaxdf,dxAx,dyAx,dzAx,irf,itf,ipf)
00065       call flgrad_4th_gridpoint(vazdf,dxAz,dyAz,dzAz,irf,itf,ipf)
00066       qq = emd(irf,itf,ipf)
00067       zfac = 1.0d0
00068       if (qq <= 1.0d-14) zfac = 0.0
00069       call peos_q2hprho(qq, hh, pre, rho, ene)
00070       Bphi = -(dxAz - dzAx)
00071 !
00072 ! meridional flow
00073       uxf(irf,itf,ipf) = 1.0d0/(rho*alphff*psiff**6)*MHDfnc_dPSI*(-dzAphi)
00074       uzf(irf,itf,ipf) = 1.0d0/(rho*alphff*psiff**6)*MHDfnc_dPSI*(+dxAphi)
00075       uxf(irf,itf,ipf) = uxf(irf,itf,ipf)*zfac
00076       uzf(irf,itf,ipf) = uzf(irf,itf,ipf)*zfac
00077 !
00078       ut = utf(irf,itf,ipf)
00079       ux = uxf(irf,itf,ipf)
00080       uy = uyf(irf,itf,ipf)
00081       uz = uzf(irf,itf,ipf)
00082       vx = ux/ut
00083       vy = uy/ut
00084       vz = uz/ut
00085       ovufc(1) = bvxuf(irf,itf,ipf) + vx
00086       ovufc(2) = bvyuf(irf,itf,ipf) + vy
00087       ovufc(3) = bvzuf(irf,itf,ipf) + vz
00088       ovdfc(1) = bvxdf(irf,itf,ipf) + gmxxdf*vx + gmxydf*vy + gmxzdf*vz
00089       ovdfc(2) = bvydf(irf,itf,ipf) + gmyxdf*vx + gmyydf*vy + gmyzdf*vz
00090       ovdfc(3) = bvzdf(irf,itf,ipf) + gmzxdf*vx + gmzydf*vy + gmzzdf*vz
00091       ovovf = ovdfc(1)*ovufc(1) + ovdfc(2)*ovufc(2) + ovdfc(3)*ovufc(3)
00092 !
00093 !write(6,*) ut
00094 !write(6,*) ux, uy, uz
00095 !write(6,*) rho, MHDfnc_dPSI
00096 !write(6,*) dxAphi, dzAphi
00097 !!write(6,*) alphff, psiff, ovovf
00098 !!write(6,*) ut, psiff, ovdfc(2),vecphif(2)
00099 !!write(6,*) Aphi, MHDfnc_Lambda, MHDfnc_dAt*uphid, hh
00100 !stop
00101 !
00102 ! ut 
00103       utf(irf,itf,ipf) = 1.0d0/sqrt(alphff**2 - psiff**4*ovovf)
00104 !
00105 !if (irf.eq.nrf.and.itf.eq.ntf/2) then
00106 !write(6,'(a3,2x,1p,3e15.7)') 'cmp', alphff,psiff, ovovf
00107 !write(6,'(a3,2x,1p,3e15.7)') 'cmp', ovufc(1), ovufc(2), ovufc(3)
00108 !write(6,'(a3,2x,1p,3e15.7)') 'cmp', ovdfc(1), ovdfc(2), ovdfc(3)
00109 !write(6,'(a3,2x,1p,3e15.7)') 'cmp', vx, vy, vz
00110 !write(6,'(a3,2x,1p,3e15.7)') 'cmp', bvxuf(irf,itf,ipf), &
00111 !                bvyuf(irf,itf,ipf), bvzuf(irf,itf,ipf)
00112 !end if
00113 !!!
00114       ut = utf(irf,itf,ipf)
00115       ux = uxf(irf,itf,ipf)
00116       uy = uyf(irf,itf,ipf)
00117       uz = uzf(irf,itf,ipf)
00118       vx = ux/ut
00119       vy = uy/ut
00120       vz = uz/ut
00121       bvdfc(1) = bvxdf(irf,itf,ipf)
00122       bvdfc(2) = bvydf(irf,itf,ipf)
00123       bvdfc(3) = bvzdf(irf,itf,ipf)
00124 !
00125 ! toroidal flow
00126 !WIND      vphi = MHDfnc_dPSI*Bphi/(rho*ut*alphff*psiff**6) - MHDfnc_dAt
00127       vphi = - MHDfnc_dAt
00128       uyf(irf,itf,ipf) = ut*vphi*vecphif(2)
00129 !
00130 ! WIND! enthalpy to emden fnc
00131 ! WIND      utd = ut*(-alphff**2 + psiff**4*(bvdfc(1)*ovufc(1) &
00132 ! WIND      &          + bvdfc(2)*ovufc(2) + bvdfc(3)*ovufc(3)))
00133 ! WIND      uphid = ut*psiff**4*ovdfc(2)*vecphif(2)
00134 ! WIND      hh = MHDfnc_Lambda/(utd - MHDfnc_dAt*uphid)
00135 !
00136 ! For GS no meridional circulation
00137       hh = ber*ut*dexp(-MHDfnc_Lambda_GS)
00138 !
00139       call peos_h2qprho(hh, qq, pre, rho, ene)
00140       emd(irf,itf,ipf) = qq
00141 !
00142 !if (itf.eq.ntfeq-1.or.itf.eq.ntfeq) then 
00143 !if (irf.ge.nrf-1) then 
00144 !write(6,*)ut,alphff, + psiff,bvdfc(1)*ovufc(1), &
00145 !      &          + bvdfc(2)*ovufc(2), + bvdfc(3)*ovufc(3),&
00146 !      & ovdfc(2),vecphif(2), qq
00147 !end if
00148 !end if
00149 !
00150     end do
00151   end do
00152 !
00153 !write(6,'(a3,2x,1p,3e15.7)') 'uxf', uxf(nrf*3/4,ntf/2,0)
00154 !write(6,'(a3,2x,1p,3e15.7)') 'uzf', uzf(nrf*3/4,ntf/2,0)
00155 !write(6,'(a3,2x,1p,3e15.7)') 'uxf', uxf(nrf,0,0), uxf(nrf,ntf/2,0), uxf(0,0,0)
00156 !write(6,'(a3,2x,1p,3e15.7)') 'uyf', uyf(nrf,0,0), uyf(nrf,ntf/2,0), uyf(0,0,0)
00157 !write(6,'(a3,2x,1p,3e15.7)') 'uzf', uzf(nrf,0,0), uzf(nrf,ntf/2,0), uzf(0,0,0)
00158 !
00159 ! Copy to phi /= 0 planes.
00160   do ipf = 1, npf
00161     do itf = 0, ntf
00162       do irf = 0, nrf
00163         emd(irf,itf,ipf) = emd(irf,itf,0)
00164         utf(irf,itf,ipf) = utf(irf,itf,0)
00165         uxf(irf,itf,ipf) = cosphig(ipf)*uxf(irf,itf,0) &
00166         &                - sinphig(ipf)*uyf(irf,itf,0)
00167         uyf(irf,itf,ipf) = sinphig(ipf)*uxf(irf,itf,0) &
00168         &                + cosphig(ipf)*uyf(irf,itf,0)
00169         uzf(irf,itf,ipf) = uzf(irf,itf,0)
00170       end do
00171     end do
00172   end do
00173 !
00174   deallocate(vaxdf)
00175   deallocate(vaydf)
00176   deallocate(vazdf)
00177   deallocate(vaphidf)
00178 !
00179 end subroutine hydrostatic_eq_WL_MHD
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/hydrostatic_eq_WL_MHD.f90
