PageRenderTime 49ms CodeModel.GetById 28ms app.highlight 20ms RepoModel.GetById 0ms app.codeStats 0ms

/tools/generate_ideal_table.f

https://bitbucket.org/NeilMiller/h5tbuilder
FORTRAN Legacy | 164 lines | 106 code | 35 blank | 23 comment | 0 complexity | f88e17aa8ee57a3059e12358c42886f5 MD5 | raw file
  1!! This simple procedure is necessary to generate the h5 table
  2
  3program generate_water_table
  4  use aneos_access
  5  use ideal_water
  6  implicit none
  7  
  8  integer :: i, j, NRho, NT, outid, outid2, ios
  9  double precision :: logRho, logT, lnRho, lnT, Rho, T, Fvect(num_free_derivs)
 10  double precision, pointer, dimension(:) :: logRhovect, logTvect
 11  double precision, pointer, dimension(:,:) :: Fr, dfdrho, dfdt, d2fdrho2, &
 12       d2fdt2, d2fdrhodt, d3fdrho2dt, d3fdrhodt2, d4fdrho2dt2
 13  double precision :: logRhoLow, logRhoHi, dlogRho
 14  double precision :: logTLow, logTHi, dlogT
 15  character (len=256) :: filename
 16  double precision :: f_val, df_drho, df_dt, d2f_drho2, d2f_dt2, d2f_drhodt, &
 17       d3f_drho3, d3f_drho2dt, d3f_drhodt2, d3f_dt3, d4f_drho2dt2
 18
 19  NRho = 81
 20  NT = 31
 21  
 22  allocate(Fr(NRho, NT))
 23  allocate(dfdrho(NRho, NT))
 24  allocate(dfdt(NRho, NT))
 25  allocate(d2fdrho2(NRho, NT))
 26  allocate(d2fdt2(NRho, NT))
 27  allocate(d2fdrhodt(NRho, NT))
 28  allocate(d3fdrho2dt(NRho, NT))
 29  allocate(d3fdrhodt2(NRho, NT))
 30  allocate(d4fdrho2dt2(NRho, NT))
 31
 32  allocate(logRhovect(NRho))
 33  allocate(logTvect(NT))
 34
 35  logRhoLow = -5.
 36  logRhoHi = 3.
 37  logTLow = 2.8
 38  logTHi = 6.8
 39  
 40  dlogRho = (logRhoHi - logRhoLow) / (NRho - 1)
 41  dlogT = (logTHi - logTLow) / (NT - 1)
 42  
 43  !! Write out both formatted and unformatted data
 44  !! I am trying to shift to the unformatted approach since it
 45  !!  1. takes up less space
 46  !!  2. Doesn't introduce more numerical errors 
 47  outid = 51
 48  outid2 = 52
 49  open(unit=outid,file='data/wateridealfe.bin', form='unformatted')  
 50  open(unit=outid2,file='data/wateridealfe.dat')
 51  
 52  do i=1,NRho
 53     do j=1,NT
 54        logRho = logRhoLow + (i-1) * dlogRho
 55        logT = logTLow + (j-1) * dlogT
 56        logRhovect(i) = logRho
 57        logTvect(j) = logT
 58        
 59        Rho = 10.d0**logRho
 60        T = 10.d0**logT
 61        
 62        lnRho = log(Rho)
 63        lnT = log(T)
 64        
 65        call get_IG_EOS(lnRho, lnT, f_val, df_drho, df_dt, &
 66             d2f_drho2, d2f_dt2, d2f_drhodt, &
 67             d3f_drho3, d3f_drho2dt, d3f_drhodt2, d3f_dt3, &
 68             d4f_drho2dt2)
 69
 70        Fvect(a_val) = f_val
 71        Fvect(a_dRho) = df_drho
 72        Fvect(a_dT) = df_dt
 73        Fvect(a_dRho2) = d2f_drho2
 74        Fvect(a_dT2) = d2f_dt2
 75        Fvect(a_dRhodT) = d2f_drhodt
 76        Fvect(a_dRho2dT) = d3f_drho2dt
 77        Fvect(a_dRhodT2) = d3f_drhodt2
 78        Fvect(a_dRho2dT2) = d4f_drho2dt2
 79        
 80        write(outid) lnRho, lnT, Fvect
 81        write(outid2,*) lnRho, lnT, Fvect
 82
 83        Fr(i,j) = Fvect(a_val)
 84        dfdrho(i,j) = Fvect(a_dRho)
 85        dfdt(i,j) = Fvect(a_dT)
 86        d2fdrho2(i,j) = Fvect(a_dRho2)
 87        d2fdt2(i,j) = Fvect(a_dT2)
 88        d2fdrhodt(i,j) = Fvect(a_dRhodT)
 89        d3fdrho2dt(i,j) = Fvect(a_dRho2dT)
 90        d3fdrhodt2(i,j) = Fvect(a_dRhodT2)
 91        d4fdrho2dt2(i,j) = Fvect(a_dRho2dT2)
 92        
 93     enddo
 94  enddo
 95
 96  filename = "h5tdata/fr_logRhovect.data"
 97  open(unit=19,file=trim(filename),status='replace', &
 98         iostat=ios,action='write',form='formatted')     
 99  write(19,'(e20.6)') logRhovect
100  close(unit=19)
101  
102  filename = "h5tdata/fr_logTvect.data"
103  open(unit=19,file=trim(filename),status='replace', &
104         iostat=ios,action='write',form='formatted')     
105  write(19,'(e20.6)') logTvect
106  close(unit=19)
107
108  filename = "h5tdata/fr.data"
109  call write_matrix(Fr, filename)
110  
111  filename = "h5tdata/dfdrho.data"
112  call write_matrix(dfdrho, filename)
113  
114  filename = "h5tdata/dfdt.data"
115  call write_matrix(dfdt, filename)  
116
117  filename = "h5tdata/d2fdrho2.data"
118  call write_matrix(d2fdrho2, filename)
119  
120  filename = "h5tdata/d2fdt2.data"
121  call write_matrix(d2fdt2, filename)
122  
123  filename = "h5tdata/d2fdRhodT.data"
124  call write_matrix(d2fdrhodt, filename)
125  
126  filename = "h5tdata/d3fdRho2dT.data"
127  call write_matrix(d3fdrho2dt, filename)
128  
129  filename = "h5tdata/d3fdRhodT2.data"
130  call write_matrix(d3fdrhodt2, filename)
131  
132  filename = "h5tdata/d4fdRho2dT2.data"
133  call write_matrix(d4fdrho2dt2, filename)
134  
135  close(unit=outid)
136  close(unit=outid2)
137  deallocate(Fr)
138  deallocate(dfdrho)
139  deallocate(dfdt)
140  deallocate(d2fdrho2)
141  deallocate(d2fdt2)
142  deallocate(d2fdrhodt)
143  deallocate(d3fdrho2dt)
144  deallocate(d3fdrhodt2)
145  deallocate(d4fdrho2dt2)
146  deallocate(logRhovect)
147  deallocate(logTvect)
148
149  
150
151
152contains
153  subroutine write_matrix(matrix, filename)
154    double precision, intent(in), dimension(:,:), pointer :: matrix
155    character (len=256) :: filename
156    
157    open(unit=19,file=trim(filename),status='replace', &
158         iostat=ios,action='write',form='formatted')
159    write(19,'(e20.6)') matrix
160    close(unit=19)
161  end subroutine write_matrix
162
163end program generate_water_table
164