<p><b>dwj07@fsu.edu</b> 2011-09-23 12:09:34 -0600 (Fri, 23 Sep 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Splitting equation_of_state subroutine into functions for each subroutine.<br>
<br>
One step in modularizing equaiton of state.<br>
Next will be separating out each subroutine into it's own module.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-22 23:12:39 UTC (rev 1017)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-23 18:09:34 UTC (rev 1018)
@@ -204,31 +204,30 @@
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND), dimension(:,:), pointer :: rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- integer :: nCells, iCell, k
+ integer :: nCells, iCell, k, indexT, indexS
type (dm_info) :: dminfo
call timer_start("equation_of_state")
if (config_eos_type.eq.'linear') then
- rho => s % rho % array
- tracers => s % tracers % array
- maxLevelCell => grid % maxLevelCell % array
- nCells = grid % nCells
+ call equation_of_state_linear(s, grid, k_displaced, displacement_type)
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- ! Linear equation of state
- rho(k,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(s % index_temperature,k,iCell) &
- + 7.6e-4*tracers(s % index_salinity,k,iCell))
- end do
- end do
-
elseif (config_eos_type.eq.'jm') then
- call equation_of_state_jm(s, grid, k_displaced, displacement_type)
+ tracers => s % tracers % array
+ indexT = s % index_temperature
+ indexS = s % index_salinity
+ if(k_displaced == 0) then
+ rho => s % rho % array
+ else
+ rho => s % rhoDisplaced % array
+ endif
+
+ call equation_of_state_jm(grid, k_displaced, displacement_type, indexT, indexS, tracers, rho)
+! call equation_of_state_jm_bak(s, grid, k_displaced, displacement_type)
+
else
print *, ' Incorrect choice of config_eos_type:',&
config_eos_type
@@ -239,11 +238,59 @@
end subroutine equation_of_state!}}}
- subroutine equation_of_state_jm(s, grid, k_displaced, displacement_type)!{{{
+ subroutine equation_of_state_linear(s, grid, k_displaced, displacement_type)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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, state % rho is returned with no displaced
+ ! If k_displaced>0,the state % rhoDisplaced is returned, and is 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
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+ integer :: k_displaced
+ character(len=8), intent(in) :: displacement_type
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: nCells, iCell, k
+ type (dm_info) :: dminfo
+
+ call timer_start("equation_of_state_linear")
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(s % index_temperature,k,iCell) &
+ + 7.6e-4*tracers(s % index_salinity,k,iCell))
+ end do
+ end do
+
+ call timer_stop("equation_of_state_linear")
+
+ end subroutine equation_of_state_linear!}}}
+
+ subroutine equation_of_state_jm_bak(s, grid, k_displaced, displacement_type)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
! The UNESCO equation of state computed using the
! potential-temperature-based bulk modulus from Jackett and
! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
@@ -371,7 +418,7 @@
integer :: k_test, k_ref
- call timer_start("equation_of_state_jm")
+ call timer_start("equation_of_state_jm_bak")
tracers => s % tracers % array
@@ -415,7 +462,7 @@
if (displacement_type == 'absolute' .and. &
(k_displaced <= 0 .or. k_displaced > nVertLevels) ) then
- write(0,*) 'Abort: In equation_of_state_jm', &
+ write(0,*) 'Abort: In equation_of_state_jm_bak', &
' k_displaced must be between 1 and nVertLevels for ', &
'displacement_type = absolute'
call dmpar_abort(dminfo)
@@ -490,6 +537,256 @@
deallocate(pRefEOS,p,p2)
+ call timer_stop("equation_of_state_jm_bak")
+
+ end subroutine equation_of_state_jm_bak!}}}
+
+ subroutine equation_of_state_jm(grid, k_displaced, displacement_type, indexT, indexS, tracers, rho)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! The UNESCO equation of state computed using the
+ ! potential-temperature-based bulk modulus from Jackett and
+ ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
+ !
+ ! 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
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+ integer :: k_displaced, indexT, indexS
+ character(len=8), intent(in) :: displacement_type
+
+ 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(:,:), intent(inout) :: &
+ rho
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: 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), dimension(:), allocatable :: &
+ 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
+
+ integer :: k_test, k_ref
+
+ call timer_start("equation_of_state_jm")
+
+ 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),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, in-situ density is returned (no displacement)
+ ! If k_displaced/=0, potential density is returned
+
+ ! if displacement_type = 'relative', potential density is calculated
+ ! referenced to level k + k_displaced
+ ! if displacement_type = 'absolute', potential density is calculated
+ ! referenced to level k_displaced for all k
+ ! NOTE: k_displaced = 0 or > nVertLevels is incompatible with 'absolute'
+ ! so abort if necessary
+
+ if (displacement_type == 'absolute' .and. &
+ (k_displaced <= 0 .or. k_displaced > nVertLevels) ) then
+ write(0,*) 'Abort: In equation_of_state_jm', &
+ ' k_displaced must be between 1 and nVertLevels for ', &
+ 'displacement_type = absolute'
+ call dmpar_abort(dminfo)
+ endif
+
+ if (k_displaced == 0) then
+ do k=1,nVertLevels
+ p(k) = pRefEOS(k)
+ p2(k) = p(k)*p(k)
+ enddo
+ else ! k_displaced /= 0
+ do k=1,nVertLevels
+ if (displacement_type == 'relative') then
+ k_test = min(k + k_displaced, nVertLevels)
+ k_ref = max(k_test, 1)
+ else
+ k_test = min(k_displaced, nVertLevels)
+ k_ref = max(k_test, 1)
+ endif
+ p(k) = pRefEOS(k_ref)
+ p2(k) = p(k)*p(k)
+ enddo
+ endif
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ SQ = max(min(tracers(indexS,k,iCell),smax),smin)
+ TQ = max(min(tracers(indexT,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(k) *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2(k)*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p(k))
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p(k) *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2(k)*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p(k))
+
+ rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS,p,p2)
+
call timer_stop("equation_of_state_jm")
end subroutine equation_of_state_jm!}}}
@@ -500,3 +797,4 @@
end module OcnEquationOfState
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
</font>
</pre>