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

00001 subroutine update_parameter_WL_MHD(convf)
00002   use phys_constant,  only  : long
00003   use grid_parameter, only : nrg, ntg, npg, nrf_deform, nrf, ntf, npf, ntgxy
00004   use def_metric,     only : psi, alph, bvxu, bvyu, bvzu
00005   use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd
00006   use def_matter,     only : utf, uxf, uyf, uzf
00007   use def_emfield,    only : vayd
00008   use def_matter_parameter
00009   use coordinate_grav_r, only : rg
00010   use trigonometry_grav_theta, only : sinthg, costhg
00011   use trigonometry_grav_phi, only   : sinphig, cosphig
00012   use def_vector_phi, only : vec_phig, vec_phif
00013   use integrability_fnc_MHD, only : MHDfnc_dAt, MHDfnc_Lambda
00014   use interface_minv
00015   implicit none
00016   integer, parameter :: neq = 3, ndim = 3
00017   real(long), intent(in) :: convf
00018   real(long) ::    sgg(neq+1),  gg(neq+1,neq+1)
00019   real(long) ::  log_a(neq),  psi2al(neq)
00020   real(long) ::     ve(neq,ndim), be(neq,ndim)
00021   real(long) ::   vphi(neq,ndim), gm(neq,ndim,ndim)
00022   real(long) :: MHDfnc_dAt_val(neq), MHDfnc_Lambda_val(neq)
00023   real(long) :: ovufc(neq,ndim), oAufc(neq,ndim)
00024   real(long) :: dovufcii, dovufcjj, doAufcjj
00025   real(long) :: gmovoA(neq), gmovov(neq), gmdovov(neq), gmdovoA(neq)
00026   real(long) :: term1(neq), dterm1do(neq), dterm1dr(neq)
00027   real(long) :: term2(neq), dterm2do(neq), dterm2dr(neq)
00028   real(long) :: log_h(neq), loghc, Aphi
00029   real(long) :: radiold, berold, omeold
00030   real(long) :: ddradi, ddber, ddome
00031   real(long) :: facfac, numero, error, hc, pre, rho, ene
00032   integer    :: nin, irg, irf, itf, ipf, it, ip, ia, ii, jj
00033 !##############################################
00034   nin = nrf_deform
00035 !
00036   sgg = 0.0d0
00037   gg  = 0.0d0
00038 !
00039   call peos_q2hprho(emdc, hc, pre, rho, ene)
00040   loghc = dlog(hc)
00041   omeold = ome
00042   berold = ber
00043   radiold = radi
00044   numero = 1
00045 !
00046   do ia = 1, 3    
00047 !   ia = 1 pole, ia = 2 equator, ia = 3 center.
00048     if (ia.eq.1) then; irg = nin; irf = nrf; it = 0    ; ip = 0; end if
00049     if (ia.eq.2) then; irg = nrf; irf = nrf; it = ntgxy; ip = 0; end if
00050     if (ia.eq.3) then; irg = 0  ; irf = 0  ; it = 0    ; ip = 0; end if
00051     log_a(ia)  = dlog(alph(irg,it,ip)**(1.0d0/radi**2))
00052     psi2al(ia) = (psi(irg,it,ip)**2/alph(irg,it,ip))**(1.0d0/radi**2)
00053     be(ia,1) = bvxu(irg,it,ip)
00054     be(ia,2) = bvyu(irg,it,ip)
00055     be(ia,3) = bvzu(irg,it,ip)
00056     vphi(ia,1) = vec_phig(irg,it,ip,1)
00057     vphi(ia,2) = vec_phig(irg,it,ip,2)
00058     vphi(ia,3) = vec_phig(irg,it,ip,3)
00059     gm(ia,1,1) = 1.0d0 + hxxd(irg,it,ip)
00060     gm(ia,1,2) =         hxyd(irg,it,ip)
00061     gm(ia,1,3) =         hxzd(irg,it,ip)
00062     gm(ia,2,2) = 1.0d0 + hyyd(irg,it,ip)
00063     gm(ia,2,3) =         hyzd(irg,it,ip)
00064     gm(ia,3,3) = 1.0d0 + hzzd(irg,it,ip)
00065     gm(ia,2,1) = gm(ia,1,2)
00066     gm(ia,3,1) = gm(ia,1,3)
00067     gm(ia,3,2) = gm(ia,2,3)
00068 !
00069     log_h(ia) = 0.0d0
00070     if (ia.eq.3) log_h(ia) = loghc
00071   end do
00072 ! ---
00073 !write (6,*) "test values"
00074 !write (6,'(1p,3e12.4)') ome, ber, radi
00075 !write (6,'(1p,3e12.4)') (log_a(ia), ia = 1, 3)
00076 !write (6,'(1p,3e12.4)') (psi2al(ia), ia = 1, 3)
00077 !write (6,'(1p,3e12.4)') (ve(ia,1), ia = 1, 3)
00078 !write (6,'(1p,3e12.4)') (ve(ia,2), ia = 1, 3)
00079 !write (6,'(1p,3e12.4)') (ve(ia,3), ia = 1, 3)
00080 !write (6,'(1p,3e12.4)') (be(ia,1), ia = 1, 3)
00081 !write (6,'(1p,3e12.4)') (be(ia,2), ia = 1, 3)
00082 !write (6,'(1p,3e12.4)') (be(ia,3), ia = 1, 3)
00083 !write (6,'(1p,3e12.4)') (vphi(ia,1), ia = 1, 3)
00084 !write (6,'(1p,3e12.4)') (vphi(ia,2), ia = 1, 3)
00085 !write (6,'(1p,3e12.4)') (vphi(ia,3), ia = 1, 3)
00086 ! write (6,'(1p,3e12.4)') (log_h(ia), ia = 1, 3)
00087 ! write (6,'(1p,3e12.4)') error
00088 !do ia = 1, 3
00089 ! write (6,'(1p,3e12.4)') (gm(ia,1,ii), ii = 1, 3)
00090 ! write (6,'(1p,3e12.4)') (gm(ia,2,ii), ii = 1, 3)
00091 ! write (6,'(1p,3e12.4)') (gm(ia,3,ii), ii = 1, 3)
00092 !end do
00093 !!!!stop
00094 ! ---
00095   do
00096     do ia = 1, 3 ! ia = pole, equator, center
00097 !
00098       if (ia.eq.1) then; irg = nin; irf = nrf; it = 0    ; ip = 0; end if
00099       if (ia.eq.2) then; irg = nrf; irf = nrf; it = ntgxy; ip = 0; end if
00100       if (ia.eq.3) then; irg = 0  ; irf = 0  ; it = 0    ; ip = 0; end if
00101       ve(ia,1) = uxf(irf,it,ip)/utf(irf,it,ip)
00102 !!      ve(ia,2) = uyf(irf,it,ip)/utf(irf,it,ip)
00103       ve(ia,2) = ome*vphi(ia,2) ! uniform rotation SEE THE END OF THIS ROUTINE
00104       ve(ia,3) = uzf(irf,it,ip)/utf(irf,it,ip)
00105       Aphi = vayd(irg,it,ip)*vec_phig(irg,it,ip,2)
00106       call calc_integrability_fnc_MHD(Aphi)
00107       MHDfnc_dAt_val(ia) = MHDfnc_dAt
00108       MHDfnc_Lambda_val(ia) = MHDfnc_Lambda
00109 !
00110       do ii = 1, 3 ! ii = x,y,z
00111         ovufc(ia,ii) = be(ia,ii) + ve(ia,ii)
00112         oAufc(ia,ii) = be(ia,ii) - MHDfnc_dAt_val(ia)*vphi(ia,ii)
00113       end do
00114       gmovoA(ia) = 0.0d0
00115       gmovov(ia) = 0.0d0
00116       gmdovoA(ia) = 0.0d0
00117       gmdovov(ia) = 0.0d0
00118       do jj = 1, 3
00119         do ii = 1, 3
00120           gmovoA(ia) = gmovoA(ia) + gm(ia,ii,jj)*ovufc(ia,ii)*oAufc(ia,jj)
00121           gmovov(ia) = gmovov(ia) + gm(ia,ii,jj)*ovufc(ia,ii)*ovufc(ia,jj)
00122           dovufcii = 0.0d0
00123           dovufcjj = 0.0d0
00124           doAufcjj = 0.0d0
00125           if (ii.eq.2) dovufcii = vphi(ia,ii)          
00126           if (jj.eq.2) then 
00127             dovufcjj = vphi(ia,jj)          
00128             doAufcjj = (-MHDfnc_dAt_val(ia)/ome)*vphi(ia,jj)
00129           end if
00130           gmdovoA(ia) = gmdovoA(ia) &
00131           &           + gm(ia,ii,jj)*dovufcii*oAufc(ia,jj) &
00132           &           + gm(ia,ii,jj)*ovufc(ia,ii)*doAufcjj
00133           gmdovov(ia) = gmdovov(ia) &
00134           &           + gm(ia,ii,jj)*dovufcii*ovufc(ia,jj) &
00135           &           + gm(ia,ii,jj)*ovufc(ia,ii)*dovufcjj
00136         end do
00137       end do
00138 !
00139       term1(ia) = 1.0d0 - psi2al(ia)**(2.d0*radi**2)*gmovoA(ia)
00140       term2(ia) = 1.0d0 - psi2al(ia)**(2.d0*radi**2)*gmovov(ia)
00141 !
00142       dterm1do(ia) = - psi2al(ia)**(2.d0*radi**2)*2.0d0*gmdovoA(ia)
00143       dterm1dr(ia) = - 4.0d0*radi*dlog(psi2al(ia)) &
00144       &              * psi2al(ia)**(2.d0*radi**2)*gmovoA(ia)
00145       dterm2do(ia) = - psi2al(ia)**(2.d0*radi**2)*2.0d0*gmdovov(ia)
00146       dterm2dr(ia) = - 4.0d0*radi*dlog(psi2al(ia)) &
00147       &              * psi2al(ia)**(2.d0*radi**2)*gmovov(ia)
00148 !
00149       sgg(ia) = -(log_h(ia) + radi**2*log_a(ia) + dlog(term1(ia)) &
00150       &           - dlog(-MHDfnc_Lambda_val(ia))- 0.5d0*dlog(term2(ia))) 
00151 !
00152       gg(ia,1) =   dterm1do(ia)/term1(ia) - 0.5d0*dterm2do(ia)/term2(ia) 
00153       gg(ia,2) = - 1.0d0/(-MHDfnc_Lambda_val(ia))
00154       gg(ia,3) = 2.0d0*radi*log_a(ia) + dterm1dr(ia)/term1(ia) &
00155       &                          -0.5d0*dterm2dr(ia)/term2(ia)
00156 !
00157     end do
00158 !    test
00159 !if (numero.eq.1) then
00160 !!ia  = 1
00161 !ia  = 2
00162 !!ia  = 3
00163 ! write(6,*)log_h(ia), radi, log_a(ia) , term1(ia) , &
00164 !     &           -MHDfnc_Lambda_val(ia),term2(ia)
00165 ! write (6,'(1p,3e12.4)') sgg(1), sgg(2), sgg(3)
00166 ! write (6,'(1p,3e12.4)') gg(1,1), gg(1,2), gg(1,3)
00167 ! write (6,'(1p,3e12.4)') gg(2,1), gg(2,2), gg(2,3)
00168 ! write (6,'(1p,3e12.4)') gg(3,1), gg(3,2), gg(3,3)
00169 !end if
00170 !!!!stop
00171 !!pause
00172 !    test
00173     call minv(gg,sgg,neq,neq+1)
00174     ddome = sgg(1)
00175     ddber = sgg(2)
00176     ddradi = sgg(3)
00177     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00178     ome = ome + ddome*facfac
00179     ber = ber + ddber*facfac
00180     radi = radi + ddradi*facfac
00181     error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00182     numero = numero + 1
00183     if (numero>1000) then
00184       write(6,*)' numero = ', numero, '   error =',error
00185     end if
00186 !      if (iter<10.and.numero>10) exit
00187     if (numero>1010) exit
00188     if (dabs(error)<1.0d-08) exit
00189   end do
00190 !    test
00191 ! write (6,*) "test values"
00192 ! write (6,*) sgg(1), sgg(2), sgg(3)
00193 ! write (6,*) ome
00194 !    test
00195 !      write(6,*)' numero = ', numero, '  error = ',error
00196 !
00197   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00198   if (radi < 0.d0) radi = - radi
00199   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00200   if (ome < 0.d0) ome = - ome
00201   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00202 !  if (ber < 0.d0) ber = - ber
00203   ome = convf*ome + (1.d0-convf)*omeold
00204   ber = convf*ber + (1.d0-convf)*berold
00205   radi = convf*radi + (1.d0-convf)*radiold
00206 !
00207 ! write (6,*) 'parameters', ome, ber, radi
00208 ! --  Improving alpha and psi.  
00209 !
00210   alph = alph**((radi/radiold)**2)
00211   psi = psi**((radi/radiold)**2)
00212 !
00213 ! modify on phi=0 plane
00214   ipf = 0
00215   do itf = 0, ntf
00216     do irf = 0, nrf
00217       uyf(irf,itf,ipf) = utf(irf,itf,ipf)*ome*vec_phif(irf,itf,ipf,2)
00218     end do
00219   end do
00220 ! Copy to phi /= 0 planes.
00221   do ipf = 1, npf
00222     do itf = 0, ntf
00223       do irf = 0, nrf
00224         uxf(irf,itf,ipf) = cosphig(ipf)*uxf(irf,itf,0) &
00225         &                - sinphig(ipf)*uyf(irf,itf,0)
00226         uyf(irf,itf,ipf) = sinphig(ipf)*uxf(irf,itf,0) &
00227         &                + cosphig(ipf)*uyf(irf,itf,0)
00228       end do
00229     end do
00230   end do
00231 !
00232 end subroutine update_parameter_WL_MHD
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_WL_MHD.f90
