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