Write the parameter file to specify the \(\mathrm{EOS}\) when computing \(\mathrm{BNS}\) or \(\mathrm{DRS}\) with \(\texttt{LORENE}\)
FT
| Type | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|
| character(len=4) | :: | eos | ||||
| character(len=4) | :: | eos_type | ||||
| character(len=:), | ALLOCATABLE | :: | err_msg | |||
| logical | :: | exist | ||||
| double precision | :: | gamma_poly | ||||
| integer | :: | ios | ||||
| double precision | :: | kappa0_lorene | ||||
| double precision | :: | kappa_poly | ||||
| double precision | :: | log10_p0_lorene | ||||
| double precision | :: | log10_rho0_lorene | ||||
| double precision | :: | log10_rho1_lorene | ||||
| double precision | :: | log10_rho2_lorene | ||||
| character(len=:), | ALLOCATABLE | :: | namefile | |||
| integer, | parameter | :: | npoly | = | 4 |
PROGRAM write_par_eos
!*****************************************************
!
!# Write the \(\mathrm{par\_eos.d}\) parameter file
! to specify the |eos| when computing |bns| or |drs|
! with |lorene|
!
! FT
!
!*****************************************************
USE utility, ONLY: density_si2cu, k_lorene2cu, k_lorene2cu_pwp, one
USE constants, ONLY: kg2g, m2cm
USE pwp_EOS, ONLY: Gamma0, K0, select_EOS_parameters, gen_pwp_eos, &
get_rho_0, get_rho_1, get_rho_2, &
get_Gamma1, get_Gamma2, get_Gamma3, &
get_K1, get_K2, get_K3
IMPLICIT NONE
INTEGER, PARAMETER:: npoly= 4
INTEGER:: ios
!DOUBLE PRECISION:: gamma1_lorene, gamma2_lorene, gamma3_lorene
DOUBLE PRECISION:: kappa0_lorene, log10_p0_lorene, log10_rho0_lorene, &
log10_rho1_lorene, log10_rho2_lorene
DOUBLE PRECISION:: kappa_poly, gamma_poly
LOGICAL:: exist
CHARACTER(LEN= :), ALLOCATABLE:: namefile
CHARACTER(LEN= :), ALLOCATABLE:: err_msg
CHARACTER(LEN= 4):: eos_type
CHARACTER(LEN= 4):: eos
namefile= "par_eos.d"
WRITE(*,'(A,/,A)') "** Do you want to produce a parameter file for a polytropic EOS or piecewise polytropic EOS?", " Please type `poly` for the first, `pwp` for the second: "
READ(*,'(A)') eos_type
select_eos_type: SELECT CASE( eos_type )
CASE( "poly" )
#ifdef __INTEL_COMPILER
WRITE(*,'("** Please write the DOUBLE PRECISION value of the polytropic exponent gamma: ",\)')
READ(*,'(F)') gamma_poly
#endif
#ifdef __GFORTRAN__
PRINT *, "** Please write the DOUBLE PRECISION value of the polytropic constant gamma, with 4 total digits and 3 decimal digits: "
READ(*,'(F4.3)') gamma_poly
#endif
#ifdef __INTEL_COMPILER
WRITE(*,'("** Please write the DOUBLE PRECISION value of the polytropic constant K in SPHINCS units: ",\)')
READ(*,'(F)') kappa_poly
#endif
#ifdef __GFORTRAN__
PRINT *, "** Please write the DOUBLE PRECISION value of the polytropic constant K in SPHINCS units, with 6 total digits and 2 decimal digits: "
READ(*,'(F6.2)') kappa_poly
#endif
kappa_poly= kappa_poly/k_lorene2cu(gamma_poly)
INQUIRE( FILE= TRIM(namefile), EXIST= exist )
IF( exist )THEN
OPEN( UNIT= 2, FILE= TRIM(namefile), STATUS= "REPLACE", &
FORM= "FORMATTED", &
POSITION= "REWIND", ACTION= "WRITE", IOSTAT= ios, &
IOMSG= err_msg )
ELSE
OPEN( UNIT= 2, FILE= TRIM(namefile), STATUS= "NEW", &
FORM= "FORMATTED", &
ACTION= "WRITE", IOSTAT= ios, IOMSG= err_msg )
ENDIF
IF( ios > 0 )THEN
PRINT *, "...error when opening " // TRIM(namefile), &
". The error message is", err_msg
STOP
ENDIF
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(A67)" ) &
"1 Polytropic EOS (cf. documentation of Eos::eos_from_file)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(A14)" ) &
"Polytropic EOS"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F11.9,A18)" ) &
gamma_poly, "gamma"!, polytropic exponent"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F19.17,A38)" ) &
kappa_poly, "kappa [rho_nuc c^2 / n_nuc^gamma]"!, polytropic constant"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F11.9,A37)" ) &
one, "mean particle mass [m_b]"
CLOSE( UNIT= 2 )
CASE( "pwp " )
#ifdef __INTEL_COMPILER
WRITE(*,'("** Please write an up-to-4-character string containing the name " &
"of the piecewise polytropic EOS: ",\)')
#endif
#ifdef __GFORTRAN__
WRITE(*,'("** Please write an up-to-4-character string containing the name " &
"of the piecewise polytropic EOS: ")')
#endif
READ(*,'(A)') eos
CALL select_EOS_parameters(eos)
kappa0_lorene= K0/k_lorene2cu_pwp(Gamma0)
log10_p0_lorene= LOG10( K0/k_lorene2cu_pwp(Gamma0) &
*( get_rho_0()/density_si2cu*kg2g/(m2cm**3.0D0) )**(Gamma0) )
log10_rho0_lorene= LOG10( get_rho_0()/density_si2cu*kg2g/(m2cm**3.0D0) )
log10_rho1_lorene= LOG10( get_rho_1()/density_si2cu*kg2g/(m2cm**3.0D0) )
log10_rho2_lorene= LOG10( get_rho_2()/density_si2cu*kg2g/(m2cm**3.0D0) )
INQUIRE( FILE= TRIM(namefile), EXIST= exist )
IF( exist )THEN
OPEN( UNIT= 2, FILE= TRIM(namefile), STATUS= "REPLACE", &
FORM= "FORMATTED", &
POSITION= "REWIND", ACTION= "WRITE", IOSTAT= ios, &
IOMSG= err_msg )
ELSE
OPEN( UNIT= 2, FILE= TRIM(namefile), STATUS= "NEW", &
FORM= "FORMATTED", &
ACTION= "WRITE", IOSTAT= ios, IOMSG= err_msg )
ENDIF
IF( ios > 0 )THEN
PRINT *, "...error when opening " // TRIM(namefile), &
". The error message is", err_msg
STOP
ENDIF
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(A69)" ) &
"110 Multipolytropic EOS (cf. document. of Eos::eos_multi_poly)"
!WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(A16,A4,A4)" ) eos
! "Multipolytropic ", eos, " EOS"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(A4)" ) &
ADJUSTL(TRIM(eos))
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(I1,A52)" ) &
npoly, "npoly, number of polytropes"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F11.9,A43)" ) &
Gamma0, "gamma0, crust (here from SLy)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F11.9,A55)" ) &
get_Gamma1(), "gamma1, adiabatic indexes (crust to core)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F11.9,A13)" ) &
get_Gamma2(), "gamma2"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F11.9,A13)" ) &
get_Gamma3(), "gamma3"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(ES16.9E3,A55)" ) &
kappa0_lorene, &
"kappa0, K for gamma0 polytrope (here from SLy)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F12.9,A73)" ) &
log10_p0_lorene, &
"log10(P0/c^2), log(p) between gamma0 and gamma1 (dyne/cm^2/c_cgs^2)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F12.9,A58)" ) &
log10_rho0_lorene, &
"log10(rho0), logs of transition densities (g/cm^3)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F12.9,A17)" ) &
log10_rho1_lorene, "log10(rho1)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F12.9,A17)" ) &
log10_rho2_lorene, "log10(rho2)"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F3.1,A77)" ) &
0.0, "decInc, array (size npoly-1) of percentages which change"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F3.1,A69)" ) &
0.0, "the transition densities by changing the"
WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = "(F3.1,A73)" ) &
0.0, "transition enthalpies (set to 0. to disable)"
CLOSE( UNIT= 2 )
END SELECT select_eos_type
PRINT *
PRINT *, "** Parameter file ", namefile," written."
PRINT *
END PROGRAM write_par_eos