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

00001 subroutine iteration_WL_MHD(iter_count)
00002   use phys_constant, only :  long, nnrg, nntg, nnpg
00003   use def_formulation, only : chgra, chope
00004   use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd
00005   use def_emfield, only : va, alva, vaxd, vayd, vazd, jtu,jxu,jyu,jzu
00006   use grid_parameter
00007   use coordinate_grav_r
00008   use coordinate_grav_phi
00009   use coordinate_grav_theta
00010   use weight_midpoint_grav
00011   use make_array_2d
00012   use make_array_3d
00013   use make_array_4d
00014   use def_metric
00015   use def_matter
00016   use def_matter_parameter
00017   use interface_interpo_fl2gr
00018   use interface_poisson_solver
00019   use interface_poisson_solver_At
00020   use interface_poisson_solver_Az
00021   use interface_compute_shift
00022   use interface_compute_alps2alph
00023   use interface_compute_fnc_multiple
00024   use interface_compute_fnc_division
00025   use interface_hydrostatic_eq_WL_MHD
00026   use interface_hydrostatic_eq_WL_peos
00027   use interface_calc_surface
00028   use interface_update_grfield
00029   use interface_update_matter
00030   use interface_update_surface
00031   use interface_error_matter
00032   use interface_source_HaC_WL_EMF
00033   use interface_source_MoC_WL_EMF
00034   use interface_source_trfreeG_WL_EMF
00035   use interface_source_trG_WL_EMF
00036   use interface_source_HaC_WL
00037   use interface_source_MoC_WL
00038   use interface_source_trfreeG_WL
00039   use interface_source_trG_WL
00040   use interface_source_MWtemp_WL
00041   use interface_source_MWspatial_WL
00042   use interface_compute_fnc_division
00043   use interface_compute_fnc_multiple
00044   use interface_error_metric_type0
00045   use interface_error_metric_type2_axisym
00046   use interface_calc_integrability_modify_At
00047   use interface_calc_integrability_modify_Aphi
00048   implicit none
00049   real(long), pointer :: sou(:,:,:), pot(:,:,:)
00050   real(long), pointer :: sou1(:,:,:), sou2(:,:,:), sou3(:,:,:)
00051   real(long), pointer :: sou4(:,:,:), sou5(:,:,:), sou6(:,:,:)
00052   real(long), pointer :: potf(:,:,:), emd_bak(:,:,:), pot_bak(:,:,:)
00053   real(long), pointer :: pottf(:,:,:), utf_bak(:,:,:)
00054   real(long), pointer :: potxf(:,:,:), uxf_bak(:,:,:)
00055   real(long), pointer :: potyf(:,:,:), uyf_bak(:,:,:)
00056   real(long), pointer :: potzf(:,:,:), uzf_bak(:,:,:)
00057   real(long), pointer :: potrs(:,:)
00058   real(long), pointer :: potx(:,:,:), poty(:,:,:), potz(:,:,:)
00059   real(long), pointer :: potx_bak(:,:,:), poty_bak(:,:,:), potz_bak(:,:,:)
00060   real(long), pointer :: souvec(:,:,:,:), souten(:,:,:,:)
00061   real(long), pointer :: pot1(:,:,:), pot2(:,:,:), pot3(:,:,:)
00062   real(long), pointer :: pot4(:,:,:), pot5(:,:,:), pot6(:,:,:)
00063   real(long), pointer :: bk1(:,:,:), bk2(:,:,:), bk3(:,:,:)
00064   real(long), pointer :: bk4(:,:,:), bk5(:,:,:), bk6(:,:,:)
00065   real(long), pointer :: work(:,:,:)
00066 !  real(long) :: work(0:nnrg,0:nntg,0:nnpg)
00067   real(long) :: error_emd, count, error_field
00068   integer    :: iter_count, flag = 0, hydro_iter = 4, flag_charge = 0
00069   integer    :: irf, itf, ipf, irg, itg, ipg, ihy
00070 !
00071   real(long) :: error_all(20)
00072   integer    :: flag_all(20) = (/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
00073   real(long) :: tic1, tic2
00074 !
00075   call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00076   call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00077   call alloc_array3d(potf,0,nrf,0,ntf,0,npf)
00078   call alloc_array3d(pottf,0,nrf,0,ntf,0,npf)
00079   call alloc_array3d(potxf,0,nrf,0,ntf,0,npf)
00080   call alloc_array3d(potyf,0,nrf,0,ntf,0,npf)
00081   call alloc_array3d(potzf,0,nrf,0,ntf,0,npf)
00082   call alloc_array3d(emd_bak,0,nrf,0,ntf,0,npf)
00083   call alloc_array3d(utf_bak,0,nrf,0,ntf,0,npf)
00084   call alloc_array3d(uxf_bak,0,nrf,0,ntf,0,npf)
00085   call alloc_array3d(uyf_bak,0,nrf,0,ntf,0,npf)
00086   call alloc_array3d(uzf_bak,0,nrf,0,ntf,0,npf)
00087   call alloc_array2d(potrs,0,ntf,0,npf)
00088   call alloc_array3d(potx,0,nrg,0,ntg,0,npg)
00089   call alloc_array3d(poty,0,nrg,0,ntg,0,npg)
00090   call alloc_array3d(potz,0,nrg,0,ntg,0,npg)
00091   call alloc_array4d(souvec,0,nrg,0,ntg,0,npg,1,3)
00092   call alloc_array4d(souten,0,nrg,0,ntg,0,npg,1,6)
00093   call alloc_array3d(sou1,0,nrg,0,ntg,0,npg)
00094   call alloc_array3d(sou2,0,nrg,0,ntg,0,npg)
00095   call alloc_array3d(sou3,0,nrg,0,ntg,0,npg)
00096   call alloc_array3d(sou4,0,nrg,0,ntg,0,npg)
00097   call alloc_array3d(sou5,0,nrg,0,ntg,0,npg)
00098   call alloc_array3d(sou6,0,nrg,0,ntg,0,npg)
00099   call alloc_array3d(pot1,0,nrg,0,ntg,0,npg)
00100   call alloc_array3d(pot2,0,nrg,0,ntg,0,npg)
00101   call alloc_array3d(pot3,0,nrg,0,ntg,0,npg)
00102   call alloc_array3d(pot4,0,nrg,0,ntg,0,npg)
00103   call alloc_array3d(pot5,0,nrg,0,ntg,0,npg)
00104   call alloc_array3d(pot6,0,nrg,0,ntg,0,npg)
00105   call alloc_array3d(pot_bak,0,nrg,0,ntg,0,npg)
00106   call alloc_array3d(potx_bak,0,nrg,0,ntg,0,npg)
00107   call alloc_array3d(poty_bak,0,nrg,0,ntg,0,npg)
00108   call alloc_array3d(potz_bak,0,nrg,0,ntg,0,npg)
00109   call alloc_array3d(bk1,0,nrg,0,ntg,0,npg)
00110   call alloc_array3d(bk2,0,nrg,0,ntg,0,npg)
00111   call alloc_array3d(bk3,0,nrg,0,ntg,0,npg)
00112   call alloc_array3d(bk4,0,nrg,0,ntg,0,npg)
00113   call alloc_array3d(bk5,0,nrg,0,ntg,0,npg)
00114   call alloc_array3d(bk6,0,nrg,0,ntg,0,npg)
00115   call alloc_array3d(work,0,nrg,0,ntg,0,npg)
00116 !
00117   iter_count = 0
00118   do
00119     iter_count = iter_count + 1      
00120     count = dble(iter_count) 
00121     conv_gra = dmin1(conv0_gra,conv_ini*count)
00122     conv_den = dmin1(conv0_den,conv_ini*count)
00123 !
00124 call cpu_time(tic1)
00125     emdg = 0.0d0
00126     call calc_vector_x_grav(1)
00127     call calc_vector_x_matter(1)
00128     call calc_vector_phi_grav(1)
00129     call calc_vector_phi_matter(1)
00130     call interpo_fl2gr(emd, emdg)
00131     call calc_shift_down2up
00132     call calc_shift2rotshift
00133 !
00134 !==============================================
00135     if (chgra == 'i') call cleargeometry
00136     if (chgra /= 'i') then 
00137       call cristoffel_midpoint
00138       call cristoffel_gridpoint
00139 !call cpu_time(tic1)
00140       call riccitensor_midpoint
00141 !call cpu_time(tic2)
00142 !write(6,*)'time riccitensor_midpoint', tic2-tic1 ; tic1=tic2
00143     end if
00144     if (chgra == 'h'.or.chgra == 'c'.or.chgra == 'C' &
00145         &.or.chgra == 'H'.or.chgra == 'W') then 
00146       call liegmab_midpoint
00147       call liegmab_gridpoint
00148     end if
00149 !
00150     call excurve_WL_midpoint
00151     call excurve_WL_gridpoint
00152     call compute_fnc_multiple(alph,psi,alps)
00153     call compute_fnc_multiple(alph,va,alva)
00154     call derivatives_EM1form
00155     call faraday
00156     call faraday_gridpoint
00157     call faraday_raise_WL
00158     call faraday_raise_gridpoint_WL
00159     call SEM_tensor_EMF
00160     call SEM_and_current_on_SFC
00161 call cpu_time(tic2)
00162 write(6,*)'time SEM_and_current_on_SFC', tic2-tic1 ; tic1=tic2
00163 !
00164 !
00165 !==============================================
00166 !$omp parallel sections num_threads(2)
00167 ! -- psi
00168 !$omp section
00169     call source_HaC_WL_EMF(sou1)
00170 !WL    call source_HaC_WL(sou1)
00171     call poisson_solver(sou1,pot1)
00172     pot1(0:nrg,0:ntg,0:npg) = pot1(0:nrg,0:ntg,0:npg) + 1.0d0
00173 ! -- alps
00174 !$omp section
00175     call source_trG_WL_EMF(sou2)
00176 !WL    call source_trG_WL(sou2)
00177     call poisson_solver(sou2,pot2)
00178     pot2(0:nrg,0:ntg,0:npg) = pot2(0:nrg,0:ntg,0:npg) + 1.0d0
00179 !$omp end parallel sections
00180 ! -- psi update
00181     pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00182     call update_grfield(pot1,psi,conv_gra)
00183 !WIND    call update_parameter_WL_MHD(conv_den)
00184     call update_parameter_WL_MHD_GS(conv_den)
00185 !WL    call update_parameter_axisym_WL_peos(conv_den)
00186 !    call error_metric_type0(psi,pot_bak,error_all(1),flag_all(1),'ns')
00187     call error_metric_type2_axisym(psi,pot_bak,error_all(1),flag_all(1),'ns')
00188 ! -- alph update
00189     pot_bak(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)
00190     call compute_alps2alph(pot2,psi)
00191     call update_grfield(pot2,alph,conv_gra)
00192 !WIND    call update_parameter_WL_MHD(conv_den)
00193     call update_parameter_WL_MHD_GS(conv_den)
00194 !WL     call update_parameter_axisym_WL_peos(conv_den)
00195 !    call error_metric_type0(alph,pot_bak,error_all(2),flag_all(2),'ns')
00196     call error_metric_type2_axisym(alph,pot_bak,error_all(2),flag_all(2),'ns')
00197 ! --
00198 !
00199     call printout_error_metric_combined_2(iter_count,error_all(1),error_all(2))
00200 !
00201 call cpu_time(tic2)
00202 write(6,*)'time psi alpha', tic2-tic1 ; tic1=tic2
00203 ! -- shift
00204     call source_MoC_WL_EMF(souvec)
00205 !WL    call source_MoC_WL(souvec)
00206 !$omp parallel sections num_threads(3)
00207 !$omp section
00208     sou1(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 1)
00209     call poisson_solver(sou1,potx)
00210 !$omp section
00211     sou2(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 2)
00212     call poisson_solver(sou2,poty)
00213 !$omp section
00214     sou3(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 3)
00215     call poisson_solver(sou3,potz)
00216 !$omp end parallel sections
00217 !
00218 !!    call poisson_solver(sou,pot)
00219 !!    work(  0:nrg, 0:ntg, 0:npg)    = potx(0:nrg, 0:ntg, 0:npg)
00220 !!    souvec(0:nrg, 0:ntg, 0:npg, 1) = work(0:nrg, 0:ntg, 0:npg)
00221 !!    work(  0:nrg, 0:ntg, 0:npg)    = poty(0:nrg, 0:ntg, 0:npg)
00222 !!    souvec(0:nrg, 0:ntg, 0:npg, 2) = work(0:nrg, 0:ntg, 0:npg)
00223 !!    work(  0:nrg, 0:ntg, 0:npg)    = potz(0:nrg, 0:ntg, 0:npg)
00224 !!    souvec(0:nrg, 0:ntg, 0:npg, 3) = work(0:nrg, 0:ntg, 0:npg)
00225 !!    call compute_shift(potx, poty, potz, souvec, pot)
00226 !
00227 !test gauge
00228 !
00229     potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00230     poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00231     potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00232     bvxd(0:nrg,0:ntg,0:npg) = potx(0:nrg,0:ntg,0:npg)
00233     bvyd(0:nrg,0:ntg,0:npg) = poty(0:nrg,0:ntg,0:npg)
00234     bvzd(0:nrg,0:ntg,0:npg) = potz(0:nrg,0:ntg,0:npg)
00235 !! CHECK
00236 !    call calc_shift_down2up
00237 !    call gauge_shift
00238 !    call calc_shift_up2down
00239 !! CHECK
00240 !   improve by iteration
00241     potx(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00242     poty(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00243     potz(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00244     bvxd(0:nrg,0:ntg,0:npg) = potx_bak(0:nrg,0:ntg,0:npg)
00245     bvyd(0:nrg,0:ntg,0:npg) = poty_bak(0:nrg,0:ntg,0:npg)
00246     bvzd(0:nrg,0:ntg,0:npg) = potz_bak(0:nrg,0:ntg,0:npg)
00247 !test gauge end
00248 !
00249 !gauge    potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00250 !gauge    poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00251 !gauge    potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00252     call update_grfield(potx,bvxd,conv_gra)
00253     call update_grfield(poty,bvyd,conv_gra)
00254     call update_grfield(potz,bvzd,conv_gra)
00255 !WIND    call update_parameter_WL_MHD(conv_den)
00256     call update_parameter_WL_MHD_GS(conv_den)
00257 !WL    call update_parameter_axisym_WL_peos(conv_den)
00258     call error_metric_type2_axisym(bvxd,potx_bak,error_all(3),flag_all(3),'ns')
00259     call error_metric_type2_axisym(bvyd,poty_bak,error_all(4),flag_all(4),'ns')
00260     call error_metric_type2_axisym(bvzd,potz_bak,error_all(5),flag_all(5),'ns')
00261     call printout_error_metric_combined_3(iter_count,error_all(3), &
00262     &                                   error_all(4),error_all(5))
00263 !
00264 call cpu_time(tic2)
00265 write(6,*)'time shift', tic2-tic1 ; tic1=tic2
00266 !==============================================
00267 !
00268 ! --  For hij
00269 !
00270 !! bvzd=0.0
00271 !! hxxd=0.0;hxyd=0.0;hxzd=0.0
00272 !! hyyd=0.0;hyzd=0.0;hzzd=0.0
00273 !! va=0.0;vaxd=0.0;vayd=0.0;vazd=0.0
00274 !! go to 9001
00275 ! go to 9000
00276 !testtest  if (chgra == 'i') go to 9000
00277 !
00278     call source_trfreeG_WL_EMF(souten)
00279 !WL    call source_trfreeG_WL(souten)
00280 !
00281     pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00282     bk1(0:nrg,0:ntg,0:npg) = hxxd(0:nrg,0:ntg,0:npg)
00283     bk2(0:nrg,0:ntg,0:npg) = hxyd(0:nrg,0:ntg,0:npg)
00284     bk3(0:nrg,0:ntg,0:npg) = hxzd(0:nrg,0:ntg,0:npg)
00285     bk4(0:nrg,0:ntg,0:npg) = hyyd(0:nrg,0:ntg,0:npg)
00286     bk5(0:nrg,0:ntg,0:npg) = hyzd(0:nrg,0:ntg,0:npg)
00287     bk6(0:nrg,0:ntg,0:npg) = hzzd(0:nrg,0:ntg,0:npg)
00288     hxxd(0:nrg,0:ntg,0:npg) = 0.0d0
00289     hxyd(0:nrg,0:ntg,0:npg) = 0.0d0
00290     hxzd(0:nrg,0:ntg,0:npg) = 0.0d0
00291     hyyd(0:nrg,0:ntg,0:npg) = 0.0d0
00292     hyzd(0:nrg,0:ntg,0:npg) = 0.0d0
00293     hzzd(0:nrg,0:ntg,0:npg) = 0.0d0
00294 !
00295 ! --  Poisson solver.  
00296 !
00297     if (chope == 'L') then
00298 !
00299       sou1(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,1)
00300       sou2(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,2)
00301       sou3(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,3)
00302       sou4(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,4)
00303       sou5(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,5)
00304       sou6(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,6)
00305 !$omp parallel sections num_threads(6)
00306 !$omp section
00307       call poisson_solver(sou1,pot1)
00308 !$omp section
00309       call poisson_solver(sou2,pot2)
00310 !$omp section
00311       call poisson_solver(sou3,pot3)
00312 !$omp section
00313       call poisson_solver(sou4,pot4)
00314 !$omp section
00315       call poisson_solver(sou5,pot5)
00316 !$omp section
00317       call poisson_solver(sou6,pot6)
00318 !$omp end parallel sections
00319 !
00320     end if
00321 !
00322 ! --  Helmholtz solver.  
00323 !
00324 !  if (chope == 'H') then
00325 !    if (ii == 1) then
00326 !      call hegr4(sou,pot,cosmpg,cosmpg,pnag,0,0,0,4)
00327 !      call hegr4sfout(ome,ssou,dssou,spot,cosmpg,cosmpg,pnag,0,0,0,4)
00328 !    end if
00329 !  end if
00330 !
00331     hxxd(0:nrg,0:ntg,0:npg) = pot1(0:nrg,0:ntg,0:npg)
00332     hxyd(0:nrg,0:ntg,0:npg) = pot2(0:nrg,0:ntg,0:npg)
00333     hxzd(0:nrg,0:ntg,0:npg) = pot3(0:nrg,0:ntg,0:npg)
00334     hyyd(0:nrg,0:ntg,0:npg) = pot4(0:nrg,0:ntg,0:npg)
00335     hyzd(0:nrg,0:ntg,0:npg) = pot5(0:nrg,0:ntg,0:npg)
00336     hzzd(0:nrg,0:ntg,0:npg) = pot6(0:nrg,0:ntg,0:npg)
00337 !
00338 ! --  Set det(gamma_{ab}) = psi^12 and Dirac gauge for new variable.
00339 !
00340 !    call gauge
00341     call invhij
00342     call gauge_hiju(conv_gra)
00343     call invhij_up2down
00344     call adjusthij
00345 !
00346 ! --  improve by iteration
00347 !
00348     pot(0:nrg,0:ntg,0:npg)  = psi(0:nrg,0:ntg,0:npg)
00349     pot1(0:nrg,0:ntg,0:npg) = hxxd(0:nrg,0:ntg,0:npg)
00350     pot2(0:nrg,0:ntg,0:npg) = hxyd(0:nrg,0:ntg,0:npg)
00351     pot3(0:nrg,0:ntg,0:npg) = hxzd(0:nrg,0:ntg,0:npg)
00352     pot4(0:nrg,0:ntg,0:npg) = hyyd(0:nrg,0:ntg,0:npg)
00353     pot5(0:nrg,0:ntg,0:npg) = hyzd(0:nrg,0:ntg,0:npg)
00354     pot6(0:nrg,0:ntg,0:npg) = hzzd(0:nrg,0:ntg,0:npg)
00355 !
00356     psi(0:nrg,0:ntg,0:npg)  = pot_bak(0:nrg,0:ntg,0:npg)
00357     hxxd(0:nrg,0:ntg,0:npg) = bk1(0:nrg,0:ntg,0:npg)
00358     hxyd(0:nrg,0:ntg,0:npg) = bk2(0:nrg,0:ntg,0:npg)
00359     hxzd(0:nrg,0:ntg,0:npg) = bk3(0:nrg,0:ntg,0:npg)
00360     hyyd(0:nrg,0:ntg,0:npg) = bk4(0:nrg,0:ntg,0:npg)
00361     hyzd(0:nrg,0:ntg,0:npg) = bk5(0:nrg,0:ntg,0:npg)
00362     hzzd(0:nrg,0:ntg,0:npg) = bk6(0:nrg,0:ntg,0:npg)
00363 !
00364     call update_grfield(pot_bak,psi,conv_gra)
00365 !###    call update_grfield(pot,psi,conv_gra)
00366     call update_grfield(pot1,hxxd,conv_gra)
00367     call update_grfield(pot2,hxyd,conv_gra)
00368     call update_grfield(pot3,hxzd,conv_gra)
00369     call update_grfield(pot4,hyyd,conv_gra)
00370     call update_grfield(pot5,hyzd,conv_gra)
00371     call update_grfield(pot6,hzzd,conv_gra)
00372 !
00373     call invhij
00374     call compute_fnc_multiple(alph,psi,alps)
00375 !
00376 !WIND    call update_parameter_WL_MHD(conv_den)
00377     call update_parameter_WL_MHD_GS(conv_den)
00378 !WL    call update_parameter_axisym_WL_peos(conv_den)
00379     call error_metric_type2_axisym(hxxd,bk1,error_all(6),flag_all(6),'ns')
00380     call error_metric_type2_axisym(hxyd,bk2,error_all(7),flag_all(7),'ns')
00381     call error_metric_type2_axisym(hxzd,bk3,error_all(8),flag_all(8),'ns')
00382     call error_metric_type2_axisym(hyyd,bk4,error_all(9),flag_all(9),'ns')
00383     call error_metric_type2_axisym(hyzd,bk5,error_all(10),flag_all(10),'ns')
00384     call error_metric_type2_axisym(hzzd,bk6,error_all(11),flag_all(11),'ns')
00385     call printout_error_metric_combined_3(iter_count,error_all(6), &
00386     &                                   error_all(7),error_all(8))
00387     call printout_error_metric_combined_3(iter_count,error_all(9), &
00388     &                                   error_all(10),error_all(11))
00389 !
00390 call cpu_time(tic2)
00391 write(6,*)'time hij', tic2-tic1 ; tic1=tic2
00392 ! --  IWM formalism
00393 !
00394  9000 continue
00395 !
00396 !GS
00397 !WL va=0.0;vaxd=0.0;vayd=0.0;vazd=0.0
00398 !WL go to 9001
00399 ! --  Maxwell equation component normal to hypersuface
00400 !
00401     call source_MWtemp_WL(sou)
00402 !    call poisson_solver(sou,pot)
00403     call poisson_solver_At(sou,pot,'axisym')
00404     call compute_fnc_division(pot,alph,work)
00405 !! Replace A_t inside of NS
00406 !!    call calc_integrability_modify_At(va,'ns')
00407 !!    call calc_integrability_modify_At(work,'ns')
00408 !! Replace A_t end
00409     pot_bak = va
00410     call update_grfield(work,va,conv_gra)
00411 !WIND    call update_parameter_WL_MHD(conv_den)
00412     call update_parameter_WL_MHD_GS(conv_den)
00413 !WL    call update_parameter_axisym_WL_peos(conv_den)
00414     call calc_integrability_modify_At(va,'ns')
00415     call compute_fnc_multiple(alph,va,alva)
00416     call error_metric_type2_axisym(va,pot_bak,error_all(12),flag_all(12),'ns')
00417     call printout_error_metric(iter_count,error_all(12))
00418 !
00419 call cpu_time(tic2)
00420 write(6,*)'time At', tic2-tic1 ; tic1=tic2
00421 !
00422 ! --  Maxwell equation components in hypersuface
00423 !
00424     call source_MWspatial_WL(souvec)
00425 !
00426 !$omp parallel sections num_threads(3)
00427 !$omp section
00428     sou1(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 1)
00429     call poisson_solver(sou1,potx)
00430 !$omp section
00431     sou2(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 2)
00432     call poisson_solver(sou2,poty)
00433 !$omp section
00434     sou3(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 3)
00435     call poisson_solver(sou3,potz)
00436 !!    call poisson_solver_Az(sou3,potz)
00437 !$omp end parallel sections
00438 !
00439     bk1(0:nrg,0:ntg,0:npg) = vaxd(0:nrg,0:ntg,0:npg)
00440     bk2(0:nrg,0:ntg,0:npg) = vayd(0:nrg,0:ntg,0:npg)
00441     bk3(0:nrg,0:ntg,0:npg) = vazd(0:nrg,0:ntg,0:npg)
00442     vaxd(0:nrg,0:ntg,0:npg) = potx(0:nrg,0:ntg,0:npg)
00443     vayd(0:nrg,0:ntg,0:npg) = poty(0:nrg,0:ntg,0:npg)
00444     vazd(0:nrg,0:ntg,0:npg) = potz(0:nrg,0:ntg,0:npg)
00445 !                                                                              
00446     call gauge_Coulomb
00447 !
00448 !   improve by iteration
00449     potx(0:nrg,0:ntg,0:npg) = vaxd(0:nrg,0:ntg,0:npg)
00450     poty(0:nrg,0:ntg,0:npg) = vayd(0:nrg,0:ntg,0:npg)
00451     potz(0:nrg,0:ntg,0:npg) = vazd(0:nrg,0:ntg,0:npg)
00452     vaxd(0:nrg,0:ntg,0:npg) = bk1(0:nrg,0:ntg,0:npg)
00453     vayd(0:nrg,0:ntg,0:npg) = bk2(0:nrg,0:ntg,0:npg)
00454     vazd(0:nrg,0:ntg,0:npg) = bk3(0:nrg,0:ntg,0:npg)
00455 !
00456     call update_grfield(potx,vaxd,conv_gra)
00457     call update_grfield(poty,vayd,conv_gra)
00458     call update_grfield(potz,vazd,conv_gra)
00459 !WIND    call update_parameter_WL_MHD(conv_den)
00460     call update_parameter_WL_MHD_GS(conv_den)
00461 !WL    call update_parameter_axisym_WL_peos(conv_den)
00462     call error_metric_type2_axisym(vaxd,bk1,error_all(13),flag_all(13),'ns')
00463     call error_metric_type2_axisym(vayd,bk2,error_all(14),flag_all(14),'ns')
00464     call error_metric_type2_axisym(vazd,bk3,error_all(15),flag_all(15),'ns')
00465     call printout_error_metric_combined_3(iter_count,error_all(13), &
00466     &                                  error_all(14),error_all(15))
00467 !
00468     call calc_Aphi_max
00469   if (mod(iter_count,30).eq.0) then
00470 !  if (mod(iter_count,20).eq.0) then
00471 !  if (mod(iter_count,10).eq.0) then
00472 !  if (mod(iter_count,4).eq.0) then
00473     call adjust_charge(flag_charge,'asympt')
00474 !    call adjust_charge_bisection(flag_charge,'asympt')
00475   end if
00476 call cpu_time(tic2)
00477 write(6,*)'time Ai', tic2-tic1 ; tic1=tic2
00478 !
00479 !!    call calc_MHDpar_charge
00480 ! Replace A_t inside of NS
00481 !call cpu_time(tic1)
00482 !    call calc_integrability_modify_Aphi(vaxd,vayd,'ns')
00483 !!!
00484 !!!    impose At(Aphi) function 
00485     call calc_integrability_modify_At(va,'ns')
00486 !!!
00487 !call cpu_time(tic2)
00488 !write(6,*)'time calc_integrability_modify_At', tic2-tic1 ; tic1=tic2
00489 ! Replace A_t end
00490 !
00491 !WL 9001 continue
00492 !!      itg = ntgeq; ipg = npgxzp
00493       itg = ntgeq-2; ipg = 2
00494 !!      itg = 73; ipg = 16
00495       open(15,file='test_vec',status='unknown')
00496         do irg = 0, nrg
00497           write(15,'(1p,9e20.12)')  rg(irg),  va(irg,itg,ipg)  &
00498               &                                , vaxd(irg,itg,ipg) &
00499               &                                , vayd(irg,itg,ipg) &
00500               &                                , vazd(irg,itg,ipg) &
00501               &                                , jtu(irg,itg,ipg) &
00502               &                                , jxu(irg,itg,ipg) &
00503               &                                , jyu(irg,itg,ipg) &
00504               &                                , jzu(irg,itg,ipg)
00505         end do
00506       close(15)
00507 !
00508       open(16,file='test_vec2',status='unknown')
00509         do irg = 0, nrg
00510           write(16,'(1p,12e20.12)') rg(irg),  psi(irg,itg,ipg)  &
00511               &                            , alph(irg,itg,ipg) &
00512               &                            , bvxd(irg,itg,ipg) &
00513               &                            , bvyd(irg,itg,ipg) &
00514               &                            , bvzd(irg,itg,ipg) &
00515               &                            , hxxd(irg,itg,ipg) &
00516               &                            , hxyd(irg,itg,ipg) &
00517               &                            , hxzd(irg,itg,ipg) &
00518               &                            , hyyd(irg,itg,ipg) &
00519               &                            , hyzd(irg,itg,ipg) &
00520               &                            , hzzd(irg,itg,ipg)
00521         end do
00522       close(16)
00523 !
00524 !      itg = ntgeq-2; ipg = 2
00525 !!      itg = 73; ipg = 16
00526       open(17,file='test_vec_3',status='unknown')
00527         do irg = 0, nrg
00528           write(17,'(1p,9e20.12)')  rg(irg), tfkij(irg,itg,ipg,1,1) &
00529               &                            , tfkij(irg,itg,ipg,1,2) &
00530               &                            , tfkij(irg,itg,ipg,1,3) &
00531               &                            , tfkij(irg,itg,ipg,2,2) &
00532               &                            , tfkij(irg,itg,ipg,2,3) &
00533               &                            , tfkij(irg,itg,ipg,3,3) &
00534               &                            , tfkijkij(irg,itg,ipg)  &
00535               &                            , trk(irg,itg,ipg)
00536         end do
00537       close(17)
00538 !==============================================
00539 !
00540 ! -- MHD equations.
00541 !testtest
00542 ! skip hydro equations
00543 !goto 9400
00544 !testtest
00545 do ihy = 1, hydro_iter
00546   if (ihy.eq.1) then 
00547     emd_bak = emd
00548     utf_bak = utf ; uxf_bak = uxf ; uyf_bak = uyf ; uzf_bak = uzf
00549   end if
00550     potf = emd 
00551     pottf = utf ; potxf = uxf ; potyf = uyf ; potzf = uzf
00552     call interpo_gr2fl_metric_WL
00553     call hydrostatic_eq_WL_MHD(potf,pottf,potxf,potyf,potzf)
00554 !WL    call hydrostatic_eq_WL_peos(potf)
00555     call calc_surface(potrs, potf)
00556     call update_matter(potf,emd,conv_den)
00557     call update_matter(pottf,utf,conv_den)
00558     call update_matter(potxf,uxf,conv_den)
00559     call update_matter(potyf,uyf,conv_den)
00560     call update_matter(potzf,uzf,conv_den)
00561     call update_surface(potrs,rs,conv_den)
00562     call reset_fluid
00563 !WIND    call update_parameter_WL_MHD(conv_den)
00564     call update_parameter_WL_MHD_GS(conv_den)
00565 !WL    call update_parameter_axisym_WL_peos(conv_den)
00566     call calc_vector_x_matter(1)
00567     call calc_vector_phi_matter(1)
00568 ! -- calculate error and print it out.
00569 if (ihy.eq.hydro_iter) then
00570     call error_matter(emd,emd_bak,error_all(16),flag_all(16))
00571     call error_matter(utf,utf_bak,error_all(17),flag_all(17))
00572     call error_matter(uxf,uxf_bak,error_all(18),flag_all(18))
00573     call error_matter(uyf,uyf_bak,error_all(19),flag_all(19))
00574     call error_matter(uzf,uzf_bak,error_all(20),flag_all(20))
00575 end if
00576 end do
00577 ! --
00578     call printout_error(iter_count,error_all(16))
00579     call printout_error(iter_count,error_all(17))
00580     call printout_error(iter_count,error_all(18))
00581     call printout_error(iter_count,error_all(19))
00582     call printout_error(iter_count,error_all(20))
00583 !
00584 call cpu_time(tic2)
00585 write(6,*)'time hydro', tic2-tic1 ; tic1=tic2
00586 !testtest
00587 !9400 continue
00588 !testtest
00589 !
00590 ! ---
00591 !
00592     if (flag_all(1)==0.and.flag_all(2)==0.and.flag_all(3)==0.and.&
00593     &   flag_all(4)==0.and.flag_all(5)==0.and.flag_all(6)==0.and.&
00594     &   flag_all(7)==0.and.flag_all(8)==0.and.flag_all(9)==0.and.&
00595     &   flag_all(10)==0.and.flag_all(11)==0.and.flag_all(12)==0.and.&
00596     &   flag_all(13)==0.and.flag_all(14)==0.and.flag_all(15)==0.and.&
00597     &   flag_all(16)==0.and.flag_all(17)==0.and.flag_all(18)==0.and.&
00598     &   flag_all(19)==0.and.flag_all(20)==0) exit
00599     if (iter_count >= iter_max) exit
00600 !
00601     flag = 0
00602     flag_all(1:20) = 0
00603 !
00604   end do
00605 !
00606   deallocate(sou)
00607   deallocate(pot)
00608   deallocate(potf)
00609   deallocate(pottf)
00610   deallocate(potxf)
00611   deallocate(potyf)
00612   deallocate(potzf)
00613   deallocate(emd_bak)
00614   deallocate(utf_bak)
00615   deallocate(uxf_bak)
00616   deallocate(uyf_bak)
00617   deallocate(uzf_bak)
00618   deallocate(potrs)
00619   deallocate(potx)
00620   deallocate(poty)
00621   deallocate(potz)
00622   deallocate(souvec)
00623   deallocate(souten)
00624   deallocate(sou1)
00625   deallocate(sou2)
00626   deallocate(sou3)
00627   deallocate(sou4)
00628   deallocate(sou5)
00629   deallocate(sou6)
00630   deallocate(pot1)
00631   deallocate(pot2)
00632   deallocate(pot3)
00633   deallocate(pot4)
00634   deallocate(pot5)
00635   deallocate(pot6)
00636   deallocate(pot_bak)
00637   deallocate(potx_bak)
00638   deallocate(poty_bak)
00639   deallocate(potz_bak)
00640   deallocate(bk1)
00641   deallocate(bk2)
00642   deallocate(bk3)
00643   deallocate(bk4)
00644   deallocate(bk5)
00645   deallocate(bk6)
00646   deallocate(work)
00647 end subroutine iteration_WL_MHD
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/iteration_WL_MHD.f90
