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

00001 subroutine update_parameter_WL_MHD_GS(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, omef
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_Lambda_GS
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 ! This routine can be used only for uniform rotation
00035 !#####################################################
00036   nin = nrf_deform
00037 !
00038   sgg = 0.0d0
00039   gg  = 0.0d0
00040 !
00041   call peos_q2hprho(emdc, hc, pre, rho, ene)
00042   loghc = dlog(hc)
00043   omeold = ome
00044   berold = ber
00045   radiold = radi
00046   numero = 1
00047 !
00048   do ia = 1, 3    
00049 !   ia = 1 pole, ia = 2 equator, ia = 3 center.
00050     if (ia.eq.1) then; irg = nin; irf = nrf; it = 0    ; ip = 0; end if
00051     if (ia.eq.2) then; irg = nrf; irf = nrf; it = ntgxy; ip = 0; end if
00052 !!    if (ia.eq.3) then; irg = 0  ; irf = 0  ; it = 0    ; ip = 0; end if
00053     if (ia.eq.3) then; irg = 0  ; irf = 0  ; it = ntgxy; ip = 0; end if
00054     log_a(ia)  = dlog(alph(irg,it,ip)**(1.0d0/radi**2))
00055     psi2al(ia) = (psi(irg,it,ip)**2/alph(irg,it,ip))**(1.0d0/radi**2)
00056     be(ia,1) = bvxu(irg,it,ip)
00057     be(ia,2) = bvyu(irg,it,ip)
00058     be(ia,3) = bvzu(irg,it,ip)
00059     vphi(ia,1) = vec_phig(irg,it,ip,1)
00060     vphi(ia,2) = vec_phig(irg,it,ip,2)
00061     vphi(ia,3) = vec_phig(irg,it,ip,3)
00062     gm(ia,1,1) = 1.0d0 + hxxd(irg,it,ip)
00063     gm(ia,1,2) =         hxyd(irg,it,ip)
00064     gm(ia,1,3) =         hxzd(irg,it,ip)
00065     gm(ia,2,2) = 1.0d0 + hyyd(irg,it,ip)
00066     gm(ia,2,3) =         hyzd(irg,it,ip)
00067     gm(ia,3,3) = 1.0d0 + hzzd(irg,it,ip)
00068     gm(ia,2,1) = gm(ia,1,2)
00069     gm(ia,3,1) = gm(ia,1,3)
00070     gm(ia,3,2) = gm(ia,2,3)
00071 !
00072     log_h(ia) = 0.0d0
00073     if (ia.eq.3) log_h(ia) = loghc
00074   end do
00075 ! ---
00076 !write (6,*) "test values"
00077 !write (6,'(1p,3e12.4)') ome, ber, radi
00078 !write (6,'(1p,3e12.4)') (log_a(ia), ia = 1, 3)
00079 !write (6,'(1p,3e12.4)') (psi2al(ia), ia = 1, 3)
00080 !write (6,'(1p,3e12.4)') (ve(ia,1), ia = 1, 3)
00081 !write (6,'(1p,3e12.4)') (ve(ia,2), ia = 1, 3)
00082 !write (6,'(1p,3e12.4)') (ve(ia,3), ia = 1, 3)
00083 !write (6,'(1p,3e12.4)') (be(ia,1), ia = 1, 3)
00084 !write (6,'(1p,3e12.4)') (be(ia,2), ia = 1, 3)
00085 !write (6,'(1p,3e12.4)') (be(ia,3), ia = 1, 3)
00086 !write (6,'(1p,3e12.4)') (vphi(ia,1), ia = 1, 3)
00087 !write (6,'(1p,3e12.4)') (vphi(ia,2), ia = 1, 3)
00088 !write (6,'(1p,3e12.4)') (vphi(ia,3), ia = 1, 3)
00089 ! write (6,'(1p,3e12.4)') (log_h(ia), ia = 1, 3)
00090 ! write (6,'(1p,3e12.4)') error
00091 !do ia = 1, 3
00092 ! write (6,'(1p,3e12.4)') (gm(ia,1,ii), ii = 1, 3)
00093 ! write (6,'(1p,3e12.4)') (gm(ia,2,ii), ii = 1, 3)
00094 ! write (6,'(1p,3e12.4)') (gm(ia,3,ii), ii = 1, 3)
00095 !end do
00096 !!!stop
00097 ! ---
00098   do
00099     do ia = 1, 3 ! ia = pole, equator, center
00100 !
00101       if (ia.eq.1) then; irg = nin; irf = nrf; it = 0    ; ip = 0; end if
00102       if (ia.eq.2) then; irg = nrf; irf = nrf; it = ntgxy; ip = 0; end if
00103 !!      if (ia.eq.3) then; irg = 0  ; irf = 0  ; it = 0    ; ip = 0; end if
00104       if (ia.eq.3) then; irg = 0  ; irf = 0  ; it = ntgxy; ip = 0; end if
00105       ve(ia,1) = uxf(irf,it,ip)/utf(irf,it,ip)
00106 !!      ve(ia,2) = uyf(irf,it,ip)/utf(irf,it,ip)
00107       ve(ia,2) = ome*vphi(ia,2) ! uniform rotation SEE THE END OF THIS ROUTINE
00108       ve(ia,3) = uzf(irf,it,ip)/utf(irf,it,ip)
00109       Aphi = vayd(irg,it,ip)*vec_phig(irg,it,ip,2)
00110       call calc_integrability_fnc_MHD(Aphi)
00111       MHDfnc_Lambda_val(ia) = MHDfnc_Lambda_GS
00112 !
00113       do ii = 1, 3 ! ii = x,y,z
00114         ovufc(ia,ii) = be(ia,ii) + ve(ia,ii)
00115       end do
00116       gmovov(ia) = 0.0d0
00117       gmdovov(ia) = 0.0d0
00118       do jj = 1, 3
00119         do ii = 1, 3
00120           gmovov(ia) = gmovov(ia) + gm(ia,ii,jj)*ovufc(ia,ii)*ovufc(ia,jj)
00121           dovufcii = 0.0d0
00122           dovufcjj = 0.0d0
00123           if (ii.eq.2) dovufcii = vphi(ia,ii)          
00124           if (jj.eq.2) dovufcjj = vphi(ia,jj)          
00125           gmdovov(ia) = gmdovov(ia) &
00126           &           + gm(ia,ii,jj)*dovufcii*ovufc(ia,jj) &
00127           &           + gm(ia,ii,jj)*ovufc(ia,ii)*dovufcjj
00128         end do
00129       end do
00130 !
00131       term2(ia) = 1.0d0 - psi2al(ia)**(2.d0*radi**2)*gmovov(ia)
00132 !
00133 !wrong   dterm2do(ia) = - psi2al(ia)**(2.d0*radi**2)*2.0d0*gmdovov(ia)
00134       dterm2do(ia) = - psi2al(ia)**(2.d0*radi**2)*gmdovov(ia)
00135       dterm2dr(ia) = - 4.0d0*radi*dlog(psi2al(ia)) &
00136       &              * psi2al(ia)**(2.d0*radi**2)*gmovov(ia)
00137 !
00138       sgg(ia) = -(log_h(ia) + radi**2*log_a(ia) + 0.5d0*dlog(term2(ia)) &
00139       &           + MHDfnc_Lambda_val(ia) - dlog(ber) ) 
00140 !
00141       gg(ia,1) =   0.5d0*dterm2do(ia)/term2(ia) 
00142       gg(ia,2) = - 1.0d0/ber
00143       gg(ia,3) = 2.0d0*radi*log_a(ia) + 0.5d0*dterm2dr(ia)/term2(ia)
00144 !
00145     end do
00146 !    test
00147 !if (numero.eq.1) then
00148 !!ia  = 1
00149 !ia  = 2
00150 !!ia  = 3
00151 ! write(6,*)log_h(ia), + radi**2, log_a(ia), + 0.5d0*dlog(term2(ia)), &
00152 !     &     MHDfnc_Lambda_val(ia), - dlog(ber), gmovov(ia)
00153 ! write (6,'(1p,3e12.4)') sgg(1), sgg(2), sgg(3)
00154 ! write (6,'(1p,3e12.4)') gg(1,1), gg(1,2), gg(1,3)
00155 ! write (6,'(1p,3e12.4)') gg(2,1), gg(2,2), gg(2,3)
00156 ! write (6,'(1p,3e12.4)') gg(3,1), gg(3,2), gg(3,3)
00157 ! write(6,*)'Aphi', Aphi,vayd(irg,it,ip)
00158 !end if
00159 !!!!stop
00160 !!pause
00161 !    test
00162     call minv(gg,sgg,neq,neq+1)
00163     ddome = sgg(1)
00164     ddber = sgg(2)
00165     ddradi = sgg(3)
00166     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00167     ome = ome + ddome*facfac
00168     ber = ber + ddber*facfac
00169     radi = radi + ddradi*facfac
00170     error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00171     numero = numero + 1
00172     if (numero>1000) then
00173       write(6,*)' numero = ', numero, '   error =',error
00174     end if
00175 !      if (iter<10.and.numero>10) exit
00176     if (numero>1010) exit
00177     if (dabs(error)<1.0d-08) exit
00178   end do
00179 !    test
00180 ! write (6,*) "test values"
00181 ! write (6,*) sgg(1), sgg(2), sgg(3)
00182 ! write (6,*) ome
00183 !    test
00184 !      write(6,*)' numero = ', numero, '  error = ',error
00185 !
00186   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00187   if (radi < 0.d0) radi = - radi
00188   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00189   if (ome < 0.d0) ome = - ome
00190   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00191 !  if (ber < 0.d0) ber = - ber
00192   ome = convf*ome + (1.d0-convf)*omeold
00193   ber = convf*ber + (1.d0-convf)*berold
00194   radi = convf*radi + (1.d0-convf)*radiold
00195 !
00196 ! write (6,*) 'parameters', ome, ber, radi
00197 ! --  Improving alpha and psi.  
00198 !
00199   alph = alph**((radi/radiold)**2)
00200   psi = psi**((radi/radiold)**2)
00201 !
00202 ! modify on phi=0 plane
00203   ipf = 0
00204   do itf = 0, ntf
00205     do irf = 0, nrf
00206       uxf(irf,itf,ipf) = 0.0d0
00207       uyf(irf,itf,ipf) = utf(irf,itf,ipf)*ome*vec_phif(irf,itf,ipf,2)
00208       uzf(irf,itf,ipf) = 0.0d0
00209       omef(irf,itf,ipf)= ome
00210     end do
00211   end do
00212 ! Copy to phi /= 0 planes.
00213   do ipf = 1, npf
00214     do itf = 0, ntf
00215       do irf = 0, nrf
00216         uxf(irf,itf,ipf) = cosphig(ipf)*uxf(irf,itf,0) &
00217         &                - sinphig(ipf)*uyf(irf,itf,0)
00218         uyf(irf,itf,ipf) = sinphig(ipf)*uxf(irf,itf,0) &
00219         &                + cosphig(ipf)*uyf(irf,itf,0)
00220         uzf(irf,itf,ipf) = uzf(irf,itf,0)
00221         omef(irf,itf,ipf)= omef(irf,itf,0)
00222       end do
00223     end do
00224   end do
00225 !!!
00226 end subroutine update_parameter_WL_MHD_GS
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_WL_MHD_GS.f90
