00001
00002 include '../Include_file/include_modulefiles_plot_BNS_CF_3mpt.f90'
00003 include '../Include_file/include_PEOS_modulefile.f90'
00004 include '../Include_file/include_modulefiles_analysis_plot_BNS_CF_3mpt.f90'
00005 include '../Include_file/include_interface_modulefiles_plot_BNS_CF_3mpt.f90'
00006 include '../Include_file/include_interface_modulefiles_analysis_plot_BNS_CF_3mpt.f90'
00007 include '../Include_file/include_subroutines_plot_BNS_CF_3mpt.f90'
00008 include '../Include_file/include_subroutines_analysis_plot_BNS_CF_3mpt.f90'
00009 include '../Include_file/include_PEOS_subroutines.f90'
00010 include '../Include_file/include_functions.f90'
00011
00012
00013
00014
00015 PROGRAM coc2cac
00016
00017 use grid_parameter_binary_excision
00018 use phys_constant
00019 use grid_parameter
00020 use coordinate_grav_r
00021 use coordinate_grav_extended
00022 use interface_modules_cartesian
00023 use interface_calc_gradvep_export
00024 use trigonometry_grav_phi
00025 use def_binary_parameter, only : dis
00026 use def_matter_parameter_mpt
00027 use interface_interpo_lag4th_2Dsurf
00028 use interface_IO_input_CF_grav_export
00029 use interface_IO_input_CF_flco_export
00030 use interface_IO_input_CF_surf_export
00031 use interface_excurve_CF_gridpoint_export
00032 use interface_interpo_gr2fl_metric_CF_export
00033 implicit none
00034 integer :: impt, impt_ex, ico, irr, isp
00035 character(30) :: char1
00036 character*400 :: dir_path
00037 character(len=1) :: np(5) = (/'1', '2','3', '4', '5'/)
00038 real(8) :: rr3, rrcm, xc,yc,zc, xc_p1, yc_p1, zc_p1, xc_p2, yc_p2, zc_p2, dis_cm
00039 real(8) :: xc_p3, yc_p3, zc_p3
00040 real(8) :: xcac, ycac, zcac
00041 real(8) :: xcoc, ycoc, zcoc
00042 real(8) :: emdca, psica, alphca, bvxdca, bvydca, bvzdca, psi4ca, psif4ca
00043 real(8) :: hca, preca, rhoca, eneca, epsca
00044 real(8) :: vxu, vyu, vzu, lam_p1, lam_p2
00045 real(8) :: bxcor, bycor, bzcor, bvxdfca, bvydfca, bvzdfca, psifca, alphfca
00046 real(8) :: gxx, gxy, gxz, gyy, gyz, gzz, kxx, kxy, kxz, kyy, kyz, kzz
00047 real(8) :: axx, axy, axz, ayy, ayz, azz
00048 real(8) :: ome_p1, ber_p1, radi_p1, r_surf_p1
00049 real(8) :: ome_p2, ber_p2, radi_p2, r_surf_p2
00050 integer :: nrg_p1, ntg_p1, npg_p1, nrf_p1, ntf_p1, npf_p1
00051 integer :: nrg_p2, ntg_p2, npg_p2, nrf_p2, ntf_p2, npf_p2
00052 integer :: nrg_p3, ntg_p3, npg_p3, nrf_p3, ntf_p3, npf_p3
00053
00054 real(long) :: rc_p1, thc_p1, phic_p1, varpic_p1, rcf_p1, rsca_p1
00055 real(long) :: rc_p2, thc_p2, phic_p2, varpic_p2, rcf_p2, rsca_p2
00056 real(long) :: rc_p3, thc_p3, phic_p3, varpic_p3
00057 real(long) :: r4(4), th4(4), phi4(4), fr4(4), ft4(4), fp4(4)
00058 real(long) :: fr4psi(4) , ft4psi(4) , fp4psi(4)
00059 real(long) :: fr4alph(4) , ft4alph(4) , fp4alph(4)
00060 real(long) :: fr4bvxd(4) , ft4bvxd(4) , fp4bvxd(4)
00061 real(long) :: fr4bvyd(4) , ft4bvyd(4) , fp4bvyd(4)
00062 real(long) :: fr4bvzd(4) , ft4bvzd(4) , fp4bvzd(4)
00063 real(long) :: fr4axx(4) , ft4axx(4) , fp4axx(4)
00064 real(long) :: fr4axy(4) , ft4axy(4) , fp4axy(4)
00065 real(long) :: fr4axz(4) , ft4axz(4) , fp4axz(4)
00066 real(long) :: fr4ayy(4) , ft4ayy(4) , fp4ayy(4)
00067 real(long) :: fr4ayz(4) , ft4ayz(4) , fp4ayz(4)
00068 real(long) :: fr4azz(4) , ft4azz(4) , fp4azz(4)
00069 real(long) :: fr4emd(4) , ft4emd(4) , fp4emd(4)
00070
00071 real(long) :: fr4psif(4) , ft4psif(4) , fp4psif(4)
00072 real(long) :: fr4alphf(4) , ft4alphf(4) , fp4alphf(4)
00073 real(long) :: fr4bvxdf(4) , ft4bvxdf(4) , fp4bvxdf(4)
00074 real(long) :: fr4bvydf(4) , ft4bvydf(4) , fp4bvydf(4)
00075 real(long) :: fr4bvzdf(4) , ft4bvzdf(4) , fp4bvzdf(4)
00076
00077 integer :: irg, itg, ipg, irgex, itgex, ipgex
00078 integer :: ir0, it0, ip0, irg0 , itg0 , ipg0, ii, jj, kk
00079 real(long), external :: lagint_4th
00080
00081 real(8), pointer :: rg_p1(:), rgex_p1(:), thgex_p1(:), phigex_p1(:)
00082 integer, pointer :: irgex_r_p1(:), itgex_th_p1(:), ipgex_phi_p1(:)
00083 integer, pointer :: itgex_r_p1(:,:), ipgex_r_p1(:,:), ipgex_th_p1(:,:)
00084 real(8), pointer :: emd_p1(:,:,:), rs_p1(:,:)
00085 real(8), pointer :: psif_p1(:,:,:), alphf_p1(:,:,:), bvxdf_p1(:,:,:), bvydf_p1(:,:,:), bvzdf_p1(:,:,:)
00086 real(8), pointer :: psi_p1(:,:,:), alph_p1(:,:,:), bvxd_p1(:,:,:), bvyd_p1(:,:,:), bvzd_p1(:,:,:)
00087 real(8), pointer :: axx_p1(:,:,:), axy_p1(:,:,:) , axz_p1(:,:,:) , ayy_p1(:,:,:) , ayz_p1(:,:,:), azz_p1(:,:,:)
00088
00089 real(8), pointer :: rg_p2(:), rgex_p2(:), thgex_p2(:), phigex_p2(:)
00090 integer, pointer :: irgex_r_p2(:), itgex_th_p2(:), ipgex_phi_p2(:)
00091 integer, pointer :: itgex_r_p2(:,:), ipgex_r_p2(:,:), ipgex_th_p2(:,:)
00092 real(8), pointer :: emd_p2(:,:,:), rs_p2(:,:)
00093 real(8), pointer :: psif_p2(:,:,:), alphf_p2(:,:,:), bvxdf_p2(:,:,:), bvydf_p2(:,:,:), bvzdf_p2(:,:,:)
00094 real(8), pointer :: psi_p2(:,:,:), alph_p2(:,:,:), bvxd_p2(:,:,:), bvyd_p2(:,:,:), bvzd_p2(:,:,:)
00095 real(8), pointer :: axx_p2(:,:,:), axy_p2(:,:,:) , axz_p2(:,:,:) , ayy_p2(:,:,:) , ayz_p2(:,:,:), azz_p2(:,:,:)
00096
00097 real(8), pointer :: rg_p3(:), rgex_p3(:), thgex_p3(:), phigex_p3(:)
00098 integer, pointer :: irgex_r_p3(:), itgex_th_p3(:), ipgex_phi_p3(:)
00099 integer, pointer :: itgex_r_p3(:,:), ipgex_r_p3(:,:), ipgex_th_p3(:,:)
00100 real(8), pointer :: psi_p3(:,:,:), alph_p3(:,:,:), bvxd_p3(:,:,:), bvyd_p3(:,:,:), bvzd_p3(:,:,:)
00101 real(8), pointer :: axx_p3(:,:,:), axy_p3(:,:,:) , axz_p3(:,:,:) , ayy_p3(:,:,:) , ayz_p3(:,:,:), azz_p3(:,:,:)
00102
00103 gxx=0.0d0; gxy=0.0d0; gxz=0.0d0; gyy=0.0d0; gyz=0.0d0; gzz=0.0d0
00104 kxx=0.0d0; kxy=0.0d0; kxz=0.0d0; kyy=0.0d0; kyz=0.0d0; kzz=0.0d0
00105 axx=0.0d0; axy=0.0d0; axz=0.0d0; ayy=0.0d0; ayz=0.0d0; azz=0.0d0
00106
00107
00108
00109
00110 dir_path='.'
00111
00112
00113 call allocate_grid_parameter_mpt
00114 call allocate_grid_parameter_binary_excision_mpt
00115 call allocate_def_matter_parameter_mpt
00116 do impt = 1, nmpt
00117 call read_parameter_mpt_cactus(impt,dir_path)
00118 if (impt .le. 2) call read_surf_parameter_mpt_cactus(impt,dir_path)
00119 call copy_grid_parameter_to_mpt(impt)
00120 call read_parameter_binary_excision_mpt_cactus(impt,dir_path)
00121 call copy_grid_parameter_binary_excision_to_mpt(impt)
00122 if (impt .le. 2) call peos_initialize_mpt_cactus(impt,dir_path)
00123 call copy_def_peos_parameter_to_mpt(impt)
00124 end do
00125
00126 call set_allocate_size_mpt
00127 call allocate_trig_grav_mphi
00128 call allocate_trigonometry_grav_phi_mpt
00129
00130 do impt = 1, nmpt
00131 call copy_grid_parameter_from_mpt(impt)
00132 call copy_grid_parameter_binary_excision_from_mpt(impt)
00133 call copy_def_peos_parameter_from_mpt(impt)
00134 call coordinate_patch_kit_grav_grid_coc2cac_mpt(3)
00135 call calc_parameter_binary_excision
00136 call copy_grid_parameter_to_mpt(impt)
00137 call copy_grid_parameter_binary_excision_to_mpt(impt)
00138 call copy_coordinate_grav_extended_to_mpt(impt)
00139 call copy_coordinate_grav_phi_to_mpt(impt)
00140 call copy_coordinate_grav_r_to_mpt(impt)
00141 call copy_coordinate_grav_theta_to_mpt(impt)
00142 call copy_def_binary_parameter_to_mpt(impt)
00143 call copy_trigonometry_grav_phi_to_mpt(impt)
00144 call copy_trigonometry_grav_theta_to_mpt(impt)
00145 end do
00146
00147 do impt = 1,nmpt
00148
00149 call copy_grid_parameter_from_mpt(impt)
00150 call copy_grid_parameter_binary_excision_from_mpt(impt)
00151 call copy_coordinate_grav_extended_from_mpt(impt)
00152 call copy_coordinate_grav_phi_from_mpt(impt)
00153 call copy_coordinate_grav_r_from_mpt(impt)
00154 call copy_coordinate_grav_theta_from_mpt(impt)
00155 call copy_def_binary_parameter_from_mpt(impt)
00156 call copy_trigonometry_grav_phi_from_mpt(impt)
00157 call copy_trigonometry_grav_theta_from_mpt(impt)
00158
00159
00160 if (impt==1) then
00161 nrg_p1=nrg; ntg_p1=ntg; npg_p1=npg; nrf_p1=nrf; ntf_p1=ntf; npf_p1=npf
00162 allocate ( rg_p1( 0:nnrg))
00163 allocate ( rgex_p1(-2:nnrg+2))
00164 allocate ( thgex_p1(-2:nntg+2))
00165 allocate ( phigex_p1(-2:nnpg+2))
00166 allocate ( irgex_r_p1(-2:nnrg+2))
00167 allocate ( itgex_th_p1(-2:nntg+2))
00168 allocate (ipgex_phi_p1(-2:nnpg+2))
00169
00170 allocate ( itgex_r_p1(0:nntg,-2:nnrg+2))
00171 allocate ( ipgex_r_p1(0:nnpg,-2:nnrg+2))
00172 allocate (ipgex_th_p1(0:nnpg,-2:nntg+2))
00173
00174 rg_p1=rg; rgex_p1=rgex; thgex_p1=thgex; phigex_p1=phigex
00175 irgex_r_p1=irgex_r; itgex_th_p1=itgex_th; ipgex_phi_p1=ipgex_phi
00176 itgex_r_p1=itgex_r; ipgex_r_p1=ipgex_r; ipgex_th_p1=ipgex_th
00177
00178 rr3 = 0.7d0*(rgout - rgmid)
00179 dis_cm = dis
00180 r_surf_p1 = r_surf
00181 write(6,'(a13,i1,a8)') "Reading COCP-", impt, " data..."
00182
00183 allocate ( emd_p1(0:nrf,0:ntf,0:npf))
00184 allocate ( psif_p1(0:nrf,0:ntf,0:npf))
00185 allocate (alphf_p1(0:nrf,0:ntf,0:npf))
00186 allocate (bvxdf_p1(0:nrf,0:ntf,0:npf))
00187 allocate (bvydf_p1(0:nrf,0:ntf,0:npf))
00188 allocate (bvzdf_p1(0:nrf,0:ntf,0:npf))
00189 allocate ( rs_p1(0:ntf,0:npf))
00190 allocate ( psi_p1(0:nrg,0:ntg,0:npg))
00191 allocate ( alph_p1(0:nrg,0:ntg,0:npg))
00192 allocate ( bvxd_p1(0:nrg,0:ntg,0:npg))
00193 allocate ( bvyd_p1(0:nrg,0:ntg,0:npg))
00194 allocate ( bvzd_p1(0:nrg,0:ntg,0:npg))
00195 allocate ( axx_p1(0:nrg,0:ntg,0:npg))
00196 allocate ( axy_p1(0:nrg,0:ntg,0:npg))
00197 allocate ( axz_p1(0:nrg,0:ntg,0:npg))
00198 allocate ( ayy_p1(0:nrg,0:ntg,0:npg))
00199 allocate ( ayz_p1(0:nrg,0:ntg,0:npg))
00200 allocate ( azz_p1(0:nrg,0:ntg,0:npg))
00201 emd_p1=0.0d0; rs_p1 =0.0d0
00202 psi_p1=0.0d0; alph_p1=0.0d0; bvxd_p1=0.0d0; bvyd_p1=0.0d0; bvzd_p1=0.0d0
00203 axx_p1=0.0d0; axy_p1 =0.0d0; axz_p1 =0.0d0; ayy_p1=0.0d0; ayz_p1=0.0d0; azz_p1=0.0d0
00204
00205 call IO_input_CF_grav_export(trim(dir_path)//"/bnsgra_3D_mpt1.las",psi_p1,alph_p1,bvxd_p1,bvyd_p1,bvzd_p1)
00206
00207 call IO_input_CF_flco_export(trim(dir_path)//"/bnsflu_3D_mpt1.las",emd_p1,ome_p1,ber_p1,radi_p1)
00208
00209 call IO_input_CF_surf_export(trim(dir_path)//"/bnssur_3D_mpt1.las",rs_p1)
00210
00211 call excurve_CF_gridpoint_export(alph_p1,bvxd_p1,bvyd_p1,bvzd_p1, &
00212 & axx_p1, axy_p1, axz_p1, ayy_p1, ayz_p1, azz_p1)
00213
00214 call interpo_gr2fl_metric_CF_export(alph_p1, psi_p1, bvxd_p1, bvyd_p1, bvzd_p1, &
00215 & alphf_p1, psif_p1, bvxdf_p1, bvydf_p1, bvzdf_p1, rs_p1)
00216
00217 end if
00218 if (impt==2) then
00219 nrg_p2=nrg; ntg_p2=ntg; npg_p2=npg; nrf_p2=nrf; ntf_p2=ntf; npf_p2=npf
00220 allocate ( rg_p2( 0:nnrg))
00221 allocate ( rgex_p2(-2:nnrg+2))
00222 allocate ( thgex_p2(-2:nntg+2))
00223 allocate ( phigex_p2(-2:nnpg+2))
00224 allocate ( irgex_r_p2(-2:nnrg+2))
00225 allocate ( itgex_th_p2(-2:nntg+2))
00226 allocate (ipgex_phi_p2(-2:nnpg+2))
00227
00228 allocate ( itgex_r_p2(0:nntg,-2:nnrg+2))
00229 allocate ( ipgex_r_p2(0:nnpg,-2:nnrg+2))
00230 allocate (ipgex_th_p2(0:nnpg,-2:nntg+2))
00231
00232 rg_p2=rg; rgex_p2=rgex; thgex_p2=thgex; phigex_p2=phigex
00233 irgex_r_p2=irgex_r; itgex_th_p2=itgex_th; ipgex_phi_p2=ipgex_phi
00234 itgex_r_p2=itgex_r; ipgex_r_p2=ipgex_r; ipgex_th_p2=ipgex_th
00235
00236 r_surf_p2 = r_surf
00237 write(6,'(a13,i1,a8)') "Reading COCP-", impt, " data..."
00238
00239 allocate ( emd_p2(0:nrf,0:ntf,0:npf))
00240 allocate ( psif_p2(0:nrf,0:ntf,0:npf))
00241 allocate (alphf_p2(0:nrf,0:ntf,0:npf))
00242 allocate (bvxdf_p2(0:nrf,0:ntf,0:npf))
00243 allocate (bvydf_p2(0:nrf,0:ntf,0:npf))
00244 allocate (bvzdf_p2(0:nrf,0:ntf,0:npf))
00245 allocate ( rs_p2(0:ntf,0:npf))
00246 allocate ( psi_p2(0:nrg,0:ntg,0:npg))
00247 allocate ( alph_p2(0:nrg,0:ntg,0:npg))
00248 allocate ( bvxd_p2(0:nrg,0:ntg,0:npg))
00249 allocate ( bvyd_p2(0:nrg,0:ntg,0:npg))
00250 allocate ( bvzd_p2(0:nrg,0:ntg,0:npg))
00251 allocate ( axx_p2(0:nrg,0:ntg,0:npg))
00252 allocate ( axy_p2(0:nrg,0:ntg,0:npg))
00253 allocate ( axz_p2(0:nrg,0:ntg,0:npg))
00254 allocate ( ayy_p2(0:nrg,0:ntg,0:npg))
00255 allocate ( ayz_p2(0:nrg,0:ntg,0:npg))
00256 allocate ( azz_p2(0:nrg,0:ntg,0:npg))
00257 emd_p2=0.0d0; rs_p2=0.0d0
00258 psi_p2=0.0d0; alph_p2=0.0d0; bvxd_p2=0.0d0; bvyd_p2=0.0d0; bvzd_p2=0.0d0
00259 axx_p2=0.0d0; axy_p2=0.0d0; axz_p2=0.0d0; ayy_p2=0.0d0; ayz_p2=0.0d0; azz_p2=0.0d0
00260
00261 call IO_input_CF_grav_export(trim(dir_path)//"/bnsgra_3D_mpt2.las",psi_p2,alph_p2,bvxd_p2,bvyd_p2,bvzd_p2)
00262
00263 call IO_input_CF_flco_export(trim(dir_path)//"/bnsflu_3D_mpt2.las",emd_p2,ome_p2,ber_p2,radi_p2)
00264
00265 call IO_input_CF_surf_export(trim(dir_path)//"/bnssur_3D_mpt2.las",rs_p2)
00266
00267 call excurve_CF_gridpoint_export(alph_p2,bvxd_p2,bvyd_p2,bvzd_p2, &
00268 & axx_p2, axy_p2, axz_p2, ayy_p2, ayz_p2, azz_p2)
00269
00270 call interpo_gr2fl_metric_CF_export(alph_p2, psi_p2, bvxd_p2, bvyd_p2, bvzd_p2, &
00271 & alphf_p2, psif_p2, bvxdf_p2, bvydf_p2, bvzdf_p2, rs_p2)
00272
00273 bvxdf_p2 = - bvxdf_p2
00274 bvydf_p2 = - bvydf_p2
00275 bvxd_p2 = - bvxd_p2
00276 bvyd_p2 = - bvyd_p2
00277 axz_p2 = - axz_p2
00278 ayz_p2 = - ayz_p2
00279 end if
00280 if (impt==3) then
00281 nrg_p3=nrg; ntg_p3=ntg; npg_p3=npg; nrf_p3=nrf; ntf_p3=ntf; npf_p3=npf
00282 allocate ( rg_p3( 0:nnrg))
00283 allocate ( rgex_p3(-2:nnrg+2))
00284 allocate ( thgex_p3(-2:nntg+2))
00285 allocate ( phigex_p3(-2:nnpg+2))
00286 allocate ( irgex_r_p3(-2:nnrg+2))
00287 allocate ( itgex_th_p3(-2:nntg+2))
00288 allocate (ipgex_phi_p3(-2:nnpg+2))
00289
00290 allocate ( itgex_r_p3(0:nntg,-2:nnrg+2))
00291 allocate ( ipgex_r_p3(0:nnpg,-2:nnrg+2))
00292 allocate (ipgex_th_p3(0:nnpg,-2:nntg+2))
00293
00294 rg_p3=rg; rgex_p3=rgex; thgex_p3=thgex; phigex_p3=phigex
00295 irgex_r_p3=irgex_r; itgex_th_p3=itgex_th; ipgex_phi_p3=ipgex_phi
00296 itgex_r_p3=itgex_r; ipgex_r_p3=ipgex_r; ipgex_th_p3=ipgex_th
00297
00298 write(6,'(a13,i1,a8)') "Reading ARCP-", impt, " data..."
00299 allocate ( psi_p3(0:nrg,0:ntg,0:npg))
00300 allocate (alph_p3(0:nrg,0:ntg,0:npg))
00301 allocate (bvxd_p3(0:nrg,0:ntg,0:npg))
00302 allocate (bvyd_p3(0:nrg,0:ntg,0:npg))
00303 allocate (bvzd_p3(0:nrg,0:ntg,0:npg))
00304 allocate ( axx_p3(0:nrg,0:ntg,0:npg))
00305 allocate ( axy_p3(0:nrg,0:ntg,0:npg))
00306 allocate ( axz_p3(0:nrg,0:ntg,0:npg))
00307 allocate ( ayy_p3(0:nrg,0:ntg,0:npg))
00308 allocate ( ayz_p3(0:nrg,0:ntg,0:npg))
00309 allocate ( azz_p3(0:nrg,0:ntg,0:npg))
00310 psi_p3=0.0d0; alph_p3=0.0d0; bvxd_p3=0.0d0; bvyd_p3=0.0d0; bvzd_p3=0.0d0
00311 axx_p3=0.0d0; axy_p3=0.0d0; axz_p3=0.0d0; ayy_p3=0.0d0; ayz_p3=0.0d0; azz_p3=0.0d0
00312
00313 call IO_input_CF_grav_export(trim(dir_path)//"/bnsgra_3D_mpt3.las",psi_p3,alph_p3,bvxd_p3,bvyd_p3,bvzd_p3)
00314
00315 call excurve_CF_gridpoint_export(alph_p3,bvxd_p3,bvyd_p3,bvzd_p3, &
00316 & axx_p3, axy_p3, axz_p3, ayy_p3, ayz_p3, azz_p3)
00317
00318 end if
00319 end do
00320 write(6,'(2e20.12)') emd_p1(0,0,0), emd_p1(58,0,0)
00321 write(6,'(3e20.12)') ome_p1, ber_p1, radi_p1
00322 write(6,'(e20.12)') dis_cm
00323
00324 write(6,'(a56)', ADVANCE = "NO") "Give cartesian coordinates (x,y,z) separated by a space:"
00325 read(5,*) xcac,ycac,zcac
00326 write(6,'(a23,3e20.12)') "Point given wrt CACTUS:", xcac,ycac,zcac
00327 write(6,'(a38,2e20.12)') "Cocal radius scale in COCP-1, COCP-2 :", radi_p1, radi_p2
00328 write(6,'(a38,2e20.12)') "Cocal surface scale in COCP-1, COCP-2:", r_surf_p1, r_surf_p2
00329 xcoc = xcac/(radi_p1)
00330 ycoc = ycac/(radi_p1)
00331 zcoc = zcac/(radi_p1)
00332 write(6,'(a23,3e20.12)') "Point given wrt COCAL:", xcoc,ycoc,zcoc
00333
00334 rrcm = sqrt(xcoc**2 + ycoc**2 + zcoc**2)
00335 write(6,*) rrcm, rr3
00336 if (rrcm > rr3) then
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 xc_p3 = xcoc
00349 yc_p3 = ycoc
00350 zc_p3 = zcoc
00351 write(6,'(a39,3e20.12)') "Point given wrt COCP-3 in Cocal scale :", xc_p3,yc_p3,zc_p3
00352
00353
00354 psica=0.0d0; alphca=0.0d0; bvxdca=0.0d0; bvydca=0.0d0; bvzdca=0.0d0
00355 axx=0.0d0 ; axy=0.0d0 ; axz=0.0d0 ; ayy=0.0d0 ; ayz=0.0d0 ; azz=0.0d0
00356
00357 rc_p3 = dsqrt(dabs(xc_p3**2 + yc_p3**2 + zc_p3**2))
00358 varpic_p3 = dsqrt(dabs(xc_p3**2 + yc_p3**2))
00359 thc_p3 = dmod(2.0d0*pi + datan2(varpic_p3,zc_p3),2.0d0*pi)
00360 phic_p3 = dmod(2.0d0*pi + datan2( yc_p3,xc_p3),2.0d0*pi)
00361
00362 do irg = 0, nrg_p3+1
00363 if (rc_p3.lt.rgex_p3(irg).and.rc_p3.ge.rgex_p3(irg-1)) ir0 = min0(irg-2,nrg_p3-3)
00364 end do
00365 do itg = 0, ntg_p3+1
00366 if (thc_p3.lt.thgex_p3(itg).and.thc_p3.ge.thgex_p3(itg-1)) it0 = itg-2
00367 end do
00368 do ipg = 0, npg_p3+1
00369 if (phic_p3.lt.phigex_p3(ipg).and.phic_p3.ge.phigex_p3(ipg-1)) ip0 = ipg-2
00370 end do
00371
00372 do ii = 1, 4
00373 irg0 = ir0 + ii - 1
00374 itg0 = it0 + ii - 1
00375 ipg0 = ip0 + ii - 1
00376 r4(ii) = rgex_p3(irg0)
00377 th4(ii) = thgex_p3(itg0)
00378 phi4(ii) = phigex_p3(ipg0)
00379 end do
00380
00381 do kk = 1, 4
00382 ipg0 = ip0 + kk - 1
00383 do jj = 1, 4
00384 itg0 = it0 + jj - 1
00385 do ii = 1, 4
00386 irg0 = ir0 + ii - 1
00387 irgex = irgex_r_p3(irg0)
00388 itgex = itgex_r_p3(itgex_th_p3(itg0),irg0)
00389 ipgex = ipgex_r_p3(ipgex_th_p3(ipgex_phi_p3(ipg0),itg0),irg0)
00390 fr4psi(ii) = psi_p3(irgex,itgex,ipgex)
00391 fr4alph(ii) = alph_p3(irgex,itgex,ipgex)
00392 fr4bvxd(ii) = bvxd_p3(irgex,itgex,ipgex)
00393 fr4bvyd(ii) = bvyd_p3(irgex,itgex,ipgex)
00394 fr4bvzd(ii) = bvzd_p3(irgex,itgex,ipgex)
00395 fr4axx(ii) = axx_p3(irgex,itgex,ipgex)
00396 fr4axy(ii) = axy_p3(irgex,itgex,ipgex)
00397 fr4axz(ii) = axz_p3(irgex,itgex,ipgex)
00398 fr4ayy(ii) = ayy_p3(irgex,itgex,ipgex)
00399 fr4ayz(ii) = ayz_p3(irgex,itgex,ipgex)
00400 fr4azz(ii) = azz_p3(irgex,itgex,ipgex)
00401 end do
00402 ft4psi(jj) = lagint_4th(r4,fr4psi ,rc_p3)
00403 ft4alph(jj) = lagint_4th(r4,fr4alph,rc_p3)
00404 ft4bvxd(jj) = lagint_4th(r4,fr4bvxd,rc_p3)
00405 ft4bvyd(jj) = lagint_4th(r4,fr4bvyd,rc_p3)
00406 ft4bvzd(jj) = lagint_4th(r4,fr4bvzd,rc_p3)
00407 ft4axx(jj) = lagint_4th(r4,fr4axx ,rc_p3)
00408 ft4axy(jj) = lagint_4th(r4,fr4axy ,rc_p3)
00409 ft4axz(jj) = lagint_4th(r4,fr4axz ,rc_p3)
00410 ft4ayy(jj) = lagint_4th(r4,fr4ayy ,rc_p3)
00411 ft4ayz(jj) = lagint_4th(r4,fr4ayz ,rc_p3)
00412 ft4azz(jj) = lagint_4th(r4,fr4azz ,rc_p3)
00413 end do
00414 fp4psi(kk) = lagint_4th(th4,ft4psi ,thc_p3)
00415 fp4alph(kk) = lagint_4th(th4,ft4alph,thc_p3)
00416 fp4bvxd(kk) = lagint_4th(th4,ft4bvxd,thc_p3)
00417 fp4bvyd(kk) = lagint_4th(th4,ft4bvyd,thc_p3)
00418 fp4bvzd(kk) = lagint_4th(th4,ft4bvzd,thc_p3)
00419 fp4axx(kk) = lagint_4th(th4,ft4axx ,thc_p3)
00420 fp4axy(kk) = lagint_4th(th4,ft4axy ,thc_p3)
00421 fp4axz(kk) = lagint_4th(th4,ft4axz ,thc_p3)
00422 fp4ayy(kk) = lagint_4th(th4,ft4ayy ,thc_p3)
00423 fp4ayz(kk) = lagint_4th(th4,ft4ayz ,thc_p3)
00424 fp4azz(kk) = lagint_4th(th4,ft4azz ,thc_p3)
00425 end do
00426 psica = lagint_4th(phi4,fp4psi ,phic_p3)
00427 alphca = lagint_4th(phi4,fp4alph,phic_p3)
00428 bvxdca = lagint_4th(phi4,fp4bvxd,phic_p3)
00429 bvydca = lagint_4th(phi4,fp4bvyd,phic_p3)
00430 bvzdca = lagint_4th(phi4,fp4bvzd,phic_p3)
00431 axx = lagint_4th(phi4,fp4axx ,phic_p3)
00432 axy = lagint_4th(phi4,fp4axy ,phic_p3)
00433 axz = lagint_4th(phi4,fp4axz ,phic_p3)
00434 ayy = lagint_4th(phi4,fp4ayy ,phic_p3)
00435 ayz = lagint_4th(phi4,fp4ayz ,phic_p3)
00436 azz = lagint_4th(phi4,fp4azz ,phic_p3)
00437
00438 psi4ca = psica**4
00439 vxu = 0.0d0
00440 vyu = 0.0d0
00441 vzu = 0.0d0
00442 emdca = 0.0d0
00443 else
00444 if (xcoc<=0.0d0) then
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456 xc_p1 = xcoc + dis_cm
00457 yc_p1 = ycoc
00458 zc_p1 = zcoc
00459 write(6,'(a39,3e20.12)') "Point given wrt COCP-1 in Cocal scale :", xc_p1,yc_p1,zc_p1
00460
00461 psica=0.0d0; alphca=0.0d0; bvxdca=0.0d0; bvydca=0.0d0; bvzdca=0.0d0
00462 axx=0.0d0 ; axy=0.0d0 ; axz=0.0d0 ; ayy=0.0d0 ; ayz=0.0d0 ; azz=0.0d0
00463
00464 rc_p1 = dsqrt(dabs(xc_p1**2 + yc_p1**2 + zc_p1**2))
00465 varpic_p1 = dsqrt(dabs(xc_p1**2 + yc_p1**2))
00466 thc_p1 = dmod(2.0d0*pi + datan2(varpic_p1,zc_p1),2.0d0*pi)
00467 phic_p1 = dmod(2.0d0*pi + datan2( yc_p1,xc_p1),2.0d0*pi)
00468
00469 do irg = 0, nrg_p1+1
00470 if (rc_p1.lt.rgex_p1(irg).and.rc_p1.ge.rgex_p1(irg-1)) ir0 = min0(irg-2,nrg_p1-3)
00471 end do
00472 do itg = 0, ntg_p1+1
00473 if (thc_p1.lt.thgex_p1(itg).and.thc_p1.ge.thgex_p1(itg-1)) it0 = itg-2
00474 end do
00475 do ipg = 0, npg_p1+1
00476 if (phic_p1.lt.phigex_p1(ipg).and.phic_p1.ge.phigex_p1(ipg-1)) ip0 = ipg-2
00477 end do
00478
00479
00480
00481 do ii = 1, 4
00482 irg0 = ir0 + ii - 1
00483 itg0 = it0 + ii - 1
00484 ipg0 = ip0 + ii - 1
00485 r4(ii) = rgex_p1(irg0)
00486 th4(ii) = thgex_p1(itg0)
00487 phi4(ii) = phigex_p1(ipg0)
00488 end do
00489
00490 do kk = 1, 4
00491 ipg0 = ip0 + kk - 1
00492 do jj = 1, 4
00493 itg0 = it0 + jj - 1
00494 do ii = 1, 4
00495 irg0 = ir0 + ii - 1
00496 irgex = irgex_r_p1(irg0)
00497 itgex = itgex_r_p1(itgex_th_p1(itg0),irg0)
00498 ipgex = ipgex_r_p1(ipgex_th_p1(ipgex_phi_p1(ipg0),itg0),irg0)
00499 fr4psi(ii) = psi_p1(irgex,itgex,ipgex)
00500 fr4alph(ii) = alph_p1(irgex,itgex,ipgex)
00501 fr4bvxd(ii) = bvxd_p1(irgex,itgex,ipgex)
00502 fr4bvyd(ii) = bvyd_p1(irgex,itgex,ipgex)
00503 fr4bvzd(ii) = bvzd_p1(irgex,itgex,ipgex)
00504 fr4axx(ii) = axx_p1(irgex,itgex,ipgex)
00505 fr4axy(ii) = axy_p1(irgex,itgex,ipgex)
00506 fr4axz(ii) = axz_p1(irgex,itgex,ipgex)
00507 fr4ayy(ii) = ayy_p1(irgex,itgex,ipgex)
00508 fr4ayz(ii) = ayz_p1(irgex,itgex,ipgex)
00509 fr4azz(ii) = azz_p1(irgex,itgex,ipgex)
00510 end do
00511 ft4psi(jj) = lagint_4th(r4,fr4psi ,rc_p1)
00512 ft4alph(jj) = lagint_4th(r4,fr4alph,rc_p1)
00513 ft4bvxd(jj) = lagint_4th(r4,fr4bvxd,rc_p1)
00514 ft4bvyd(jj) = lagint_4th(r4,fr4bvyd,rc_p1)
00515 ft4bvzd(jj) = lagint_4th(r4,fr4bvzd,rc_p1)
00516 ft4axx(jj) = lagint_4th(r4,fr4axx ,rc_p1)
00517 ft4axy(jj) = lagint_4th(r4,fr4axy ,rc_p1)
00518 ft4axz(jj) = lagint_4th(r4,fr4axz ,rc_p1)
00519 ft4ayy(jj) = lagint_4th(r4,fr4ayy ,rc_p1)
00520 ft4ayz(jj) = lagint_4th(r4,fr4ayz ,rc_p1)
00521 ft4azz(jj) = lagint_4th(r4,fr4azz ,rc_p1)
00522 end do
00523 fp4psi(kk) = lagint_4th(th4,ft4psi ,thc_p1)
00524 fp4alph(kk) = lagint_4th(th4,ft4alph,thc_p1)
00525 fp4bvxd(kk) = lagint_4th(th4,ft4bvxd,thc_p1)
00526 fp4bvyd(kk) = lagint_4th(th4,ft4bvyd,thc_p1)
00527 fp4bvzd(kk) = lagint_4th(th4,ft4bvzd,thc_p1)
00528 fp4axx(kk) = lagint_4th(th4,ft4axx ,thc_p1)
00529 fp4axy(kk) = lagint_4th(th4,ft4axy ,thc_p1)
00530 fp4axz(kk) = lagint_4th(th4,ft4axz ,thc_p1)
00531 fp4ayy(kk) = lagint_4th(th4,ft4ayy ,thc_p1)
00532 fp4ayz(kk) = lagint_4th(th4,ft4ayz ,thc_p1)
00533 fp4azz(kk) = lagint_4th(th4,ft4azz ,thc_p1)
00534 end do
00535 psica = lagint_4th(phi4,fp4psi ,phic_p1)
00536 alphca = lagint_4th(phi4,fp4alph,phic_p1)
00537 bvxdca = lagint_4th(phi4,fp4bvxd,phic_p1)
00538 bvydca = lagint_4th(phi4,fp4bvyd,phic_p1)
00539 bvzdca = lagint_4th(phi4,fp4bvzd,phic_p1)
00540 axx = lagint_4th(phi4,fp4axx ,phic_p1)
00541 axy = lagint_4th(phi4,fp4axy ,phic_p1)
00542 axz = lagint_4th(phi4,fp4axz ,phic_p1)
00543 ayy = lagint_4th(phi4,fp4ayy ,phic_p1)
00544 ayz = lagint_4th(phi4,fp4ayz ,phic_p1)
00545 azz = lagint_4th(phi4,fp4azz ,phic_p1)
00546
00547 psi4ca = psica**4
00548
00549
00550 call interpo_lag4th_2Dsurf(rsca_p1,rs_p1,thc_p1,phic_p1)
00551 rcf_p1 = rc_p1/rsca_p1
00552
00553 if (rcf_p1.le.rg_p1(nrf_p1)) then
00554 do irg = 0, nrf_p1+1
00555 if (rcf_p1.lt.rgex_p1(irg).and.rcf_p1.ge.rgex_p1(irg-1)) ir0 = min0(irg-2,nrf_p1-3)
00556 end do
00557
00558 do ii = 1, 4
00559 irg0 = ir0 + ii - 1
00560 r4(ii) = rgex_p1(irg0)
00561 end do
00562
00563 do kk = 1, 4
00564 ipg0 = ip0 + kk - 1
00565 do jj = 1, 4
00566 itg0 = it0 + jj - 1
00567 do ii = 1, 4
00568 irg0 = ir0 + ii - 1
00569 irgex = irgex_r_p1(irg0)
00570 itgex = itgex_r_p1(itgex_th_p1(itg0),irg0)
00571 ipgex = ipgex_r_p1(ipgex_th_p1(ipgex_phi_p1(ipg0),itg0),irg0)
00572 fr4emd(ii) = emd_p1(irgex,itgex,ipgex)
00573
00574 fr4psif(ii) = psif_p1(irgex,itgex,ipgex)
00575 fr4alphf(ii) = alphf_p1(irgex,itgex,ipgex)
00576 fr4bvxdf(ii) = bvxdf_p1(irgex,itgex,ipgex)
00577 fr4bvydf(ii) = bvydf_p1(irgex,itgex,ipgex)
00578 fr4bvzdf(ii) = bvzdf_p1(irgex,itgex,ipgex)
00579 end do
00580 ft4emd(jj) = lagint_4th(r4,fr4emd ,rcf_p1)
00581
00582 ft4psif(jj) = lagint_4th(r4,fr4psif ,rcf_p1)
00583 ft4alphf(jj) = lagint_4th(r4,fr4alphf,rcf_p1)
00584 ft4bvxdf(jj) = lagint_4th(r4,fr4bvxdf,rcf_p1)
00585 ft4bvydf(jj) = lagint_4th(r4,fr4bvydf,rcf_p1)
00586 ft4bvzdf(jj) = lagint_4th(r4,fr4bvzdf,rcf_p1)
00587 end do
00588 fp4emd(kk) = lagint_4th(th4,ft4emd ,thc_p1)
00589
00590 fp4psif(kk) = lagint_4th(th4,ft4psif ,thc_p1)
00591 fp4alphf(kk) = lagint_4th(th4,ft4alphf,thc_p1)
00592 fp4bvxdf(kk) = lagint_4th(th4,ft4bvxdf,thc_p1)
00593 fp4bvydf(kk) = lagint_4th(th4,ft4bvydf,thc_p1)
00594 fp4bvzdf(kk) = lagint_4th(th4,ft4bvzdf,thc_p1)
00595 end do
00596 emdca = lagint_4th(phi4,fp4emd ,phic_p1)
00597
00598 psifca = lagint_4th(phi4,fp4psif ,phic_p1)
00599 alphfca = lagint_4th(phi4,fp4alphf,phic_p1)
00600 bvxdfca = lagint_4th(phi4,fp4bvxdf,phic_p1)
00601 bvydfca = lagint_4th(phi4,fp4bvydf,phic_p1)
00602 bvzdfca = lagint_4th(phi4,fp4bvzdf,phic_p1)
00603
00604 psif4ca = psifca**4
00605 bxcor = bvxdfca + ome_p1*(-ycoc)
00606 bycor = bvydfca + ome_p1*(xcoc)
00607 bzcor = bvzdfca
00608
00609 vxu = bxcor/alphfca
00610 vyu = bycor/alphfca
00611 vzu = bzcor/alphfca
00612 else
00613 emdca=0.0d0
00614 vxu=0.0d0; vyu=0.0d0; vzu=0.0d0
00615 end if
00616 else
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628 xc_p2 = -(xcoc - dis_cm)
00629 yc_p2 = -ycoc
00630 zc_p2 = zcoc
00631 write(6,'(a39,3e20.12)') "Point given wrt COCP-2 in Cocal scale :", xc_p2,yc_p2,zc_p2
00632
00633 psica=0.0d0; alphca=0.0d0; bvxdca=0.0d0; bvydca=0.0d0; bvzdca=0.0d0
00634 axx=0.0d0 ; axy=0.0d0 ; axz=0.0d0 ; ayy=0.0d0 ; ayz=0.0d0 ; azz=0.0d0
00635
00636 rc_p2 = dsqrt(dabs(xc_p2**2 + yc_p2**2 + zc_p2**2))
00637 varpic_p2 = dsqrt(dabs(xc_p2**2 + yc_p2**2))
00638 thc_p2 = dmod(2.0d0*pi + datan2(varpic_p2,zc_p2),2.0d0*pi)
00639 phic_p2 = dmod(2.0d0*pi + datan2( yc_p2,xc_p2),2.0d0*pi)
00640
00641 do irg = 0, nrg_p2+1
00642 if (rc_p2.lt.rgex_p2(irg).and.rc_p2.ge.rgex_p2(irg-1)) ir0 = min0(irg-2,nrg_p2-3)
00643 end do
00644 do itg = 0, ntg_p2+1
00645 if (thc_p2.lt.thgex_p2(itg).and.thc_p2.ge.thgex_p2(itg-1)) it0 = itg-2
00646 end do
00647 do ipg = 0, npg_p2+1
00648 if (phic_p2.lt.phigex_p2(ipg).and.phic_p2.ge.phigex_p2(ipg-1)) ip0 = ipg-2
00649 end do
00650
00651 do ii = 1, 4
00652 irg0 = ir0 + ii - 1
00653 itg0 = it0 + ii - 1
00654 ipg0 = ip0 + ii - 1
00655 r4(ii) = rgex_p2(irg0)
00656 th4(ii) = thgex_p2(itg0)
00657 phi4(ii) = phigex_p2(ipg0)
00658 end do
00659
00660 do kk = 1, 4
00661 ipg0 = ip0 + kk - 1
00662 do jj = 1, 4
00663 itg0 = it0 + jj - 1
00664 do ii = 1, 4
00665 irg0 = ir0 + ii - 1
00666 irgex = irgex_r_p2(irg0)
00667 itgex = itgex_r_p2(itgex_th_p2(itg0),irg0)
00668 ipgex = ipgex_r_p2(ipgex_th_p2(ipgex_phi_p2(ipg0),itg0),irg0)
00669 fr4psi(ii) = psi_p2(irgex,itgex,ipgex)
00670 fr4alph(ii) = alph_p2(irgex,itgex,ipgex)
00671 fr4bvxd(ii) = bvxd_p2(irgex,itgex,ipgex)
00672 fr4bvyd(ii) = bvyd_p2(irgex,itgex,ipgex)
00673 fr4bvzd(ii) = bvzd_p2(irgex,itgex,ipgex)
00674 fr4axx(ii) = axx_p2(irgex,itgex,ipgex)
00675 fr4axy(ii) = axy_p2(irgex,itgex,ipgex)
00676 fr4axz(ii) = axz_p2(irgex,itgex,ipgex)
00677 fr4ayy(ii) = ayy_p2(irgex,itgex,ipgex)
00678 fr4ayz(ii) = ayz_p2(irgex,itgex,ipgex)
00679 fr4azz(ii) = azz_p2(irgex,itgex,ipgex)
00680 end do
00681 ft4psi(jj) = lagint_4th(r4,fr4psi ,rc_p2)
00682 ft4alph(jj) = lagint_4th(r4,fr4alph,rc_p2)
00683 ft4bvxd(jj) = lagint_4th(r4,fr4bvxd,rc_p2)
00684 ft4bvyd(jj) = lagint_4th(r4,fr4bvyd,rc_p2)
00685 ft4bvzd(jj) = lagint_4th(r4,fr4bvzd,rc_p2)
00686 ft4axx(jj) = lagint_4th(r4,fr4axx ,rc_p2)
00687 ft4axy(jj) = lagint_4th(r4,fr4axy ,rc_p2)
00688 ft4axz(jj) = lagint_4th(r4,fr4axz ,rc_p2)
00689 ft4ayy(jj) = lagint_4th(r4,fr4ayy ,rc_p2)
00690 ft4ayz(jj) = lagint_4th(r4,fr4ayz ,rc_p2)
00691 ft4azz(jj) = lagint_4th(r4,fr4azz ,rc_p2)
00692 end do
00693 fp4psi(kk) = lagint_4th(th4,ft4psi ,thc_p2)
00694 fp4alph(kk) = lagint_4th(th4,ft4alph,thc_p2)
00695 fp4bvxd(kk) = lagint_4th(th4,ft4bvxd,thc_p2)
00696 fp4bvyd(kk) = lagint_4th(th4,ft4bvyd,thc_p2)
00697 fp4bvzd(kk) = lagint_4th(th4,ft4bvzd,thc_p2)
00698 fp4axx(kk) = lagint_4th(th4,ft4axx ,thc_p2)
00699 fp4axy(kk) = lagint_4th(th4,ft4axy ,thc_p2)
00700 fp4axz(kk) = lagint_4th(th4,ft4axz ,thc_p2)
00701 fp4ayy(kk) = lagint_4th(th4,ft4ayy ,thc_p2)
00702 fp4ayz(kk) = lagint_4th(th4,ft4ayz ,thc_p2)
00703 fp4azz(kk) = lagint_4th(th4,ft4azz ,thc_p2)
00704 end do
00705 psica = lagint_4th(phi4,fp4psi ,phic_p2)
00706 alphca = lagint_4th(phi4,fp4alph,phic_p2)
00707 bvxdca = lagint_4th(phi4,fp4bvxd,phic_p2)
00708 bvydca = lagint_4th(phi4,fp4bvyd,phic_p2)
00709 bvzdca = lagint_4th(phi4,fp4bvzd,phic_p2)
00710 axx = lagint_4th(phi4,fp4axx ,phic_p2)
00711 axy = lagint_4th(phi4,fp4axy ,phic_p2)
00712 axz = lagint_4th(phi4,fp4axz ,phic_p2)
00713 ayy = lagint_4th(phi4,fp4ayy ,phic_p2)
00714 ayz = lagint_4th(phi4,fp4ayz ,phic_p2)
00715 azz = lagint_4th(phi4,fp4azz ,phic_p2)
00716
00717 psi4ca = psica**4
00718 call interpo_lag4th_2Dsurf(rsca_p2,rs_p2,thc_p2,phic_p2)
00719 rcf_p2 = rc_p2/rsca_p2
00720
00721 if (rcf_p2.le.rg_p2(nrf_p2)) then
00722 do irg = 0, nrf_p2+1
00723 if (rcf_p2.lt.rgex_p2(irg).and.rcf_p2.ge.rgex_p2(irg-1)) ir0 = min0(irg-2,nrf_p2-3)
00724 end do
00725
00726 do ii = 1, 4
00727 irg0 = ir0 + ii - 1
00728 r4(ii) = rgex_p2(irg0)
00729 end do
00730
00731 do kk = 1, 4
00732 ipg0 = ip0 + kk - 1
00733 do jj = 1, 4
00734 itg0 = it0 + jj - 1
00735 do ii = 1, 4
00736 irg0 = ir0 + ii - 1
00737 irgex = irgex_r_p2(irg0)
00738 itgex = itgex_r_p2(itgex_th_p2(itg0),irg0)
00739 ipgex = ipgex_r_p2(ipgex_th_p2(ipgex_phi_p2(ipg0),itg0),irg0)
00740 fr4emd(ii) = emd_p2(irgex,itgex,ipgex)
00741
00742 fr4psif(ii) = psif_p2(irgex,itgex,ipgex)
00743 fr4alphf(ii) = alphf_p2(irgex,itgex,ipgex)
00744 fr4bvxdf(ii) = bvxdf_p2(irgex,itgex,ipgex)
00745 fr4bvydf(ii) = bvydf_p2(irgex,itgex,ipgex)
00746 fr4bvzdf(ii) = bvzdf_p2(irgex,itgex,ipgex)
00747 end do
00748 ft4emd(jj) = lagint_4th(r4,fr4emd ,rcf_p2)
00749
00750 ft4psif(jj) = lagint_4th(r4,fr4psif ,rcf_p2)
00751 ft4alphf(jj) = lagint_4th(r4,fr4alphf,rcf_p2)
00752 ft4bvxdf(jj) = lagint_4th(r4,fr4bvxdf,rcf_p2)
00753 ft4bvydf(jj) = lagint_4th(r4,fr4bvydf,rcf_p2)
00754 ft4bvzdf(jj) = lagint_4th(r4,fr4bvzdf,rcf_p2)
00755 end do
00756 fp4emd(kk) = lagint_4th(th4,ft4emd ,thc_p2)
00757
00758 fp4psif(kk) = lagint_4th(th4,ft4psif ,thc_p2)
00759 fp4alphf(kk) = lagint_4th(th4,ft4alphf,thc_p2)
00760 fp4bvxdf(kk) = lagint_4th(th4,ft4bvxdf,thc_p2)
00761 fp4bvydf(kk) = lagint_4th(th4,ft4bvydf,thc_p2)
00762 fp4bvzdf(kk) = lagint_4th(th4,ft4bvzdf,thc_p2)
00763 end do
00764 emdca = lagint_4th(phi4,fp4emd ,phic_p2)
00765
00766 psifca = lagint_4th(phi4,fp4psif ,phic_p2)
00767 alphfca = lagint_4th(phi4,fp4alphf,phic_p2)
00768 bvxdfca = lagint_4th(phi4,fp4bvxdf,phic_p2)
00769 bvydfca = lagint_4th(phi4,fp4bvydf,phic_p2)
00770 bvzdfca = lagint_4th(phi4,fp4bvzdf,phic_p2)
00771
00772 psif4ca = psifca**4
00773 bxcor = bvxdfca + ome_p2*(-ycoc)
00774 bycor = bvydfca + ome_p2*(xcoc)
00775 bzcor = bvzdfca
00776
00777 vxu = bxcor/alphfca
00778 vyu = bycor/alphfca
00779 vzu = bzcor/alphfca
00780 else
00781 emdca=0.0d0
00782 vxu=0.0d0; vyu=0.0d0; vzu=0.0d0
00783 end if
00784
00785 end if
00786 end if
00787
00788 gxx = psi4ca
00789 gxy = 0.0d0
00790 gxz = 0.0d0
00791 gyy = psi4ca
00792 gyz = 0.0d0
00793 gzz = psi4ca
00794
00795 kxx = psi4ca*axx/(radi_p1)
00796 kxy = psi4ca*axy/(radi_p1)
00797 kxz = psi4ca*axz/(radi_p1)
00798 kyy = psi4ca*ayy/(radi_p1)
00799 kyz = psi4ca*ayz/(radi_p1)
00800 kzz = psi4ca*azz/(radi_p1)
00801
00802 call peos_q2hprho(emdca, hca, preca, rhoca, eneca)
00803
00804 epsca = eneca/rhoca - 1.0d0
00805
00806 write(6,'(a6,e20.12)') "psi =", psica
00807 write(6,'(a6,e20.12)') "alph =", alphca
00808 write(6,'(a6,e20.12)') "bvxd =", bvxdca
00809 write(6,'(a6,e20.12)') "bvyd =", bvydca
00810 write(6,'(a6,e20.12)') "bvzd =", bvzdca
00811 write(6,'(a6,e20.12)') "Radi =", r_surf_p1*radi_p1
00812 write(6,'(a6,e20.12)') "Omeg =", ome_p1/radi_p1
00813 write(6,'(a6,e20.12)') "emd =", emdca
00814 write(6,'(a6,e20.12)') "h =", hca
00815 write(6,'(a6,e20.12)') "pre =", preca
00816 write(6,'(a6,e20.12)') "rho =", rhoca
00817 write(6,'(a6,e20.12)') "ene =", eneca
00818 write(6,'(a6,e20.12)') "eps =", epsca
00819
00820 write(6,'(a18)') "Kij at gridpoints:"
00821 write(6,'(3e20.12)') kxx, kxy, kxz
00822 write(6,'(3e20.12)') kxy, kyy, kyz
00823 write(6,'(3e20.12)') kxz, kyz, kzz
00824
00825 write(6,'(a13)') "v^i eulerian:"
00826 write(6,'(a6,e20.12)') "vxu =", vxu
00827 write(6,'(a6,e20.12)') "vyu =", vyu
00828 write(6,'(a6,e20.12)') "vzu =", vzu
00829
00830 write(6,'(a16)') "Deallocating...."
00831 deallocate( rg_p1); deallocate( rg_p2); deallocate( rg_p3)
00832 deallocate( rgex_p1); deallocate( rgex_p2); deallocate( rgex_p3)
00833 deallocate( thgex_p1); deallocate( thgex_p2); deallocate( thgex_p3)
00834 deallocate( phigex_p1); deallocate( phigex_p2); deallocate( phigex_p3)
00835 deallocate( irgex_r_p1); deallocate( irgex_r_p2); deallocate( irgex_r_p3)
00836 deallocate( itgex_th_p1); deallocate( itgex_th_p2); deallocate( itgex_th_p3)
00837 deallocate(ipgex_phi_p1); deallocate(ipgex_phi_p2); deallocate(ipgex_phi_p3)
00838 deallocate( itgex_r_p1); deallocate( itgex_r_p2); deallocate( itgex_r_p3)
00839 deallocate( ipgex_r_p1); deallocate( ipgex_r_p2); deallocate( ipgex_r_p3)
00840 deallocate( ipgex_th_p1); deallocate( ipgex_th_p2); deallocate( ipgex_th_p3)
00841
00842 deallocate( emd_p1); deallocate( emd_p2);
00843 deallocate( psif_p1); deallocate( psif_p2);
00844 deallocate(alphf_p1); deallocate(alphf_p2);
00845 deallocate(bvxdf_p1); deallocate(bvxdf_p2);
00846 deallocate(bvydf_p1); deallocate(bvydf_p2);
00847 deallocate(bvzdf_p1); deallocate(bvzdf_p2);
00848 deallocate( rs_p1); deallocate( rs_p2);
00849 deallocate( psi_p1); deallocate( psi_p2); deallocate(psi_p3)
00850 deallocate( alph_p1); deallocate( alph_p2); deallocate(alph_p3)
00851 deallocate( bvxd_p1); deallocate( bvxd_p2); deallocate(bvxd_p3)
00852 deallocate( bvyd_p1); deallocate( bvyd_p2); deallocate(bvyd_p3)
00853 deallocate( bvzd_p1); deallocate( bvzd_p2); deallocate(bvzd_p3)
00854 deallocate( axx_p1); deallocate( axx_p2); deallocate(axx_p3)
00855 deallocate( axy_p1); deallocate( axy_p2); deallocate(axy_p3)
00856 deallocate( axz_p1); deallocate( axz_p2); deallocate(axz_p3)
00857 deallocate( ayy_p1); deallocate( ayy_p2); deallocate(ayy_p3)
00858 deallocate( ayz_p1); deallocate( ayz_p2); deallocate(ayz_p3)
00859 deallocate( azz_p1); deallocate( azz_p2); deallocate(azz_p3)
00860
00861 END PROGRAM coc2cac
