<p><b>dwj07@fsu.edu</b> 2011-09-29 09:06:23 -0600 (Thu, 29 Sep 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Splitting equation of state module into separate modules.<br>
<br>
Also added some folds and vim mode lines to ThickVadv.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-28 23:30:51 UTC (rev 1040)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-29 15:06:23 UTC (rev 1041)
@@ -34,11 +34,13 @@
module_OcnVmixCoefsRich.o \
module_OcnVmixCoefsTanh.o \
module_OcnRestoring.o \
- module_OcnEquationOfState.o \
module_OcnTendency.o \
module_OcnTimeIntegration.o \
module_OcnTimeIntegrationRK4.o \
module_OcnTimeIntegrationSplit.o \
+ module_OcnEquationOfState.o \
+ module_OcnEquationOfStateJM.o \
+ module_OcnEquationOfStateLinear.o \
module_global_diagnostics.o
@@ -123,8 +125,13 @@
module_OcnVmixCoefsTanh.o:
-module OcnEquationOfState.o:
+module_OcnEquationOfState.o: module_OcnEquationOfStateJM.o module_OcnEquationOfStateLinear.o
+module_OcnEquationOfStateJM.o:
+
+module_OcnEquationOfStateLinear.o:
+
+
module_mpas_core.o: module_advection.o \
module_OcnThickHadv.o \
module_OcnThickVadv.o \
@@ -160,6 +167,8 @@
module_OcnVmixCoefsRich.o \
module_OcnVmixCoefsTanh.o \
module_OcnEquationOfState.o \
+ module_OcnEquationOfStateJM.o \
+ module_OcnEquationOfStateLinear.o \
module_OcnTendency.o \
module_OcnTimeIntegration.o \
module_OcnTimeIntegrationRK4.o \
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-28 23:30:51 UTC (rev 1040)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-29 15:06:23 UTC (rev 1041)
@@ -4,7 +4,7 @@
!
!> \brief MPAS ocean equation of state driver
!> \author Doug Jacobsen
-!> \date 19 September 2011
+!> \date 28 September 2011
!> \version SVN:$Id:$
!> \details
!> This module contains the main driver routine for calling
@@ -17,8 +17,8 @@
use grid_types
use configure
use timer
-! use OcnEquationOfStateLinear
-! use OcnEquationOfStateJM
+ use OcnEquationOfStateJM
+ use OcnEquationOfStateLinear
implicit none
private
@@ -37,8 +37,7 @@
!--------------------------------------------------------------------
public :: OcnEquationOfStateRho, &
- OcnEquationOfStateInit, &
- equation_of_state
+ OcnEquationOfStateInit
!--------------------------------------------------------------------
!
@@ -55,131 +54,18 @@
!***********************************************************************
!
-! routine OcnEquationOfState
+! routine OcnEquationOfStateRho
!
!> \brief Calls equation of state
!> \author Doug Jacobsen
-!> \date 19 September 2011
+!> \date 28 September 2011
!> \version SVN:$Id$
!> \details
!> This routine calls the equation of state to update the density
!
!-----------------------------------------------------------------------
- subroutine OcnEquationOfStateRho(grid, displacement_type, k_displaced, indexT, indexS, tracers, rho, err)
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- character(len=8), intent(in) :: displacement_type
-
- real (kind=RKIND), dimension(:,:,:), intent(in) :: &
- tracers !< Input: tracers
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- integer, intent(in) :: indexT, indexS, k_displaced
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(inout) :: &
- rho !< Output: density
-
- integer, intent(out) :: err
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: err1, err2
-
- !-----------------------------------------------------------------
- !
- ! call relevant routines for computing tendencies
- ! note that the user can choose multiple options and the
- ! tendencies will be added together
- !
- !-----------------------------------------------------------------
-
- err = 0
- if(.not.eosON) return
-
-! call OcnEquationOfStateLinearRho(grid, indexT, indexS, tracers, rho, err1)
-! call OcnEquationOfStateJMRho(grid, displacement_type, k_displaced, indexT, indexS, tracers, rho, err2)
-
- err = err1 .or. err2
-
- !--------------------------------------------------------------------
-
- end subroutine OcnEquationOfStateRho
-
-!***********************************************************************
-!
-! routine OcnEquationOfStateInit
-!
-!> \brief Initializes ocean momentum horizontal mixing quantities
-!> \author Doug Jacobsen
-!> \date 19 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine initializes a variety of quantities related to
-!> horizontal velocity mixing in the ocean. Since a variety of
-!> parameterizations are available, this routine primarily calls the
-!> individual init routines for each parameterization.
-!
-!-----------------------------------------------------------------------
-
-
- subroutine OcnEquationOfStateInit(err)
-
- !--------------------------------------------------------------------
-
- !-----------------------------------------------------------------
- !
- ! call individual init routines for each parameterization
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err
-
- integer :: err1, err2
-
- err = 0
- ! For an isopycnal model, density should remain constant.
- ! For zlevel, calculate in-situ density
- eosON = .false.
-
- if(config_vert_grid_type.eq.'zlevel') then
- eosON = .true.
-! call OcnEquationOfStateLinearInit(err1)
-! call OcnEquationOfStateJMInit(err2)
-
- err = err1 .or. err2
- endif
-
-
- !--------------------------------------------------------------------
-
- end subroutine OcnEquationOfStateInit
-
- subroutine equation_of_state(s, grid, k_displaced, displacement_type)!{{{
+ subroutine OcnEquationOfStateRho(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.
@@ -211,7 +97,7 @@
if (config_eos_type.eq.'linear') then
- call equation_of_state_linear(s, grid, k_displaced, displacement_type)
+ call OcnEquationOfStateLinearRho(s, grid, k_displaced, displacement_type)
elseif (config_eos_type.eq.'jm') then
@@ -225,7 +111,7 @@
rho => s % rhoDisplaced % array
endif
- call equation_of_state_jm(grid, k_displaced, displacement_type, indexT, indexS, tracers, rho)
+ call OcnEquationOfStateJMRho(grid, k_displaced, displacement_type, indexT, indexS, tracers, rho)
! call equation_of_state_jm_bak(s, grid, k_displaced, displacement_type)
else
@@ -236,562 +122,57 @@
call timer_stop("equation_of_state")
- end subroutine equation_of_state!}}}
+ end subroutine OcnEquationOfStateRho!}}}
- 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.
- !
- ! 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 (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
- integer :: k_displaced
- 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(:,:), pointer :: &
- rhoPointer
- 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), dimension(:), allocatable :: &
- p, p2 ! temporary pressure scalars
-
-!-----------------------------------------------------------------------
+!***********************************************************************
!
-! UNESCO EOS constants and JMcD bulk modulus constants
+! routine OcnEquationOfStateInit
!
-!-----------------------------------------------------------------------
-
- !*** 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_bak")
-
- 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),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_bak', &
- ' k_displaced must be between 1 and nVertLevels for ', &
- 'displacement_type = absolute'
- call dmpar_abort(dminfo)
- endif
-
- if (k_displaced == 0) then
- rhoPointer => s % rho % array
- do k=1,nVertLevels
- p(k) = pRefEOS(k)
- p2(k) = p(k)*p(k)
- enddo
- else ! k_displaced /= 0
- rhoPointer => s % rhoDisplaced % array
- 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(s%index_salinity,k,iCell),smax),smin)
- TQ = max(min(tracers(s%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(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))
-
- rhoPointer(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
-
- end do
- end do
-
- 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
-
-!-----------------------------------------------------------------------
+!> \brief Initializes ocean momentum horizontal mixing quantities
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
!
-! 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
+ subroutine OcnEquationOfStateInit(err)!{{{
- !*** 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
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
- !*** from Table A1 of Jackett and McDougall
+ integer, intent(out) :: err
- real (kind=RKIND), parameter :: &
- bup0s0t0 = 1.965933e+4, &
- bup0s0t1 = 1.444304e+2, &
- bup0s0t2 = -1.706103 , &
- bup0s0t3 = 9.648704e-3, &
- bup0s0t4 = -4.190253e-5
+ integer :: err1, err2
- real (kind=RKIND), parameter :: &
- bup0s1t0 = 5.284855e+1, &
- bup0s1t1 = -3.101089e-1, &
- bup0s1t2 = 6.283263e-3, &
- bup0s1t3 = -5.084188e-5
+ err = 0
+ ! For an isopycnal model, density should remain constant.
+ ! For zlevel, calculate in-situ density
+ eosON = .false.
- real (kind=RKIND), parameter :: &
- bup0sqt0 = 3.886640e-1, &
- bup0sqt1 = 9.085835e-3, &
- bup0sqt2 = -4.619924e-4
+ if(config_vert_grid_type.eq.'zlevel') then
+ eosON = .true.
+! call OcnEquationOfStateLinearInit(err1)
+! call OcnEquationOfStateJMInit(err2)
- real (kind=RKIND), parameter :: &
- bup1s0t0 = 3.186519 , &
- bup1s0t1 = 2.212276e-2, &
- bup1s0t2 = -2.984642e-4, &
- bup1s0t3 = 1.956415e-6
+ err = err1 .or. err2
+ endif
- 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
+ end subroutine OcnEquationOfStateInit!}}}
- 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!}}}
-
-
!***********************************************************************
end module OcnEquationOfState
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateJM.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateJM.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateJM.F 2011-09-29 15:06:23 UTC (rev 1041)
@@ -0,0 +1,351 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnEquationOfStateJM
+!
+!> \brief MPAS ocean equation of state driver
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for calling
+!> the equation of state.
+!
+!-----------------------------------------------------------------------
+
+module OcnEquationOfStateJM
+
+ use grid_types
+ use configure
+ use timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnEquationOfStateJMRho, &
+ OcnEquationOfStateJMInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateJMRho
+!
+!> \brief Calls JM equation of state
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine uses a JM equation of state to update the density
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnEquationOfStateJMRho(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 OcnEquationOfStateJMRho!}}}
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateJMInit
+!
+!> \brief Initializes ocean momentum horizontal mixing quantities
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnEquationOfStateJMInit(err)!{{{
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ integer :: err1, err2
+ !--------------------------------------------------------------------
+
+ end subroutine OcnEquationOfStateJMInit!}}}
+
+!***********************************************************************
+
+end module OcnEquationOfStateJM
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateLinear.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateLinear.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateLinear.F 2011-09-29 15:06:23 UTC (rev 1041)
@@ -0,0 +1,152 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnEquationOfStateLinear
+!
+!> \brief MPAS ocean equation of state driver
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for calling
+!> the equation of state.
+!
+!-----------------------------------------------------------------------
+
+module OcnEquationOfStateLinear
+
+ use grid_types
+ use configure
+ use timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnEquationOfStateLinearRho, &
+ OcnEquationOfStateLinearInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateLinearRho
+!
+!> \brief Calls equation of state
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine uses a linear equation of state to update the density
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnEquationOfStateLinearRho(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 OcnEquationOfStateLinearRho!}}}
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateLinearInit
+!
+!> \brief Initializes ocean momentum horizontal mixing quantities
+!> \author Doug Jacobsen
+!> \date 28 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal velocity mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnEquationOfStateLinearInit(err)!{{{
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ integer :: err1, err2
+
+ err = 0
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnEquationOfStateLinearInit!}}}
+
+!***********************************************************************
+
+end module OcnEquationOfStateLinear
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F 2011-09-28 23:30:51 UTC (rev 1040)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F 2011-09-29 15:06:23 UTC (rev 1041)
@@ -395,26 +395,9 @@
!
if (.not.config_implicit_vertical_mix) then
call timer_start("OcnTendU-explicit vert mix")
+
call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
-! allocate(fluxVertTop(nVertLevels+1))
-! fluxVertTop(1) = 0.0
-! do iEdge=1,grid % nEdgesSolve
-!
-! do k=2,maxLevelEdgeTop(iEdge)
-! fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
-! * ( u(k-1,iEdge) - u(k,iEdge) ) &
-! * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
-! enddo
-! fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
-! do k=1,maxLevelEdgeTop(iEdge)
-! tend_u(k,iEdge) = tend_u(k,iEdge) &
-! + (fluxVertTop(k) - fluxVertTop(k+1)) &
-! / h_edge(k,iEdge)
-! enddo
-
-! end do
-! deallocate(fluxVertTop)
call timer_stop("OcnTendU-explicit vert mix")
endif
call timer_stop("OcnTendU")
@@ -571,36 +554,9 @@
!
if (.not.config_implicit_vertical_mix) then
call timer_start("OcnTendScalar-explicit vert diff")
+
call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
-! allocate(fluxVertTop(num_tracers,nVertLevels+1))
-! fluxVertTop(:,1) = 0.0
-! do iCell=1,nCellsSolve
-! do k=2,maxLevelCell(iCell)
-! do iTracer=1,num_tracers
-! ! compute \kappa_v d\phi/dz
-! fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
-! * ( tracers(iTracer,k-1,iCell) &
-! - tracers(iTracer,k ,iCell) ) &
-! * 2 / (h(k-1,iCell) + h(k,iCell))
-
-! enddo
-! enddo
-! fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
-
-! do k=1,maxLevelCell(iCell)
-! do iTracer=1,num_tracers
-! ! This is h d/dz( fluxVertTop) but h and dz cancel, so
-! ! reduces to delta( fluxVertTop)
-! tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
-! + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
-! enddo
-! enddo
-!print '(a,50e12.2)', 'fluxVertTop',fluxVertTop(3,1:maxLevelCell(iCell)+1)
-!print '(a,50e12.2)', 'tend_tr ',tend_tr(3,1,1:maxLevelCell(iCell))
-
-! enddo ! iCell loop
-! deallocate(fluxVertTop)
call timer_stop("OcnTendScalar-explicit vert diff")
endif
@@ -1053,9 +1009,9 @@
! For an isopycnal model, density should remain constant.
! For zlevel, calculate in-situ density
if (config_vert_grid_type.eq.'zlevel') then
- call equation_of_state(s,grid,0,'relative')
+ call OcnEquationOfStateRho(s, grid, 0, 'relative')
! mrp 110324 In order to visualize rhoDisplaced, include the following
- call equation_of_state(s, grid, 1,'relative')
+ call OcnEquationOfStateRho(s, grid, 1, 'relative')
endif
!
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnThickVadv.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnThickVadv.F 2011-09-28 23:30:51 UTC (rev 1040)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnThickVadv.F 2011-09-29 15:06:23 UTC (rev 1041)
@@ -61,7 +61,7 @@
!
!-----------------------------------------------------------------------
- subroutine OcnThickVadvTend(grid, wTop, tend, err)
+ subroutine OcnThickVadvTend(grid, wTop, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -121,7 +121,7 @@
!--------------------------------------------------------------------
- end subroutine OcnThickVadvTend
+ end subroutine OcnThickVadvTend!}}}
!***********************************************************************
!
@@ -139,9 +139,8 @@
!
!-----------------------------------------------------------------------
+ subroutine OcnThickVadvInit(err)!{{{
- subroutine OcnThickVadvInit(err)
-
!--------------------------------------------------------------------
!-----------------------------------------------------------------
@@ -156,7 +155,7 @@
!--------------------------------------------------------------------
- end subroutine OcnThickVadvInit
+ end subroutine OcnThickVadvInit!}}}
!***********************************************************************
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F 2011-09-28 23:30:51 UTC (rev 1040)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F 2011-09-29 15:06:23 UTC (rev 1041)
@@ -139,8 +139,8 @@
rhoDisplaced => s % rhoDisplaced % array
tracers => s % tracers % array
- call equation_of_state(s, grid, 0, 'relative')
- call equation_of_state(s, grid, 1, 'relative')
+ call OcnEquationOfStateRho(s, grid, 0, 'relative')
+ call OcnEquationOfStateRho(s, grid, 1, 'relative')
call OcnVmixGetRichNumbers(grid, indexT, indexS, u, h, h_edge, &
rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, err1)
</font>
</pre>