module grid_structure use constants !## SPECIFICATION ################################################################################################################## implicit none ! Path of the directory where the grid is located character(len = *), parameter :: grid_dir = "edit this string with your path to the directory with the Zref folders" ! The following block define the structure of the grid at a directory-tree level integer, parameter :: ngZ = 8, & ! number of grid nodes: reference metallicity ngX = 3, & ! number of grid nodes: hydrogen mass fraction ngCO = 10, & ! number of grid nodes: carbon-to-oxygen ratio (actually fCO) ngC = 10, & ! number of grid nodes: carbon abundance (actually fC) ngN = 1, & ! number of grid nodes: nitrogen abundance (actually fN) ngML = 3, & ! number of grid nodes: mixing length parameter ngM = 44, & ! number of grid nodes: total mass ngMc = 2 ! number of grid nodes: core mass (actually core mass - luminosity functions) ! Values of the grid nodes for each parameter of the grid real(8), parameter, dimension(ngZ) :: gZ = (/1.D-3, 2.D-3, 4.D-3, 6.D-3, 8.D-3, 1.D-2, 1.4D-2, 1.7D-2/) real(8), parameter, dimension(ngX) :: gX = (/6.D-1, 7.D-1, 8.D-1/) real(8), parameter, dimension(ngCO) :: gCO = (/-2.63D-1, 0.D0, 1.05D-1, 2.38D-1, 2.6D-1, & 2.81D-1, 3.74D-1, 5.15D-1, 7.37D-1, 9.59D-1/) real(8), parameter, dimension(ngC) :: gC = (/-2.63D-1, 0.D0, 1.05D-1, 2.38D-1, 2.6D-1, & 2.81D-1, 3.74D-1, 5.15D-1, 7.37D-1, 9.59D-1/) real(8), parameter, dimension(ngN) :: gN = (/0.D0/) real(8), parameter, dimension(ngML) :: gML = (/1.5D0, 2.D0, 2.5D0/) real(8), parameter, dimension(ngM) :: gM = (/6.D-1, 6.5D-1, 7.D-1, 7.5D-1, 8.D-1, 8.5D-1, 9.D-1, 9.5D-1, 1.D0, & 1.1D0, 1.2D0, 1.3D0, 1.4D0, 1.5D0, 1.6D0, 1.7D0, 1.8D0, 1.9D0, 2.D0, & 2.2D0, 2.4D0, 2.6D0, 2.8D0, 3.D0, 3.2D0, 3.4D0, 3.6D0, 3.8D0, 4.D0, & 4.2D0, 4.4D0, 4.6D0, 4.8D0, 5.D0, 5.2D0, 5.4D0, 5.6D0, 5.8D0, 6.D0, & 6.2D0, 6.4D0, 6.6D0, 6.8D0, 7.D0/) character(len = 5), parameter, dimension(ngMc) :: gMc = (/'lower', 'upper'/) ! the following block describe the structure of the grid at a single table-file level integer, parameter :: jlogL = 1, & ! index of log(L/Lsun) column in table files jlogTe = 2, & ! index of log(Teff) column in table files jM = 3, & ! index of M/Msun column in table files jMc = 4, & ! index of Mc/Msun column in table files jX = 5, & ! index of X column in table files jZ = 6, & ! index of Z column in table files jXC = 7, & ! index of XC column in table files jXN = 8, & ! index of XN column in table files jXO = 9, & ! index of XO column in table files joa = 1, & ! index of adiabatic frequency column among the three columns for pulsation properties jor = 2, & ! index of real part of the frequency column among the three columns for pulsation properties joi = 3 ! index of imag part of the frequency column among the three columns for pulsation properties integer, parameter :: npar = 9, & ! number of columns in the tables not related to pulsation modes (logL, logTe, M, ...) npul = 3, & ! number of columns with pulsation properties for each mode nmodes = 5, & ! number of modes per pulsation model ntot = npar + npul * nmodes ! the following block is for other variables used during execution real(8), parameter :: UNDEF = -9.99D0 ! "undefined" value ! Variable arrays to store selected interpolation nodes for each case, initialized to UNDEF real(8), dimension(2) :: snZ = (/UNDEF, UNDEF/), & sNX = (/UNDEF, UNDEF/), & snCO = (/UNDEF, UNDEF/), & snC = (/UNDEF, UNDEF/), & snN = (/UNDEF, UNDEF/), & snM = (/UNDEF, UNDEF/), & snMc = (/UNDEF, UNDEF/) real(8), dimension(3) :: snML = (/UNDEF, UNDEF, UNDEF/) integer :: conversion = 0 !## SUB-PROGRAMS ################################################################################################################### contains !=================================================================================================================================== ! Subroutine: set_interp_nodes ! Arguments: ! (input) ! value: real ! lg: integer ! grid: array of real ! (output) ! nodes: array of 2 real ! Description: ! Searches the sorted array "grid" of length "lg" for 2 contiguous values "nodes(1)" and "nodes(2)" ! such that nodes(1) <= value <= nodes(2), given the argument "value", in order to perform linear interpolation in such nodes; ! if value <= min(grid), then nodes(1) and nodes(2) are set as the two smallest values in grid, for extrapolation downwards ! if value >= max(grid), then nodes(1) and nodes(2) are set as the two largest values in grid, for extrapolation upwards subroutine set_interp_nodes(value, grid, lg, nodes, ret) real(8), intent(in) :: value ! value to interpolate for integer, intent(in) :: lg ! total number of nodes for the selected parameter real(8), dimension(lg), intent(in) :: grid ! values of each node for the selected parameter real(8), dimension(2), intent(out) :: nodes ! the 2 chosen interpolation (possibly extrapolation) nodes integer, intent(out) :: ret integer :: ig ! running index 10 format('Something went wrong while searching for interpolation nodes...') ! if the value to interpolate for is smaller than all nodes, ! the two smallest nodes are returned for downwards extrapolation if (value < grid(1)) then nodes(1) = grid(1) nodes(2) = grid(2) ret = 1 return ! if the value to interpolate for is larger than all nodes, ! the two largest nodes are returned for upwards extrapolation else if (value > grid(lg)) then nodes(1) = grid(lg - 1) nodes(2) = grid(lg) ret = 2 return endif ! scan all nodes until two are found such that node1 <= value <= node2 do ig = 2, lg if (value <= grid(ig)) then nodes(1) = grid(ig - 1) nodes(2) = grid(ig) ret = 0 return endif enddo ! There should be no error while searching for the interpolation nodes, ! however in such a case both nodes are set to UNDEF nodes(1) = UNDEF nodes(2) = UNDEF ret = 9 write(*,10) return end subroutine set_interp_nodes !=================================================================================================================================== ! Subroutine: interpolate_table_luminosity ! Arguments: ! (input) ! logL: real ! input_name: charcater ! (output) ! table_data: array of real ! err: integer ! Description: ! Reads a table file from a grid of pulsation models. A table file represents a sequence of models with increasing luminosity. ! The file is specified by its path through the input argument "input_name". The file is read line by line until two lines are ! found which luminosity values bound the input argument "logL", meaning that logL1 <= logL <= logL2, where logL1 and logL2 are ! the luminosity values from the two lines. Then the subroutine interpolates linearly in logL between logL1 and logL2 ! for each parameter in the columns of the grid. Note that such parameters are defined by the variables in the specification ! part of the present module. Interpolated data are returned via the output argument "table_data". A normal execution returns ! also the output argument "err" with value 0. If err is 1, it means the execution was ok, but logL was larger than the maximum ! luminosity in the file, so the subroutine extrapolated linearly upwards. If err is -1, it extrapolated downwards. ! If err has another value, some error occurred: ! err = 2: file not found ! err = 3: there was no data in the file (the sequence was "empty") ! err = 4: the file contained only one line (one model), neither interpolation or extrapolation are possible ! err = 9: some other unidentified error subroutine interpolate_table_luminosity(logL, input_name, conversion, table_data, err) implicit none real(8), intent(in) :: logL ! luminosity value to interpolate for character(len = *), intent(in) :: input_name ! table file name integer, intent(in) :: conversion ! specify how to convert the frequency data, recommendend = 2: ! conversion = 0: use omega_ad[1/rad], omega_R[1/rad], omega_I[1/rad] ! conversion = 1: use P_ad[d], P[d], GR ! conversion = 2: use log(P_ad[d]), log(P[d]), GR real(8), dimension(ntot), intent(out) :: table_data ! global and pulsation parameters after interpolation in logL integer, intent(out) :: err ! error variable: ! err = 0: subroutine worked fine ! err = -1: subroutine worked fine, but extrapolated downward ! err = +1: subroutine worked fine, but extrapolated upward ! err = 2: file not found error ! err = 3: empty sequence error ! err = 4: single model error ! err = 9: unknown error integer :: lTN, nmodes_in, i, ja, jr, ji, ierr real(8), dimension(ntot) :: row0, row1, row2 real(8) :: dx, w1, w2 ! linear interpolation quantities: ! interpolating linearly in x between nodes x1 and x2 with values y1 and y2 ! y(x) = w1 * y1 + w2 * y2 ! where weights w1, w2 are ! w1 = (x2 - x) / dx ! w2 = (x - x1) / dx ! and dx = (x2 - x1) != 0 real(8) :: oA, oR, oI character(len = 1000) :: line logical :: input_exist !----------------------------------------------------------------------------------------------------------------------------------- 10 format('File not found: ' ,a) 11 format('Error while reading file: ',a) 12 format('Sequence is empty in file: ',a) 13 format('Sequence has only 1 model in file: ',a) 14 format('Unknown conversion code: ',i5) 15 format('Mismatch in number of nodes: expected: ',i5,', got: 'i5) !----------------------------------------------------------------------------------------------------------------------------------- ! initialize error value to 9 = unknown error err = 9 ! get length of input_name lTN = len_trim(input_name) ! check that table file exists, otherwise send stdout message, set error value to 2 and return inquire(file = input_name(1:lTN), exist = input_exist) if (.not.input_exist) then write(*,10) input_name(1:lTN) err = 2 return endif ! open table file via unit 42 open(unit = 42, file = input_name(1:lTN), form = 'formatted', status = 'old') ! read first line of table file read(42, '(a500)', iostat = ierr) line ! check for errors, if so set error value to 9, send stdout message and return if (ierr /= 0) then write(*,11) input_name(1:lTN) err = 9 return endif ! check if line is comment (means there is a legend header), ! if so read the second line and check for errors as did for the first line if (line(1:1) == '#') then read(42, '(a500)', iostat = ierr) line if (ierr /= 0) then write(*,11) input_name(1:lTN) err = 9 return endif endif ! read number of modes from first (possibly second) line of the table file read(line, '(i4)') nmodes_in ! check that number of modes is consistent with what expected from the defined grid structure if (nmodes_in /= nmodes) then write(*,15) nmodes, nmodes_in return endif ! read first row of the table and check for reading errors, ! in which case send stdout message, set error value to 3 and return read(42, *, iostat = ierr) (row1(i), i = 1, ntot) if (ierr /= 0) then write(*,12) input_name(1:lTN) err = 3 return endif ! read second row of the table and check for reading errors, ! in which case send stdout message, set error value to 4 and return read(42, *, iostat = ierr) (row2(i), i = 1, ntot) if (ierr /= 0) then write(*,13) input_name(1:lTN) err = 4 return endif ! read the rest of the table until two rows with luminosity values suitable for interpolation are found do if (logL <= row2(jlogL)) exit read(42, *, iostat = ierr) (row0(i), i = 1, ntot) if (ierr < 0) then exit else if (ierr == 0) then row1(:) = row2(:) row2(:) = row0(:) else write(*,11) err = 9 return endif enddo if (conversion == 1) then ! convert frequencies for all modes into P_ad[days], P[days] and GR do i = 1, nmodes ja = npar + npul * (i - 1) + joa jr = npar + npul * (i - 1) + jor ji = npar + npul * (i - 1) + joi oA = row1(ja) oR = row1(jr) oI = row1(ji) row1(ja) = TWOPI / oA / SECONDS_PER_DAY row1(jr) = EXP(TWOPI * oR / oI) - 1.D0 row1(ji) = TWOPI / oI / SECONDS_PER_DAY oA = row2(ja) oR = row2(jr) oI = row2(ji) row2(ja) = TWOPI / oA / SECONDS_PER_DAY row2(jr) = EXP(TWOPI * oR / oI) - 1.D0 row2(ji) = TWOPI / oI / SECONDS_PER_DAY enddo else if (conversion == 2) then ! convert frequencies for all modes into log(P_ad[days]), log(P[days]) and GR do i = 1, nmodes ja = npar + npul * (i - 1) + joa jr = npar + npul * (i - 1) + jor ji = npar + npul * (i - 1) + joi oA = row1(ja) oR = row1(jr) oI = row1(ji) row1(ja) = LOG10(TWOPI / oA / SECONDS_PER_DAY) row1(jr) = EXP(TWOPI * oR / oI) - 1.D0 row1(ji) = LOG10(TWOPI / oI / SECONDS_PER_DAY) oA = row2(ja) oR = row2(jr) oI = row2(ji) row2(ja) = LOG10(TWOPI / oA / SECONDS_PER_DAY) row2(jr) = EXP(TWOPI * oR / oI) - 1.D0 row2(ji) = LOG10(TWOPI / oI / SECONDS_PER_DAY) enddo else if (conversion /= 0) then ! in case of invalid conversion code, send stdout message and return write(*,14) conversion err = 9 return endif ! compute linear interpolation weights dx = row2(jlogL) - row1(jlogL) w1 = (row2(jlogL) - logL) / dx w2 = (logL - row1(jlogL)) / dx ! interpolate all parameters and store temporarily in row0 row0(jlogL) = logL if (conversion == 2) then do i = 2, npar row0(i) = w1 * row1(i) + w2 * row2(i) enddo do i = 1, nmodes ja = npar + npul * (i - 1) + joa jr = npar + npul * (i - 1) + jor ji = npar + npul * (i - 1) + joi if ((ABS(row1(jr) + 1.D0) <= 1.D-16).or.(ABS(row2(jr) + 1.D0) <= 1.D-16)) then row0(ja) = 1.D-10 row0(jr) = -1.D0 row0(ji) = 1.D-10 else row0(ja) = w1 * row1(ja) + w2 * row2(ja) row0(jr) = w1 * row1(jr) + w2 * row2(jr) row0(ji) = w1 * row1(ji) + w2 * row2(ji) endif enddo else do i = 2, ntot row0(i) = w1 * row1(i) + w2 * row2(i) enddo endif ! copy row0 onto table_data table_data(:) = row0(:) ! set successful error value, according to whether data were interpolated or extrapolated if (logL < row1(jlogL)) then err = -1 ! downward extrapolation else if (logL > row2(jlogL)) then err = 1 ! upward extrapolation else err = 0 ! interpolation endif end subroutine end module grid_structure