<p><b>mpetersen@lanl.gov</b> 2011-03-10 11:14:37 -0700 (Thu, 10 Mar 2011)</p><p>BRANCH COMMIT: Add adiabatic displacement to Jackett and McDougall EOS.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-09 22:35:35 UTC (rev 754)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-10 18:14:37 UTC (rev 755)
@@ -1636,7 +1636,7 @@
!
! For an isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
- call equation_of_state(s,grid)
+ call equation_of_state(s,grid,-1)
endif
!
@@ -1779,13 +1779,18 @@
end subroutine enforce_boundaryEdge
- subroutine equation_of_state(s, grid)
+ subroutine equation_of_state(s, grid,k_displaced)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This module contains routines necessary for computing the density
! from model temperature and salinity using an equation of state.
!
! Input: grid - grid metadata
! s - state: tracers
+ ! k_displaced
+ ! If k_displaced<=0, density is returned with no displaced
+ ! If k_displaced>0,the density returned is that for a parcel
+ ! adiabatically displaced from its original level to level
+ ! k_displaced. This does not effect the linear EOS.
!
! Output: s - state: computed density
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1793,6 +1798,7 @@
type (state_type), intent(inout) :: s
type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND), dimension(:,:), pointer :: rho
@@ -1819,7 +1825,7 @@
elseif (config_eos_type.eq.'jm') then
- call equation_of_state_jm(s, grid)
+ call equation_of_state_jm(s, grid,k_displaced)
else
print *, ' Incorrect choice of config_eos_type:',&
@@ -1830,7 +1836,7 @@
end subroutine equation_of_state
- subroutine equation_of_state_jm(s, grid)
+ subroutine equation_of_state_jm(s, grid,k_displaced)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This module contains routines necessary for computing the density
! from model temperature and salinity using an equation of state.
@@ -1841,6 +1847,13 @@
!
! Input: grid - grid metadata
! s - state: tracers
+ ! k_displaced
+
+ ! If k_displaced<=0, density is returned with no displaced
+ ! If k_displaced>0,the density returned is that for a parcel
+ ! adiabatically displaced from its original level to level
+ ! k_displaced.
+
!
! Output: s - state: computed density
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1849,6 +1862,7 @@
type (state_type), intent(in) :: s
type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
type (dm_info) :: dminfo
@@ -1880,7 +1894,8 @@
tmin, tmax, &! valid temperature range for level k
smin, smax ! valid salinity range for level k
- real (kind=RKIND) :: p, p2 ! temporary pressure scalars
+ real (kind=RKIND), dimension(:), allocatable :: &
+ p, p2 ! temporary pressure scalars
!-----------------------------------------------------------------------
!
@@ -1975,19 +1990,36 @@
! average temperatures and salinities from Levitus 1994, and
! integrating using hydrostatic balance.
- allocate(pRefEOS(nVertLevels))
+ allocate(pRefEOS(nVertLevels),p(nVertLevels),p2(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
+ ! If k_displaced<=0, density is returned with no displaced
+ ! If k_displaced>0,the density returned is that for a parcel
+ ! adiabatically displaced from its original level to level
+ ! k_displaced.
+ if (k_displaced.le.0) then
+ do k=1,nVertLevels
+ p(k) = pRefEOS(k)
+ p2(k) = p(k)*p(k)
+ enddo
+ elseif (k_displaced.le.nVertLevels) then
+ do k=1,nVertLevels
+ p(k) = pRefEOS(k_displaced)
+ p2(k) = p(k)*p(k)
+ enddo
+ else
+ write(0,*) 'Abort: In equation_of_state_jm', &
+ ' k_displaced must not be larger than nVertLevels'
+ call dmpar_abort(dminfo)
+ endif
+
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
- p = pRefEOS(k)
- p2 = pRefEOS(k)*pRefEOS(k)
-
SQ = max(min(tracers(s%index_salinity,k,iCell),smax),smin)
TQ = max(min(tracers(s%index_temperature,k,iCell),tmax),tmin)
@@ -2012,26 +2044,26 @@
WORK3 = bup0s1t0 + bup0s1t1*TQ + &
(bup0s1t2 + bup0s1t3*TQ)*T2 + &
- p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
- p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ p(k) *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2(k)*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
- bup1sqt0*p)
+ bup1sqt0*p(k))
BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
(bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
- p *(bup1s0t0 + bup1s0t1*TQ + &
+ p(k) *(bup1s0t0 + bup1s0t1*TQ + &
(bup1s0t2 + bup1s0t3*TQ)*T2) + &
- p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ p2(k)*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
SQ*(WORK3 + WORK4)
- DENOMK = 1.0/(BULK_MOD - p)
+ DENOMK = 1.0/(BULK_MOD - p(k))
rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
end do
end do
- deallocate(pRefEOS)
+ deallocate(pRefEOS,p,p2)
end subroutine equation_of_state_jm
</font>
</pre>