<p><b>dwj07@fsu.edu</b> 2011-09-19 13:44:04 -0600 (Mon, 19 Sep 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding modules for equation of state.<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-19 18:46:12 UTC (rev 1009)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-19 19:44:04 UTC (rev 1010)
@@ -30,6 +30,9 @@
module_OcnTracerHmixDel2.o \
module_OcnTracerHmixDel4.o \
module_OcnRestoring.o \
+ module_OcnEquationOfState.o \
+ module_OcnEquationOfStateLinear.o \
+ module_OcnEquationOfStateJM.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -100,6 +103,12 @@
module_OcnRestoring.o:
+module_OcnEquationOfState.o: module_OcnEquationOfStateLinear.o module_OcnEquationOfStateJM.o
+
+module_OcnEquationOfStateLinear.o:
+
+module_OcnEquationOfStateJM.o:
+
module_mpas_core.o: module_advection.o \
module_OcnThickHadv.o \
module_OcnThickVadv.o \
@@ -130,6 +139,9 @@
module_OcnTracerHmixDel2.o \
module_OcnTracerHmixDel4.o \
module_OcnRestoring.o \
+ module_OcnEquationOfState.o \
+ module_OcnEquationOfStateLinear.o \
+ module_OcnEquationOfStateJM.o \
module_time_integration.o
clean:
Added: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-19 19:44:04 UTC (rev 1010)
@@ -0,0 +1,184 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnEquationOfState
+!
+!> \brief MPAS ocean equation of state driver
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for calling
+!> the equation of state.
+!
+!-----------------------------------------------------------------------
+
+module OcnEquationOfState
+
+ use grid_types
+ use configure
+ use OcnEquationOfStateLinear
+ use OcnEquationOfStateJM
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnEquationOfStateRho, &
+ OcnEquationOfStateInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: eosON
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnEquationOfState
+!
+!> \brief Calls equation of state
+!> \author Doug Jacobsen
+!> \date 19 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
+
+!***********************************************************************
+
+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-19 19:44:04 UTC (rev 1010)
@@ -0,0 +1,372 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnEquationOfStateJM
+!
+!> \brief MPAS ocean Jackett and McDougall equation of state
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the Jackett and McDougall
+!> 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
+ !
+ !--------------------------------------------------------------------
+
+ logical :: eosJMOn
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateJM
+!
+!> \brief Calls the Jackett and McDougall equation of state
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine updates the density with a JM equation of state
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnEquationOfStateJMRho(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(out) :: &
+ rho !< Output: density
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ integer :: k_test, k_ref
+ integer :: iEdge, iCell, iVertex, k
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+
+ 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
+
+ !-----------------------------------------------------------------
+ !
+ ! 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.eosJMOn) return
+
+ call timer_start("equation_of_state-jm")
+
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+ maxLevelCell => grid % maxLevelCell % array
+ 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'
+ err = 1
+ 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 Jackett and McDougall equation of state
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a JM equation of state.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnEquationOfStateJMInit(err)
+
+ !--------------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+ eosJMOn = .false.
+
+ if(config_eos_type.eq.'JM') then
+ eosJMOn = .true.
+ endif
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnEquationOfStateJMInit
+
+!***********************************************************************
+
+end module OcnEquationOfStateJM
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
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-19 19:44:04 UTC (rev 1010)
@@ -0,0 +1,176 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnEquationOfStateLinear
+!
+!> \brief MPAS ocean equation of state driver
+!> \author Doug Jacobsen
+!> \date 19 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
+ !
+ !--------------------------------------------------------------------
+
+ logical :: linearOn
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateLinear
+!
+!> \brief Calls equation of state
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine updates the density with a linear equation of state
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnEquationOfStateLinearRho(grid, indexT, indexS, tracers, rho, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: tracers
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ integer, intent(in) :: indexT, indexS
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ rho !< Output: density
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iCell, k, nCells
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ !-----------------------------------------------------------------
+ !
+ ! 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.linearOn) return
+
+ call timer_start("equation_of_state-linear")
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(indexT,k,iCell) &
+ + 7.6e-4*tracers(indexS,k,iCell))
+ end do
+ end do
+ call timer_start("equation_of_state-linear")
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnEquationOfStateLinearRho
+
+!***********************************************************************
+!
+! routine OcnEquationOfStateLinearInit
+!
+!> \brief Initializes ocean linear equation of state
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a linear equation of state.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnEquationOfStateLinearInit(err)
+
+ !--------------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ err = 0
+ linearOn = .false.
+
+ if(config_eos_type.eq.'linear') then
+ linearOn = .true.
+ endif
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnEquationOfStateLinearInit
+
+!***********************************************************************
+
+end module OcnEquationOfStateLinear
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F 2011-09-19 18:46:12 UTC (rev 1009)
+++ branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F 2011-09-19 19:44:04 UTC (rev 1010)
@@ -15,6 +15,8 @@
use OcnTracerHmix
use OcnRestoring
+ use OcnEquationOfState
+
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -94,6 +96,8 @@
call OcnTracerHmixInit(err)
call OcnRestoringInit(err)
+ call OcnEquationOfStateInit(err)
+
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
! input arguement into mpas_init. Ask about that later. For now, there will be
! no initial statistics write.
Modified: branches/ocean_projects/performance/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_time_integration.F 2011-09-19 18:46:12 UTC (rev 1009)
+++ branches/ocean_projects/performance/src/core_ocean/module_time_integration.F 2011-09-19 19:44:04 UTC (rev 1010)
@@ -22,6 +22,8 @@
use OcnTracerHmix
use OcnRestoring
+ use OcnEquationOfState
+
contains
subroutine timestep(domain, dt, timeStamp)!{{{
@@ -2446,7 +2448,7 @@
type (mesh_type), intent(in) :: grid
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov, err
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
@@ -2459,7 +2461,7 @@
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- rho, temperature, salinity
+ rho, temperature, salinity, rhoDisplaced
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:,:), allocatable:: div_u
@@ -2494,6 +2496,7 @@
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
rho => s % rho % array
+ rhoDisplaced => s % rhoDisplaced % array
tracers => s % tracers % array
MontPot => s % MontPot % array
pressure => s % pressure % array
@@ -2844,13 +2847,11 @@
!
! equation of state
!
- ! 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(grid, 'relative', 0, s%index_temperature, s%index_salinity, &
+ tracers, rho, err)
! mrp 110324 In order to visualize rhoDisplaced, include the following
- call equation_of_state(s, grid, 1,'relative')
- endif
+ call OcnEquationOfStateRho(grid, 'relative', 1, s%index_temperature, s%index_salinity, &
+ tracers, rhoDisplaced, err)
!
! Pressure
@@ -3024,7 +3025,7 @@
type (dm_info) :: dminfo
integer :: iEdge, iCell, iVertex, k, cell1, cell2, i, j
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, err
real (kind=RKIND), dimension(:,:), allocatable:: &
drhoTopOfCell, drhoTopOfEdge, &
@@ -3046,7 +3047,7 @@
integer, dimension(:,:), pointer :: &
cellsOnEdge
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two, tracers
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: r, h1, h2
@@ -3057,6 +3058,7 @@
u => s % u % array
h => s % h % array
h_edge => s % h_edge % array
+ tracers => s % tracers % array
vertViscTopOfEdge => d % vertViscTopOfEdge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
@@ -3089,8 +3091,11 @@
! compute density of parcel displaced to next deeper z-level,
! in state % rhoDisplaced
!maltrud make sure rho is current--check this for redundancy
- call equation_of_state(s, grid, 0, 'relative')
- call equation_of_state(s, grid, 1, 'relative')
+ call OcnEquationOfStateRho(grid, 'relative', 0, s%index_temperature, s%index_salinity, &
+ tracers, rho, err)
+ ! mrp 110324 In order to visualize rhoDisplaced, include the following
+ call OcnEquationOfStateRho(grid, 'relative', 1, s%index_temperature, s%index_salinity, &
+ tracers, rhoDisplaced, err)
! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
drhoTopOfCell = 0.0
@@ -3348,321 +3353,6 @@
end subroutine enforce_boundaryEdge!}}}
- subroutine equation_of_state(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")
-
- if (config_eos_type.eq.'linear') then
-
- 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
-
- elseif (config_eos_type.eq.'jm') then
-
- call equation_of_state_jm(s, grid, k_displaced, displacement_type)
-
- else
- print *, ' Incorrect choice of config_eos_type:',&
- config_eos_type
- call dmpar_abort(dminfo)
- endif
-
- call timer_stop("equation_of_state")
-
- end subroutine equation_of_state!}}}
-
- subroutine equation_of_state_jm(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
-!
-!-----------------------------------------------------------------------
-
- !*** 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")
-
- 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', &
- ' 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")
-
- end subroutine equation_of_state_jm!}}}
-
subroutine tridiagonal_solve(a,b,c,r,x,n)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solve the matrix equation Ax=r for x, where A is tridiagonal.
</font>
</pre>