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_calc_surface
00027 use interface_update_grfield
00028 use interface_update_matter
00029 use interface_update_surface
00030 use interface_error_matter
00031 use interface_source_HaC_WL_EMF
00032 use interface_source_MoC_WL_EMF
00033 use interface_source_trfreeG_WL_EMF
00034 use interface_source_trG_WL_EMF
00035 use interface_source_MWtemp_WL
00036 use interface_source_MWspatial_WL
00037 use interface_compute_fnc_division
00038 use interface_compute_fnc_multiple
00039 use interface_error_metric_type0
00040 use interface_error_metric_type2_axisym
00041 use interface_calc_integrability_modify_At
00042 use interface_calc_integrability_modify_Aphi
00043 implicit none
00044 real(long), pointer :: sou(:,:,:), pot(:,:,:)
00045 real(long), pointer :: sou1(:,:,:), sou2(:,:,:), sou3(:,:,:)
00046 real(long), pointer :: sou4(:,:,:), sou5(:,:,:), sou6(:,:,:)
00047 real(long), pointer :: potf(:,:,:), emd_bak(:,:,:), pot_bak(:,:,:)
00048 real(long), pointer :: pottf(:,:,:), utf_bak(:,:,:)
00049 real(long), pointer :: potxf(:,:,:), uxf_bak(:,:,:)
00050 real(long), pointer :: potyf(:,:,:), uyf_bak(:,:,:)
00051 real(long), pointer :: potzf(:,:,:), uzf_bak(:,:,:)
00052 real(long), pointer :: potrs(:,:)
00053 real(long), pointer :: potx(:,:,:), poty(:,:,:), potz(:,:,:)
00054 real(long), pointer :: potx_bak(:,:,:), poty_bak(:,:,:), potz_bak(:,:,:)
00055 real(long), pointer :: souvec(:,:,:,:), souten(:,:,:,:)
00056 real(long), pointer :: pot1(:,:,:), pot2(:,:,:), pot3(:,:,:)
00057 real(long), pointer :: pot4(:,:,:), pot5(:,:,:), pot6(:,:,:)
00058 real(long), pointer :: bk1(:,:,:), bk2(:,:,:), bk3(:,:,:)
00059 real(long), pointer :: bk4(:,:,:), bk5(:,:,:), bk6(:,:,:)
00060 real(long), pointer :: work(:,:,:)
00061
00062 real(long) :: error_emd, count, error_field
00063 integer :: iter_count, flag = 0, hydro_iter = 4, flag_charge = 0
00064 integer :: irf, itf, ipf, irg, itg, ipg, ihy
00065
00066 real(long) :: error_all(20)
00067 integer :: flag_all(20) = (/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
00068 real(long) :: tic1, tic2
00069
00070 call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00071 call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00072 call alloc_array3d(potf,0,nrf,0,ntf,0,npf)
00073 call alloc_array3d(pottf,0,nrf,0,ntf,0,npf)
00074 call alloc_array3d(potxf,0,nrf,0,ntf,0,npf)
00075 call alloc_array3d(potyf,0,nrf,0,ntf,0,npf)
00076 call alloc_array3d(potzf,0,nrf,0,ntf,0,npf)
00077 call alloc_array3d(emd_bak,0,nrf,0,ntf,0,npf)
00078 call alloc_array3d(utf_bak,0,nrf,0,ntf,0,npf)
00079 call alloc_array3d(uxf_bak,0,nrf,0,ntf,0,npf)
00080 call alloc_array3d(uyf_bak,0,nrf,0,ntf,0,npf)
00081 call alloc_array3d(uzf_bak,0,nrf,0,ntf,0,npf)
00082 call alloc_array2d(potrs,0,ntf,0,npf)
00083 call alloc_array3d(potx,0,nrg,0,ntg,0,npg)
00084 call alloc_array3d(poty,0,nrg,0,ntg,0,npg)
00085 call alloc_array3d(potz,0,nrg,0,ntg,0,npg)
00086 call alloc_array4d(souvec,0,nrg,0,ntg,0,npg,1,3)
00087 call alloc_array4d(souten,0,nrg,0,ntg,0,npg,1,6)
00088 call alloc_array3d(sou1,0,nrg,0,ntg,0,npg)
00089 call alloc_array3d(sou2,0,nrg,0,ntg,0,npg)
00090 call alloc_array3d(sou3,0,nrg,0,ntg,0,npg)
00091 call alloc_array3d(sou4,0,nrg,0,ntg,0,npg)
00092 call alloc_array3d(sou5,0,nrg,0,ntg,0,npg)
00093 call alloc_array3d(sou6,0,nrg,0,ntg,0,npg)
00094 call alloc_array3d(pot1,0,nrg,0,ntg,0,npg)
00095 call alloc_array3d(pot2,0,nrg,0,ntg,0,npg)
00096 call alloc_array3d(pot3,0,nrg,0,ntg,0,npg)
00097 call alloc_array3d(pot4,0,nrg,0,ntg,0,npg)
00098 call alloc_array3d(pot5,0,nrg,0,ntg,0,npg)
00099 call alloc_array3d(pot6,0,nrg,0,ntg,0,npg)
00100 call alloc_array3d(pot_bak,0,nrg,0,ntg,0,npg)
00101 call alloc_array3d(potx_bak,0,nrg,0,ntg,0,npg)
00102 call alloc_array3d(poty_bak,0,nrg,0,ntg,0,npg)
00103 call alloc_array3d(potz_bak,0,nrg,0,ntg,0,npg)
00104 call alloc_array3d(bk1,0,nrg,0,ntg,0,npg)
00105 call alloc_array3d(bk2,0,nrg,0,ntg,0,npg)
00106 call alloc_array3d(bk3,0,nrg,0,ntg,0,npg)
00107 call alloc_array3d(bk4,0,nrg,0,ntg,0,npg)
00108 call alloc_array3d(bk5,0,nrg,0,ntg,0,npg)
00109 call alloc_array3d(bk6,0,nrg,0,ntg,0,npg)
00110 call alloc_array3d(work,0,nrg,0,ntg,0,npg)
00111
00112 iter_count = 0
00113 do
00114 iter_count = iter_count + 1
00115 count = dble(iter_count)
00116 conv_gra = dmin1(conv0_gra,conv_ini*count)
00117 conv_den = dmin1(conv0_den,conv_ini*count)
00118
00119 call cpu_time(tic1)
00120 call calc_vector_x_grav(1)
00121 call calc_vector_x_matter(1)
00122 call calc_vector_phi_grav(1)
00123 call calc_vector_phi_matter(1)
00124 call calc_shift_down2up
00125 call calc_shift2rotshift
00126
00127
00128
00129
00130
00131 if (chgra == 'i') call cleargeometry
00132 if (chgra /= 'i') then
00133 call cristoffel_midpoint
00134 call cristoffel_gridpoint
00135
00136 call riccitensor_midpoint
00137
00138
00139 end if
00140 if (chgra == 'h'.or.chgra == 'c'.or.chgra == 'C' &
00141 &.or.chgra == 'H'.or.chgra == 'W') then
00142 call liegmab_midpoint
00143 call liegmab_gridpoint
00144 end if
00145
00146 call excurve_WL_midpoint
00147 call excurve_WL_gridpoint
00148 call compute_fnc_multiple(alph,psi,alps)
00149 call compute_fnc_multiple(alph,va,alva)
00150 call derivatives_EM1form
00151 call faraday
00152 call faraday_gridpoint
00153 call faraday_raise_WL
00154 call faraday_raise_gridpoint_WL
00155 call SEM_tensor_EMF
00156 call SEM_and_current_on_SFC
00157 call cpu_time(tic2)
00158 write(6,*)'time SEM_and_current_on_SFC', tic2-tic1 ; tic1=tic2
00159
00160
00161
00162
00163
00164
00165
00166 call source_HaC_WL_EMF(sou1)
00167 call poisson_solver(sou1,pot1)
00168 pot1 = pot1 + 1.0d0
00169
00170
00171
00172 call source_trG_WL_EMF(sou2)
00173 call poisson_solver(sou2,pot2)
00174 pot2 = pot2 + 1.0d0
00175
00176
00177
00178 pot_bak = psi
00179 call update_grfield(pot1,psi,conv_gra)
00180
00181 call update_parameter_WL_MHD_GS(conv_den)
00182
00183 call error_metric_type2_axisym(psi,pot_bak,error_all(1),flag_all(1),'ns')
00184
00185
00186
00187
00188
00189
00190
00191
00192 call compute_alps2alph(pot2,psi)
00193 pot_bak = alph
00194 call update_grfield(pot2,alph,conv_gra)
00195
00196 call update_parameter_WL_MHD_GS(conv_den)
00197
00198 call error_metric_type2_axisym(alph,pot_bak,error_all(2),flag_all(2),'ns')
00199
00200 call printout_error_metric_combined_2(iter_count,error_all(1),error_all(2))
00201
00202 call cpu_time(tic2)
00203 write(6,*)'time psi alpha', tic2-tic1 ; tic1=tic2
00204
00205
00206 call source_MoC_WL_EMF(souvec)
00207
00208
00209 sou1(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 1)
00210 call poisson_solver(sou1,potx)
00211
00212 sou2(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 2)
00213 call poisson_solver(sou2,poty)
00214
00215 sou3(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 3)
00216 call poisson_solver(sou3,potz)
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00229 poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00230 potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00231 call update_grfield(potx,bvxd,conv_gra)
00232 call update_grfield(poty,bvyd,conv_gra)
00233 call update_grfield(potz,bvzd,conv_gra)
00234
00235 call update_parameter_WL_MHD_GS(conv_den)
00236 call error_metric_type2_axisym(bvxd,potx_bak,error_all(3),flag_all(3),'ns')
00237 call error_metric_type2_axisym(bvyd,poty_bak,error_all(4),flag_all(4),'ns')
00238 call error_metric_type2_axisym(bvzd,potz_bak,error_all(5),flag_all(5),'ns')
00239 call printout_error_metric_combined_3(iter_count,error_all(3), &
00240 & error_all(4),error_all(5))
00241
00242 call cpu_time(tic2)
00243 write(6,*)'time shift', tic2-tic1 ; tic1=tic2
00244
00245
00246
00247
00248
00249
00250
00251 call source_trfreeG_WL_EMF(souten)
00252
00253 pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00254 bk1(0:nrg,0:ntg,0:npg) = hxxd(0:nrg,0:ntg,0:npg)
00255 bk2(0:nrg,0:ntg,0:npg) = hxyd(0:nrg,0:ntg,0:npg)
00256 bk3(0:nrg,0:ntg,0:npg) = hxzd(0:nrg,0:ntg,0:npg)
00257 bk4(0:nrg,0:ntg,0:npg) = hyyd(0:nrg,0:ntg,0:npg)
00258 bk5(0:nrg,0:ntg,0:npg) = hyzd(0:nrg,0:ntg,0:npg)
00259 bk6(0:nrg,0:ntg,0:npg) = hzzd(0:nrg,0:ntg,0:npg)
00260 hxxd(0:nrg,0:ntg,0:npg) = 0.0d0
00261 hxyd(0:nrg,0:ntg,0:npg) = 0.0d0
00262 hxzd(0:nrg,0:ntg,0:npg) = 0.0d0
00263 hyyd(0:nrg,0:ntg,0:npg) = 0.0d0
00264 hyzd(0:nrg,0:ntg,0:npg) = 0.0d0
00265 hzzd(0:nrg,0:ntg,0:npg) = 0.0d0
00266
00267
00268
00269 if (chope == 'L') then
00270
00271 sou1(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,1)
00272 sou2(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,2)
00273 sou3(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,3)
00274 sou4(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,4)
00275 sou5(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,5)
00276 sou6(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,6)
00277
00278
00279 call poisson_solver(sou1,pot1)
00280
00281 call poisson_solver(sou2,pot2)
00282
00283 call poisson_solver(sou3,pot3)
00284
00285 call poisson_solver(sou4,pot4)
00286
00287 call poisson_solver(sou5,pot5)
00288
00289 call poisson_solver(sou6,pot6)
00290
00291
00292 end if
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 hxxd(0:nrg,0:ntg,0:npg) = pot1(0:nrg,0:ntg,0:npg)
00304 hxyd(0:nrg,0:ntg,0:npg) = pot2(0:nrg,0:ntg,0:npg)
00305 hxzd(0:nrg,0:ntg,0:npg) = pot3(0:nrg,0:ntg,0:npg)
00306 hyyd(0:nrg,0:ntg,0:npg) = pot4(0:nrg,0:ntg,0:npg)
00307 hyzd(0:nrg,0:ntg,0:npg) = pot5(0:nrg,0:ntg,0:npg)
00308 hzzd(0:nrg,0:ntg,0:npg) = pot6(0:nrg,0:ntg,0:npg)
00309
00310
00311
00312
00313 call invhij
00314 call gauge_hiju(conv_gra)
00315 call invhij_up2down
00316 call adjusthij
00317
00318
00319
00320 pot(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00321 pot1(0:nrg,0:ntg,0:npg) = hxxd(0:nrg,0:ntg,0:npg)
00322 pot2(0:nrg,0:ntg,0:npg) = hxyd(0:nrg,0:ntg,0:npg)
00323 pot3(0:nrg,0:ntg,0:npg) = hxzd(0:nrg,0:ntg,0:npg)
00324 pot4(0:nrg,0:ntg,0:npg) = hyyd(0:nrg,0:ntg,0:npg)
00325 pot5(0:nrg,0:ntg,0:npg) = hyzd(0:nrg,0:ntg,0:npg)
00326 pot6(0:nrg,0:ntg,0:npg) = hzzd(0:nrg,0:ntg,0:npg)
00327
00328 psi(0:nrg,0:ntg,0:npg) = pot_bak(0:nrg,0:ntg,0:npg)
00329 hxxd(0:nrg,0:ntg,0:npg) = bk1(0:nrg,0:ntg,0:npg)
00330 hxyd(0:nrg,0:ntg,0:npg) = bk2(0:nrg,0:ntg,0:npg)
00331 hxzd(0:nrg,0:ntg,0:npg) = bk3(0:nrg,0:ntg,0:npg)
00332 hyyd(0:nrg,0:ntg,0:npg) = bk4(0:nrg,0:ntg,0:npg)
00333 hyzd(0:nrg,0:ntg,0:npg) = bk5(0:nrg,0:ntg,0:npg)
00334 hzzd(0:nrg,0:ntg,0:npg) = bk6(0:nrg,0:ntg,0:npg)
00335
00336 call update_grfield(pot_bak,psi,conv_gra)
00337 call update_grfield(pot1,hxxd,conv_gra)
00338 call update_grfield(pot2,hxyd,conv_gra)
00339 call update_grfield(pot3,hxzd,conv_gra)
00340 call update_grfield(pot4,hyyd,conv_gra)
00341 call update_grfield(pot5,hyzd,conv_gra)
00342 call update_grfield(pot6,hzzd,conv_gra)
00343
00344 call invhij
00345 call compute_fnc_multiple(alph,psi,alps)
00346
00347
00348 call update_parameter_WL_MHD_GS(conv_den)
00349 call error_metric_type2_axisym(hxxd,bk1,error_all(6),flag_all(6),'ns')
00350 call error_metric_type2_axisym(hxyd,bk2,error_all(7),flag_all(7),'ns')
00351 call error_metric_type2_axisym(hxzd,bk3,error_all(8),flag_all(8),'ns')
00352 call error_metric_type2_axisym(hyyd,bk4,error_all(9),flag_all(9),'ns')
00353 call error_metric_type2_axisym(hyzd,bk5,error_all(10),flag_all(10),'ns')
00354 call error_metric_type2_axisym(hzzd,bk6,error_all(11),flag_all(11),'ns')
00355 call printout_error_metric_combined_3(iter_count,error_all(6), &
00356 & error_all(7),error_all(8))
00357 call printout_error_metric_combined_3(iter_count,error_all(9), &
00358 & error_all(10),error_all(11))
00359
00360 call cpu_time(tic2)
00361 write(6,*)'time hij', tic2-tic1 ; tic1=tic2
00362
00363
00364 9000 continue
00365
00366
00367
00368
00369 call source_MWtemp_WL(sou)
00370
00371 call poisson_solver_At(sou,pot,'axisym')
00372 call compute_fnc_division(pot,alph,work)
00373
00374
00375
00376
00377 pot_bak = va
00378 call update_grfield(work,va,conv_gra)
00379
00380 call update_parameter_WL_MHD_GS(conv_den)
00381 call compute_fnc_multiple(alph,va,alva)
00382 call error_metric_type2_axisym(va,pot_bak,error_all(12),flag_all(12),'ns')
00383 call printout_error_metric(iter_count,error_all(12))
00384
00385 call cpu_time(tic2)
00386 write(6,*)'time At', tic2-tic1 ; tic1=tic2
00387
00388
00389
00390 call source_MWspatial_WL(souvec)
00391
00392
00393
00394 sou1(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 1)
00395 call poisson_solver(sou1,potx)
00396
00397 sou2(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 2)
00398 call poisson_solver(sou2,poty)
00399
00400 sou3(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 3)
00401 call poisson_solver(sou3,potz)
00402
00403
00404
00405 bk1(0:nrg,0:ntg,0:npg) = vaxd(0:nrg,0:ntg,0:npg)
00406 bk2(0:nrg,0:ntg,0:npg) = vayd(0:nrg,0:ntg,0:npg)
00407 bk3(0:nrg,0:ntg,0:npg) = vazd(0:nrg,0:ntg,0:npg)
00408 vaxd(0:nrg,0:ntg,0:npg) = potx(0:nrg,0:ntg,0:npg)
00409 vayd(0:nrg,0:ntg,0:npg) = poty(0:nrg,0:ntg,0:npg)
00410 vazd(0:nrg,0:ntg,0:npg) = potz(0:nrg,0:ntg,0:npg)
00411
00412 call gauge_Coulomb
00413
00414
00415 potx(0:nrg,0:ntg,0:npg) = vaxd(0:nrg,0:ntg,0:npg)
00416 poty(0:nrg,0:ntg,0:npg) = vayd(0:nrg,0:ntg,0:npg)
00417 potz(0:nrg,0:ntg,0:npg) = vazd(0:nrg,0:ntg,0:npg)
00418 vaxd(0:nrg,0:ntg,0:npg) = bk1(0:nrg,0:ntg,0:npg)
00419 vayd(0:nrg,0:ntg,0:npg) = bk2(0:nrg,0:ntg,0:npg)
00420 vazd(0:nrg,0:ntg,0:npg) = bk3(0:nrg,0:ntg,0:npg)
00421
00422 call update_grfield(potx,vaxd,conv_gra)
00423 call update_grfield(poty,vayd,conv_gra)
00424 call update_grfield(potz,vazd,conv_gra)
00425
00426 call update_parameter_WL_MHD_GS(conv_den)
00427 call error_metric_type2_axisym(vaxd,bk1,error_all(13),flag_all(13),'ns')
00428 call error_metric_type2_axisym(vayd,bk2,error_all(14),flag_all(14),'ns')
00429 call error_metric_type2_axisym(vazd,bk3,error_all(15),flag_all(15),'ns')
00430 call printout_error_metric_combined_3(iter_count,error_all(13), &
00431 & error_all(14),error_all(15))
00432
00433 call calc_Aphi_max
00434 if (mod(iter_count,30).eq.0) then
00435
00436
00437
00438 call adjust_charge(flag_charge,'asympt')
00439
00440 end if
00441 call cpu_time(tic2)
00442 write(6,*)'time Ai', tic2-tic1 ; tic1=tic2
00443
00444
00445
00446
00447
00448
00449
00450 call calc_integrability_modify_At(va,'ns')
00451
00452
00453
00454
00455
00456
00457 itg = ntgeq-2; ipg = 2
00458
00459 open(15,file='test_vec',status='unknown')
00460 do irg = 0, nrg
00461 write(15,'(1p,9e20.12)') rg(irg), va(irg,itg,ipg) &
00462 & , vaxd(irg,itg,ipg) &
00463 & , vayd(irg,itg,ipg) &
00464 & , vazd(irg,itg,ipg) &
00465 & , jtu(irg,itg,ipg) &
00466 & , jxu(irg,itg,ipg) &
00467 & , jyu(irg,itg,ipg) &
00468 & , jzu(irg,itg,ipg)
00469 end do
00470 close(15)
00471
00472 open(16,file='test_vec2',status='unknown')
00473 do irg = 0, nrg
00474 write(16,'(1p,12e20.12)') rg(irg), psi(irg,itg,ipg) &
00475 & , alph(irg,itg,ipg) &
00476 & , bvxd(irg,itg,ipg) &
00477 & , bvyd(irg,itg,ipg) &
00478 & , bvzd(irg,itg,ipg) &
00479 & , hxxd(irg,itg,ipg) &
00480 & , hxyd(irg,itg,ipg) &
00481 & , hxzd(irg,itg,ipg) &
00482 & , hyyd(irg,itg,ipg) &
00483 & , hyzd(irg,itg,ipg) &
00484 & , hzzd(irg,itg,ipg)
00485 end do
00486 close(16)
00487
00488
00489
00490
00491
00492
00493
00494 do ihy = 1, hydro_iter
00495 if (ihy.eq.1) then
00496 emd_bak = emd
00497 utf_bak = utf ; uxf_bak = uxf ; uyf_bak = uyf ; uzf_bak = uzf
00498 end if
00499 potf = emd
00500 pottf = utf ; potxf = uxf ; potyf = uyf ; potzf = uzf
00501
00502 call hydrostatic_eq_WL_MHD(potf,pottf,potxf,potyf,potzf)
00503 call calc_surface(potrs, potf)
00504 call update_matter(potf,emd,conv_den)
00505 call update_matter(pottf,utf,conv_den)
00506 call update_matter(potxf,uxf,conv_den)
00507 call update_matter(potyf,uyf,conv_den)
00508 call update_matter(potzf,uzf,conv_den)
00509 call update_surface(potrs,rs,conv_den)
00510 call reset_fluid
00511
00512 call update_parameter_WL_MHD_GS(conv_den)
00513 call calc_vector_x_matter(1)
00514 call calc_vector_phi_matter(1)
00515
00516 if (ihy.eq.hydro_iter) then
00517 call error_matter(emd,emd_bak,error_all(16),flag_all(16))
00518 call error_matter(utf,utf_bak,error_all(17),flag_all(17))
00519 call error_matter(uxf,uxf_bak,error_all(18),flag_all(18))
00520 call error_matter(uyf,uyf_bak,error_all(19),flag_all(19))
00521 call error_matter(uzf,uzf_bak,error_all(20),flag_all(20))
00522 end if
00523 end do
00524
00525 call printout_error(iter_count,error_all(16))
00526 call printout_error(iter_count,error_all(17))
00527 call printout_error(iter_count,error_all(18))
00528 call printout_error(iter_count,error_all(19))
00529 call printout_error(iter_count,error_all(20))
00530
00531 call cpu_time(tic2)
00532 write(6,*)'time hydro', tic2-tic1 ; tic1=tic2
00533
00534
00535
00536
00537
00538
00539 if (flag_all(1)==0.and.flag_all(2)==0.and.flag_all(3)==0.and.&
00540 & flag_all(4)==0.and.flag_all(5)==0.and.flag_all(6)==0.and.&
00541 & flag_all(7)==0.and.flag_all(8)==0.and.flag_all(9)==0.and.&
00542 & flag_all(10)==0.and.flag_all(11)==0.and.flag_all(12)==0.and.&
00543 & flag_all(13)==0.and.flag_all(14)==0.and.flag_all(15)==0.and.&
00544 & flag_all(16)==0.and.flag_all(17)==0.and.flag_all(18)==0.and.&
00545 & flag_all(19)==0.and.flag_all(20)==0) exit
00546 if (iter_count >= iter_max) exit
00547
00548 flag = 0
00549 flag_all(1:20) = 0
00550
00551 end do
00552
00553 deallocate(sou)
00554 deallocate(pot)
00555 deallocate(potf)
00556 deallocate(pottf)
00557 deallocate(potxf)
00558 deallocate(potyf)
00559 deallocate(potzf)
00560 deallocate(emd_bak)
00561 deallocate(utf_bak)
00562 deallocate(uxf_bak)
00563 deallocate(uyf_bak)
00564 deallocate(uzf_bak)
00565 deallocate(potrs)
00566 deallocate(potx)
00567 deallocate(poty)
00568 deallocate(potz)
00569 deallocate(souvec)
00570 deallocate(souten)
00571 deallocate(sou1)
00572 deallocate(sou2)
00573 deallocate(sou3)
00574 deallocate(sou4)
00575 deallocate(sou5)
00576 deallocate(sou6)
00577 deallocate(pot1)
00578 deallocate(pot2)
00579 deallocate(pot3)
00580 deallocate(pot4)
00581 deallocate(pot5)
00582 deallocate(pot6)
00583 deallocate(pot_bak)
00584 deallocate(potx_bak)
00585 deallocate(poty_bak)
00586 deallocate(potz_bak)
00587 deallocate(bk1)
00588 deallocate(bk2)
00589 deallocate(bk3)
00590 deallocate(bk4)
00591 deallocate(bk5)
00592 deallocate(bk6)
00593 deallocate(work)
00594 end subroutine iteration_WL_MHD
