<p><b>mpetersen@lanl.gov</b> 2010-09-01 17:21:01 -0600 (Wed, 01 Sep 2010)</p><p>Updating topography branch to reflect my local copy. Added Jacket and McDougal EOS and pgrad for diagnostics and testing.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/Registry 2010-09-01 22:54:29 UTC (rev 491)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/Registry 2010-09-01 23:21:01 UTC (rev 492)
@@ -76,9 +76,6 @@
var integer nEdgesOnEdge ( nEdges ) iro nEdgesOnEdge - -
var integer edgesOnCell ( maxEdges nCells ) iro edgesOnCell - -
var integer edgesOnEdge ( maxEdges2 nEdges ) iro edgesOnEdge - -
-var integer maxLevelCell ( nCells ) iro maxLevelCell - -
-var integer maxLevelEdge ( nEdges ) ro maxLevelEdge - -
-var integer maxLevelVertex ( nVertices ) ro maxLevelVertex - -
var real weightsOnEdge ( maxEdges2 nEdges ) iro weightsOnEdge - -
var real dvEdge ( nEdges ) iro dvEdge - -
@@ -105,8 +102,9 @@
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
# Arrays for z-level version of mpas-ocean
-var integer maxLevelsCell ( nCells ) iro kmaxCell - -
-var integer maxLevelsEdge ( nEdges ) iro kmaxEdge - -
+var integer maxLevelCell ( nCells ) iro maxLevelCell - -
+var integer maxLevelEdge ( nEdges ) ro maxLevelEdge - -
+var integer maxLevelVertex ( nVertices ) ro maxLevelVertex - -
var real hZLevel ( nVertLevels ) iro hZLevel - -
var real zMidZLevel ( nVertLevels ) iro zMidZLevel - -
var real zTopZLevel ( nVertLevelsP1 ) iro zTopZLevel - -
@@ -150,6 +148,10 @@
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
var real wTop ( nVertLevelsP1 nCells Time ) o wTop - -
var real ssh ( nCells Time ) o ssh - -
+ ! mrp 100817 for testing
+var real pgrad ( nVertLevels nCells Time ) o pgrad - -
+var real pgrad_edge ( nVertLevels nEdges Time ) o pgrad_edge - -
+ ! mrp 100817 for testing
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-09-01 22:54:29 UTC (rev 491)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-09-01 23:21:01 UTC (rev 492)
@@ -32,12 +32,21 @@
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho
+ real (kind=RKIND), dimension(:), allocatable :: T_temp
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdge, maxLevelVertex
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, boundaryEdge
+ ! from POP, for linear eos
+ real (kind=RKIND), parameter :: &
+ T_leos_ref = 19.0, &! reference T for linear eos (deg C)
+ S_leos_ref = 0.035, &! reference S for linear eos (msu)
+ rho_leos_ref = 1.025022, &! ref dens (g/cm3) at ref T,S and 0 bar
+ alf = 2.55e-4, &! expansion coeff -(drho/dT) (gr/cm^3/K)
+ bet = 7.64e-1 ! expansion coeff (drho/dS) (gr/cm^3/msu)
+
! mrp 100507 end: for diagnostic output
if (config_test_case == 0) then
@@ -201,8 +210,16 @@
zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
enddo
+
+ !!!!!!!!!!!! mrp 100817 setting wind forcing to zero
+ u_src=0.0
+ !!!!!!!!!!!! mrp 100817 setting wind forcing to zero
+
! Set tracers, if not done in grid.nc file
-! tracers = -2.0
+ allocate(t_temp(nVertLevels))
+
+ t_temp = (/ 30.0, 28.0, 26.0, 24.0, 22.0, 20.0, 19.0, 18.0, 17.0, 16.0, 15.0, 14.0, 13.0, 12.0, 11.0, 10.0, 9.0, 8.5, 8.0, 7.5, 7.0, 6.5, 6.0, 5.5, 5.0 /)
+ !tracers = -1.0e34
do iCell = 1,nCells
! tracers(index_tracer2,1,iCell) = maxLevelCell(iCell)
! tracers(index_tracer1,1,iCell) = latCell(iCell)
@@ -212,21 +229,26 @@
! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
- ! for x3, 25 layer test
-! tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
-! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
+ ! for x3, 25 layer test with T constant
+ !tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
+ !tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-! tracers(index_tracer1,iLevel,iCell) = 1.0
+ ! for x3, 25 layer test with S constant
+ !tracers(index_temperature,iLevel,iCell) = t_temp(iLevel) ! temperature
+ !tracers(index_salinity,iLevel,iCell) = 35.0 ! salinity
+
+ tracers(index_tracer1,iLevel,iCell) = 1.0
! tracers(index_tracer1,iLevel,iCell) = maxLevelCell(iCell)
! tracers(index_tracer2,iLevel,iCell) = &
! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
- rho(iLevel,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(index_temperature,iLevel,iCell) &
- + 7.6e-4*tracers(index_salinity,iLevel,iCell))
+ rho(iLevel,iCell) = 1000*(rho_leos_ref + &
+ bet*(tracers(index_salinity,iLevel,iCell)/1000 - S_leos_ref) - &
+ alf*(tracers(index_temperature,iLevel,iCell) - T_leos_ref))
enddo
enddo
+ deallocate(t_temp)
endif
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-09-01 22:54:29 UTC (rev 491)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-09-01 23:21:01 UTC (rev 492)
@@ -115,6 +115,8 @@
! BEGIN RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do rk_step = 1, 4
+! mrp print
+!print '(a,i2)', 'begin rk_step=',rk_step
! --- update halos for diagnostic variables
@@ -265,6 +267,10 @@
divergence, MontPot, pZLevel, zMidEdge, zTopEdge
type (dm_info) :: dminfo
+ ! mrp 100817 for testing
+ real (kind=RKIND), dimension(25) :: maxp, minp, meanp, tempp
+ ! mrp 100817 for testing
+
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdge, maxLevelVertex
integer, dimension(:,:), pointer :: &
@@ -285,7 +291,7 @@
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
- ke_edge => s % ke_edge % array
+ ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
pZLevel => s % pZLevel % array
@@ -404,6 +410,7 @@
deallocate(w_dudzTopEdge)
endif
+
!
! velocity tendency: pressure gradient
!
@@ -429,6 +436,7 @@
enddo
endif
+
!
! velocity tendency: -q(h u^\perp) - </font>
<font color="black">abla K
! +</font>
<font color="black">u_h(</font>
<font color="black">abla \delta + {\bf k}\times </font>
<font color="gray">abla \xi)
@@ -482,9 +490,11 @@
! bottom drag is the same as POP:
! -c |u| u where c is unitless and 1.0e-3.
! see POP Reference guide, section 3.4.4.
+
+! mrp 100729 added /h_edge(k,iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- 1.0e-3*u(k,iEdge) &
- *sqrt(2.0*ke_edge(k,iEdge))
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
! mrp 100603 The following method is more straight forward,
! that the above method of computing ke_edge, but I have
@@ -841,7 +851,7 @@
type (dm_info) :: dminfo
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, rho0Inv
integer :: nCells, nEdges, nVertices, nVertLevels
@@ -851,7 +861,7 @@
hZLevel, ssh, zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
- circulation, vorticity, ke, ke_edge, &
+ circulation, vorticity, ke, ke_edge, pgrad, pgrad_edge, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
zMid, zTop, zMidEdge, zTopEdge, p, pTop, MontPot, rho, temperature, salinity, pZLevel
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -863,7 +873,15 @@
maxLevelCell, maxLevelEdge, maxLevelVertex
real (kind=RKIND) :: r, h1, h2
+ ! from POP, for linear eos
+ real (kind=RKIND), parameter :: &
+ T_leos_ref = 19.0, &! reference T for linear eos (deg C)
+ S_leos_ref = 0.035, &! reference S for linear eos (msu)
+ rho_leos_ref = 1.025022, &! ref dens (g/cm3) at ref T,S and 0 bar
+ alf = 2.55e-4, &! expansion coeff -(drho/dT) (gr/cm^3/K)
+ bet = 7.64e-1 ! expansion coeff (drho/dS) (gr/cm^3/msu)
+
h => s % h % array
u => s % u % array
v => s % v % array
@@ -890,6 +908,9 @@
pTop => s % pTop % array
MontPot => s % MontPot % array
ssh => s % ssh % array
+! mrp
+ pgrad => s % pgrad % array
+ pgrad_edge => s % pgrad_edge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1109,11 +1130,13 @@
!
! Modify PV edge with upstream bias.
!
+!mrp 100720 turn off enstrophy dissipation
do iEdge = 1,nEdges
do k = 1,maxLevelEdge(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
enddo
enddo
+!mrp 100720 turn off enstrophy dissipation end
!
@@ -1149,11 +1172,13 @@
! Modify PV edge with upstream bias.
!
+!mrp 100720 turn off enstrophy dissipation
do iEdge = 1,nEdges
do k = 1,maxLevelEdge(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
enddo
enddo
+!mrp 100720 turn off enstrophy dissipation end
!
! Compute sea surface height
@@ -1167,15 +1192,31 @@
!
! For a isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
+ !call equation_of_state(s,grid)
+
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
! Linear equation of state, for the time being
- rho(k,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(index_temperature,k,iCell) &
- + 7.6e-4*tracers(index_salinity,k,iCell))
+ ! mrp before 100723
+ !rho(k,iCell) = 1000.0*( 1.0 &
+ ! - 2.5e-4*tracers(index_temperature,k,iCell) &
+ ! + 7.6e-4*tracers(index_salinity,k,iCell))
+ ! mrp now directly from POP
+ ! POP:
+! real (r8), parameter :: &
+! T_leos_ref = 19.0_r8, &! reference T for linear eos (deg C)
+! S_leos_ref = 0.035_r8, &! reference S for linear eos (msu)
+! rho_leos_ref = 1.025022_r8, &! ref dens (g/cm3) at ref T,S and 0 bar
+! alf = 2.55e-4_r8, &! expansion coeff -(drho/dT) (gr/cm^3/K)
+! bet = 7.64e-1_r8 ! expansion coeff (drho/dS) (gr/cm^3/msu)
+! if (present(RHOFULL)) RHOFULL = rho_leos_ref + &
+ ! bet*(SALTK - S_leos_ref) - &
+ ! alf*(TEMPK - T_leos_ref)
+ rho(k,iCell) = 1000*(rho_leos_ref + &
+ bet*(tracers(index_salinity,k,iCell)/1000 - S_leos_ref) - &
+ alf*(tracers(index_temperature,k,iCell) - T_leos_ref))
end do
-
end do
endif
@@ -1282,18 +1323,21 @@
do iCell=1,nCells
! compute pressure for z-level coordinates
! assume atmospheric pressure at the surface is zero for now.
+
pZLevel(1,iCell) = rho(1,iCell)*gravity &
* (h(1,iCell)-0.5*hZLevel(1))
+
do k=2,maxLevelCell(iCell)
- delta_p = rho(k,iCell)*gravity*hZLevel(k)
- pZLevel(k,iCell) = pZLevel(k-1,iCell) + 0.5*delta_p
+ pZLevel(k,iCell) = pZLevel(k-1,iCell) &
+ + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ + rho(k ,iCell)*hZLevel(k ))
end do
end do
endif
+
-
! compute vertical velocity
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
@@ -1344,6 +1388,35 @@
endif
+ ! mrp 100817 for testing
+ !
+ ! Compute pressure gradient in each cell
+ !
+ rho0Inv = 1.0/config_rho0
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdge(iEdge)
+ pgrad_edge(k,iEdge) = - rho0Inv*( pZLevel(k,cell2) &
+ - pZLevel(k,cell1) )/dcEdge(iEdge)
+ end do
+ end do
+
+ pgrad(:,:) = 0.0
+ do iCell=1,nCells
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ do k=1,maxLevelEdge(iEdge)
+ pgrad(k,iCell) = pgrad(k,iCell) + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * abs(pgrad_edge(k,iEdge))
+ end do
+ end do
+ do k=1,maxLevelCell(iCell)
+ pgrad(k,iCell) = pgrad(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+ ! mrp 100817 for testing end
+
end subroutine compute_solve_diagnostics
@@ -1389,4 +1462,202 @@
end subroutine enforce_boundaryEdge
+ subroutine equation_of_state(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(inout) :: s
+ type (grid_meta), intent(in) :: grid
+
+
+ type (dm_info) :: dminfo
+ integer :: iEdge, iCell, iVertex, k
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: &
+ TQ,SQ, &! adjusted T,S
+ BULK_MOD, &! Bulk modulus
+ RHO_S, &! density at the surface
+ DRDT0, &! d(density)/d(temperature), for surface
+ DRDS0, &! d(density)/d(salinity ), for surface
+ DKDT, &! d(bulk modulus)/d(pot. temp.)
+ DKDS, &! d(bulk modulus)/d(salinity )
+ SQR,DENOMK, &! work arrays
+ WORK1, WORK2, WORK3, WORK4, T2, depth
+
+ real (kind=RKIND) :: &
+ tmin, tmax, &! valid temperature range for level k
+ smin, smax ! valid salinity range for level k
+
+ real (kind=RKIND) :: p, p2 ! temporary pressure scalars
+
+!-----------------------------------------------------------------------
+!
+! UNESCO EOS constants and JMcD bulk modulus constants
+!
+!-----------------------------------------------------------------------
+
+ !*** for density of fresh water (standard UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ unt0 = 999.842594, &
+ unt1 = 6.793952e-2, &
+ unt2 = -9.095290e-3, &
+ unt3 = 1.001685e-4, &
+ unt4 = -1.120083e-6, &
+ unt5 = 6.536332e-9
+
+ !*** for dependence of surface density on salinity (UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ uns1t0 = 0.824493 , &
+ uns1t1 = -4.0899e-3, &
+ uns1t2 = 7.6438e-5, &
+ uns1t3 = -8.2467e-7, &
+ uns1t4 = 5.3875e-9, &
+ unsqt0 = -5.72466e-3, &
+ unsqt1 = 1.0227e-4, &
+ unsqt2 = -1.6546e-6, &
+ uns2t0 = 4.8314e-4
+
+ !*** from Table A1 of Jackett and McDougall
+
+ real (kind=RKIND), parameter :: &
+ bup0s0t0 = 1.965933e+4, &
+ bup0s0t1 = 1.444304e+2, &
+ bup0s0t2 = -1.706103 , &
+ bup0s0t3 = 9.648704e-3, &
+ bup0s0t4 = -4.190253e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0s1t0 = 5.284855e+1, &
+ bup0s1t1 = -3.101089e-1, &
+ bup0s1t2 = 6.283263e-3, &
+ bup0s1t3 = -5.084188e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0sqt0 = 3.886640e-1, &
+ bup0sqt1 = 9.085835e-3, &
+ bup0sqt2 = -4.619924e-4
+
+ real (kind=RKIND), parameter :: &
+ bup1s0t0 = 3.186519 , &
+ bup1s0t1 = 2.212276e-2, &
+ bup1s0t2 = -2.984642e-4, &
+ bup1s0t3 = 1.956415e-6
+
+ real (kind=RKIND), parameter :: &
+ bup1s1t0 = 6.704388e-3, &
+ bup1s1t1 = -1.847318e-4, &
+ bup1s1t2 = 2.059331e-7, &
+ bup1sqt0 = 1.480266e-4
+
+ real (kind=RKIND), parameter :: &
+ bup2s0t0 = 2.102898e-4, &
+ bup2s0t1 = -1.202016e-5, &
+ bup2s0t2 = 1.394680e-7, &
+ bup2s1t0 = -2.040237e-6, &
+ bup2s1t1 = 6.128773e-8, &
+ bup2s1t2 = 6.207323e-10
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nVertLevels = grid % nVertLevels
+ zMidZLevel => grid % zMidZLevel % array
+
+
+! Jackett and McDougall
+ tmin = -2.0 ! valid pot. temp. range
+ tmax = 40.0
+ smin = 0.0 ! valid salinity, in psu
+ smax = 42.0
+
+ ! This could be put in a startup routine.
+ ! Note I am using zMidZlevel, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters.
+
+! This function computes pressure in bars from depth in meters
+! using a mean density derived from depth-dependent global
+! average temperatures and salinities from Levitus 1994, and
+! integrating using hydrostatic balance.
+
+ allocate(pRefEOS(nVertLevels))
+ do k = 1,nVertLevels
+ depth = -zMidZLevel(k)
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ enddo
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ p = pRefEOS(k)
+ p2 = pRefEOS(k)*pRefEOS(k)
+
+ SQ = max(min(tracers(index_salinity,k,iCell),smax),smin)
+ TQ = max(min(tracers(index_temperature,k,iCell),tmax),tmin)
+
+ SQR = sqrt(SQ)
+ T2 = TQ*TQ
+
+ !***
+ !*** first calculate surface (p=0) values from UNESCO eqns.
+ !***
+
+ WORK1 = uns1t0 + uns1t1*TQ + &
+ (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
+ WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)
+
+ RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
+ + (uns2t0*SQ + WORK1 + WORK2)*SQ
+
+ !***
+ !*** now calculate bulk modulus at pressure p from
+ !*** Jackett and McDougall formula
+ !***
+
+ WORK3 = bup0s1t0 + bup0s1t1*TQ + &
+ (bup0s1t2 + bup0s1t3*TQ)*T2 + &
+ p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p)
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p)
+
+ rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS)
+ end subroutine equation_of_state
+
end module time_integration
</font>
</pre>