<p><b>dwj07@fsu.edu</b> 2011-09-22 10:54:00 -0600 (Thu, 22 Sep 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Replacing equation of state module (temporarily) with older<br>
equation_of_state funcitons. Need to modularize these, but for now<br>
they are in place for vertical mixing coefficients.<br>
<br>
Fixing an issue with richardson vertical mixing coefficients.<br>
<br>
Commented out the call to the coriolis module, as something causes statistics to be slightly different.<br>
Need to fix this issue as well.<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-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile 2011-09-22 16:54:00 UTC (rev 1016)
@@ -29,17 +29,16 @@
module_OcnTracerHmix.o \
module_OcnTracerHmixDel2.o \
module_OcnTracerHmixDel4.o \
- module_OcnRestoring.o \
- module_OcnEquationOfState.o \
- module_OcnEquationOfStateLinear.o \
- module_OcnEquationOfStateJM.o \
module_OcnVmix.o \
module_OcnVmixCoefsConst.o \
module_OcnVmixCoefsRich.o \
module_OcnVmixCoefsTanh.o \
+ module_OcnRestoring.o \
+ module_OcnEquationOfState.o \
module_time_integration.o \
module_global_diagnostics.o
+
all: core_hyd
core_hyd: $(OBJS)
@@ -49,7 +48,38 @@
module_advection.o:
-module_time_integration.o:
+module_time_integration.o: module_OcnThickHadv.o \
+ module_OcnThickVadv.o \
+ module_OcnVelCoriolis.o \
+ module_OcnVelVadv.o \
+ module_OcnVelHmix.o \
+ module_OcnVelHmixDel2.o \
+ module_OcnVelHmixDel4.o \
+ module_OcnVelForcing.o \
+ module_OcnVelForcingWindStress.o \
+ module_OcnVelForcingBottomDrag.o \
+ module_OcnVelPressureGrad.o \
+ module_OcnTracerVadv.o \
+ module_OcnTracerVadvSpline.o \
+ module_OcnTracerVadvSpline2.o \
+ module_OcnTracerVadvSpline3.o \
+ module_OcnTracerVadvStencil.o \
+ module_OcnTracerVadvStencil2.o \
+ module_OcnTracerVadvStencil3.o \
+ module_OcnTracerVadvStencil4.o \
+ module_OcnTracerHadv.o \
+ module_OcnTracerHadv2.o \
+ module_OcnTracerHadv3.o \
+ module_OcnTracerHadv4.o \
+ module_OcnTracerHmix.o \
+ module_OcnTracerHmixDel2.o \
+ module_OcnTracerHmixDel4.o \
+ module_OcnVmix.o \
+ module_OcnVmixCoefsConst.o \
+ module_OcnVmixCoefsRich.o \
+ module_OcnVmixCoefsTanh.o \
+ module_OcnRestoring.o \
+ module_OcnEquationOfState.o
module_global_diagnostics.o:
@@ -107,20 +137,16 @@
module_OcnRestoring.o:
-module_OcnEquationOfState.o: module_OcnEquationOfStateLinear.o module_OcnEquationOfStateJM.o
-
-module_OcnEquationOfStateLinear.o:
-
-module_OcnEquationOfStateJM.o:
-
module_OcnVmix.o: module_OcnVmixCoefsConst.o module_OcnVmixCoefsRich.o module_OcnVmixCoefsTanh.o
module_OcnVmixCoefsConst.o:
-module_OcnVmixCoefsRich.o:
+module_OcnVmixCoefsRich.o: module_OcnEquationOfState.o
module_OcnVmixCoefsTanh.o:
+module OcnEquationOfState.o:
+
module_mpas_core.o: module_advection.o \
module_OcnThickHadv.o \
module_OcnThickVadv.o \
@@ -151,13 +177,11 @@
module_OcnTracerHmixDel2.o \
module_OcnTracerHmixDel4.o \
module_OcnRestoring.o \
- module_OcnEquationOfState.o \
- module_OcnEquationOfStateLinear.o \
- module_OcnEquationOfStateJM.o \
module_OcnVmix.o \
module_OcnVmixCoefsConst.o \
module_OcnVmixCoefsRich.o \
module_OcnVmixCoefsTanh.o \
+ module_OcnEquationOfState.o \
module_time_integration.o
clean:
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfState.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -16,8 +16,9 @@
use grid_types
use configure
- use OcnEquationOfStateLinear
- use OcnEquationOfStateJM
+ use timer
+! use OcnEquationOfStateLinear
+! use OcnEquationOfStateJM
implicit none
private
@@ -36,7 +37,8 @@
!--------------------------------------------------------------------
public :: OcnEquationOfStateRho, &
- OcnEquationOfStateInit
+ OcnEquationOfStateInit, &
+ equation_of_state
!--------------------------------------------------------------------
!
@@ -119,8 +121,8 @@
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)
+! call OcnEquationOfStateLinearRho(grid, indexT, indexS, tracers, rho, err1)
+! call OcnEquationOfStateJMRho(grid, displacement_type, k_displaced, indexT, indexS, tracers, rho, err2)
err = err1 .or. err2
@@ -166,8 +168,8 @@
if(config_vert_grid_type.eq.'zlevel') then
eosON = .true.
- call OcnEquationOfStateLinearInit(err1)
- call OcnEquationOfStateJMInit(err2)
+! call OcnEquationOfStateLinearInit(err1)
+! call OcnEquationOfStateJMInit(err2)
err = err1 .or. err2
endif
@@ -177,6 +179,322 @@
end subroutine OcnEquationOfStateInit
+ 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!}}}
+
+
!***********************************************************************
end module OcnEquationOfState
Deleted: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateJM.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateJM.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateJM.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -1,372 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! 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
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Deleted: branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateLinear.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateLinear.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnEquationOfStateLinear.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -1,176 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! 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_OcnVelCoriolis.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVelCoriolis.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVelCoriolis.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -119,7 +119,7 @@
nEdgesSolve = grid % nEdgesSolve
- do iEdge=1,nEdgesSolve
+ do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -139,7 +139,6 @@
end do
end do
-
!--------------------------------------------------------------------
end subroutine OcnVelCoriolisTend
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnVmix.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmix.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmix.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -721,4 +721,4 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!vim: foldmethod=marker
+! vim: foldmethod=marker
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsConst.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsConst.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsConst.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -303,4 +303,4 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!vim: foldmethod=marker
+! vim: foldmethod=marker
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsRich.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -38,10 +38,6 @@
!
!--------------------------------------------------------------------
- private :: OcnVelVmixCoefsRich, &
- OcnTracerVmixCoefsRich, &
- OcnVmixGetRichNumbers
-
public :: OcnVmixCoefsRichBuild, &
OcnVmixCoefsRichInit
@@ -146,7 +142,7 @@
call equation_of_state(s, grid, 0, 'relative')
call equation_of_state(s, grid, 1, 'relative')
- call OcnVmixGetRichNumbers(grid, indexT, indexS, u, h, h_edge, &
+ call OcnVmixGetRichNumbers(grid, indexT, indexS, u, h, h_edge, &
rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, err1)
call OcnVelVmixCoefsRich(grid, RiTopOfEdge, h_edge, vertViscTopOfEdge, err2)
@@ -476,6 +472,7 @@
end do
! interpolate drhoTopOfCell to drhoTopOfEdge
+ drhoTopOfEdge = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -495,7 +492,7 @@
end do
! interpolate du2TopOfEdge to du2TopOfCell
- du2TopOfCell(:,:) = 0.0
+ du2TopOfCell = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -590,326 +587,10 @@
end subroutine OcnVmixCoefsRichInit!}}}
- !!--- TEMPORARY EQUATION OF STATE ROUTINES
- 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!}}}
-
!***********************************************************************
end module OcnVmixCoefsRich
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!vim: foldmethod=marker
+! vim: foldmethod=marker
Modified: branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsTanh.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsTanh.F 2011-09-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnVmixCoefsTanh.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -329,4 +329,4 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!vim: foldmethod=marker
+! vim: foldmethod=marker
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-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -15,7 +15,6 @@
use OcnTracerHmix
use OcnRestoring
- use OcnEquationOfState
use OcnVmix
type (io_output_object) :: restart_obj
@@ -97,13 +96,8 @@
call OcnTracerHmixInit(err)
call OcnRestoringInit(err)
- call OcnEquationOfStateInit(err)
call OcnVmixInit(err)
- if(err) then
- call dmpar_abort(dminfo)
- endif
-
! 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-21 21:23:00 UTC (rev 1015)
+++ branches/ocean_projects/performance/src/core_ocean/module_time_integration.F 2011-09-22 16:54:00 UTC (rev 1016)
@@ -180,7 +180,7 @@
block => domain % blocklist
do while (associated(block))
if (.not.config_implicit_vertical_mix) then
- call OcnVmixCoefs(block%mesh, provis, block%diagnostics, err)
+ call OcnVmixCoefs(block % mesh, provis, block % diagnostics, err)
end if
call compute_tend_h(block % tend, provis, block % diagnostics, block % mesh)
call compute_tend_u(block % tend, provis, block % diagnostics, block % mesh)
@@ -309,12 +309,50 @@
if (config_implicit_vertical_mix) then
call timer_start("RK4-implicit vert mix")
+ allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels), &
+ tracersTemp(num_tracers,nVertLevels))
- call OcnVmixCoefs(block%mesh, provis, block%diagnostics, err)
+ call OcnVmixCoefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
- call OcnVelVmixTendImplicit(block%mesh, dt, ke_edge, vertViscTopOfEdge,h, &
- h_edge, u, err)
-
+ !
+ ! Implicit vertical solve for momentum
+ !
+ call OcnVelVmixTendImplicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
+! do iEdge=1,nEdges
+! if (maxLevelEdgeTop(iEdge).gt.0) then
+
+ ! Compute A(k), C(k) for momentum
+ ! mrp 110315 efficiency note: for z-level, could precompute
+ ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
+ ! h_edge is computed in compute_solve_diag, and is not available yet.
+ ! This could be removed if hZLevel used instead.
+! cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+! cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+! do k=1,maxLevelEdgeTop(iEdge)
+! h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+! end do
+
+! do k=1,maxLevelEdgeTop(iEdge)-1
+! A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
+! / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
+! / h_edge(k,iEdge)
+! enddo
+! A(maxLevelEdgeTop(iEdge)) = -dt*config_bottom_drag_coeff &
+! *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+
+! C(1) = 1 - A(1)
+! do k=2,maxLevelEdgeTop(iEdge)
+! C(k) = 1 - A(k) - A(k-1)
+! enddo
+
+! call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
+
+! u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
+! u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
+
+! end if
+! end do
+
! mrp 110718 filter btr mode out of u
if (config_rk_filter_btr_mode) then
call filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
@@ -324,9 +362,34 @@
!
! Implicit vertical solve for tracers
!
-
- call OcnTracerVmixTendImplicit(block%mesh, dt, vertDiffTopOfCell, h, &
- tracers, err)
+
+ call OcnTracerVmixTendImplicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
+! do iCell=1,nCells
+ ! Compute A(k), C(k) for tracers
+ ! mrp 110315 efficiency note: for z-level, could precompute
+ ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
+! do k=1,maxLevelCell(iCell)-1
+! A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
+! / (h(k,iCell) + h(k+1,iCell)) &
+! / h(k,iCell)
+! enddo
+
+! A(maxLevelCell(iCell)) = 0.0
+
+! C(1) = 1 - A(1)
+! do k=2,maxLevelCell(iCell)
+! C(k) = 1 - A(k) - A(k-1)
+! enddo
+
+! call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
+! tracersTemp, maxLevelCell(iCell), &
+! nVertLevels,num_tracers)
+
+! tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
+! tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
+! end do
+! deallocate(A,C,uTemp,tracersTemp)
+! call timer_stop("RK4-implicit vert mix")
end if
! mrp 110725 momentum decay term
@@ -395,9 +458,9 @@
uPerp, uCorr, tracerTemp, coef
real (kind=RKIND), dimension(:), pointer :: sshNew
- integer :: num_tracers, ucorr_coef
+ integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
- u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell
+ u, h, h_edge, ke_edge, vertViscTopOfEdge, vertDiffTopOfCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
@@ -492,7 +555,7 @@
block => domain % blocklist
do while (associated(block))
if (.not.config_implicit_vertical_mix) then
- call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+ call OcnVmixCoefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
end if
call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
call enforce_boundaryEdge(block % tend, block % mesh)
@@ -1391,6 +1454,7 @@
tracers => block % state % time_levs(2) % state % tracers % array
h => block % state % time_levs(2) % state % h % array
h_edge => block % state % time_levs(2) % state % h_edge % array
+ ke_edge => block % state % time_levs(2) % state % ke_edge % array
num_tracers = block % state % time_levs(2) % state % num_tracers
vertViscTopOfEdge => block % diagnostics % vertViscTopOfEdge % array
vertDiffTopOfCell => block % diagnostics % vertDiffTopOfCell % array
@@ -1401,73 +1465,77 @@
allocate(A(block % mesh % nVertLevels),C(block % mesh % nVertLevels),uTemp(block % mesh % nVertLevels), &
tracersTemp(num_tracers,block % mesh % nVertLevels))
- call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+ call OcnVmixCoefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
!
! Implicit vertical solve for momentum
!
- do iEdge=1,block % mesh % nEdges
- if (maxLevelEdgeTop(iEdge).gt.0) then
+ call OcnVelVmixTendImplicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
+
+! do iEdge=1,block % mesh % nEdges
+! if (maxLevelEdgeTop(iEdge).gt.0) then
+
! Compute A(k), C(k) for momentum
! mrp 110315 efficiency note: for z-level, could precompute
! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
! h_edge is computed in compute_solve_diag, and is not available yet.
! This could be removed if hZLevel used instead.
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
+! cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+! cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+! do k=1,maxLevelEdgeTop(iEdge)
+! h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+! end do
- do k=1,maxLevelEdgeTop(iEdge)-1
- A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
- / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
- / h_edge(k,iEdge)
- enddo
- A(maxLevelEdgeTop(iEdge)) = 0.0
+! do k=1,maxLevelEdgeTop(iEdge)-1
+! A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
+! / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
+! / h_edge(k,iEdge)
+! enddo
+! A(maxLevelEdgeTop(iEdge)) = 0.0
- C(1) = 1 - A(1)
- do k=2,maxLevelEdgeTop(iEdge)
- C(k) = 1 - A(k) - A(k-1)
- enddo
+! C(1) = 1 - A(1)
+! do k=2,maxLevelEdgeTop(iEdge)
+! C(k) = 1 - A(k) - A(k-1)
+! enddo
- call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
+! call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
- u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
- u(maxLevelEdgeTop(iEdge)+1:block % mesh % nVertLevels,iEdge) = 0.0
+! u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
+! u(maxLevelEdgeTop(iEdge)+1:block % mesh % nVertLevels,iEdge) = 0.0
- end if
- end do
+! end if
+! end do
!
! Implicit vertical solve for tracers
!
- do iCell=1,block % mesh % nCells
+ call OcnTracerVmixTendImplicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
+! do iCell=1,block % mesh % nCells
! Compute A(k), C(k) for tracers
! mrp 110315 efficiency note: for z-level, could precompute
! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
- do k=1,maxLevelCell(iCell)-1
- A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
- / (h(k,iCell) + h(k+1,iCell)) &
- / h(k,iCell)
- enddo
+! do k=1,maxLevelCell(iCell)-1
+! A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
+! / (h(k,iCell) + h(k+1,iCell)) &
+! / h(k,iCell)
+! enddo
- A(maxLevelCell(iCell)) = 0.0
+! A(maxLevelCell(iCell)) = 0.0
- C(1) = 1 - A(1)
- do k=2,maxLevelCell(iCell)
- C(k) = 1 - A(k) - A(k-1)
- enddo
+! C(1) = 1 - A(1)
+! do k=2,maxLevelCell(iCell)
+! C(k) = 1 - A(k) - A(k-1)
+! enddo
- call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
- tracersTemp, maxLevelCell(iCell), &
- block % mesh % nVertLevels,num_tracers)
+! call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
+! tracersTemp, maxLevelCell(iCell), &
+! block % mesh % nVertLevels,num_tracers)
- tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
- tracers(:,maxLevelCell(iCell)+1:block % mesh % nVertLevels,iCell) = -1e34
- end do
- deallocate(A,C,uTemp,tracersTemp)
+! tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
+! tracers(:,maxLevelCell(iCell)+1:block % mesh % nVertLevels,iCell) = -1e34
+! end do
+! deallocate(A,C,uTemp,tracersTemp)
end if
! mrp 110725 adding momentum decay term
@@ -1746,8 +1814,28 @@
call timer_start("compute_tend_u-coriolis")
- call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
+! call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ q = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
+ q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+ end do
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + q &
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+
+ end do
+ end do
+
call timer_stop("compute_tend_u-coriolis")
!
@@ -1798,8 +1886,30 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
+ if (.not.config_implicit_vertical_mix) then
+ call timer_start("compute_tend_u-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("compute_tend_u-explicit vert mix")
+ endif
call timer_stop("compute_tend_u")
end subroutine compute_tend_u!}}}
@@ -2299,9 +2409,41 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
- call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
+ if (.not.config_implicit_vertical_mix) then
+ call timer_start("compute_scalar_tend-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("compute_scalar_tend-explicit vert diff")
+ endif
+
! mrp 110516 printing
!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&
! maxval(tend_tr(3,1,1:nCells))
@@ -2337,7 +2479,7 @@
type (mesh_type), intent(in) :: grid
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov, err
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
@@ -2350,7 +2492,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, rhoDisplaced
+ rho, temperature, salinity
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:,:), allocatable:: div_u
@@ -2385,7 +2527,6 @@
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
@@ -2736,18 +2877,14 @@
!
! equation of state
!
-! 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)
+ ! 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')
! mrp 110324 In order to visualize rhoDisplaced, include the following
call equation_of_state(s, grid, 1,'relative')
endif
-
!
! Pressure
! This section must be after computing rho
@@ -2902,310 +3039,6 @@
end subroutine compute_wtop!}}}
- subroutine compute_vertical_mix_coefficients(s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: s
- type (diagnostics_type), intent(inout) :: d
- type (mesh_type), intent(in) :: grid
-
- type (dm_info) :: dminfo
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, i, j
- integer :: nCells, nEdges, nVertices, nVertLevels, err
-
- real (kind=RKIND), dimension(:,:), allocatable:: &
- drhoTopOfCell, drhoTopOfEdge, &
- du2TopOfCell, du2TopOfEdge
- real (kind=RKIND) :: coef
- real (kind=RKIND), dimension(:,:), pointer :: &
- vertViscTopOfEdge, vertDiffTopOfCell, &
- RiTopOfEdge, RiTopOfCell, rhoDisplaced, rho, &
- kiteAreasOnVertex, h_edge, h, u
-
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel, zTopZLevel
- character :: c1*6
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexBot, maxLevelVertexTop
- integer, dimension(:,:), pointer :: &
- cellsOnEdge
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two, tracers
- real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND) :: r, h1, h2
-
- call timer_start("compute_vertical_mix_coefficients")
-
- rho => s % rho % array
- rhoDisplaced => s % rhoDisplaced % array
- 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
- RiTopOfEdge => d % RiTopOfEdge % array
- RiTopOfCell => d % RiTopOfCell % array
-
- zTopZLevel => grid % zTopZLevel % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelEdgeBot => grid % maxLevelEdgeBot % array
- cellsOnEdge => grid % cellsOnEdge % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- !
- ! Compute Richardson Number
- !
- if (config_vert_visc_type.eq.'rich'.or. &
- config_vert_diff_type.eq.'rich') then
- allocate( &
- drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges+1), &
- du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges+1))
-
- ! compute density of parcel displaced to next deeper z-level,
- ! in state % rhoDisplaced
-!maltrud make sure rho is current--check this for redundancy
-! 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)
-
- call equation_of_state(s, grid, 0, 'relative')
- call equation_of_state(s, grid, 1, 'relative')
-
- ! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
- drhoTopOfCell = 0.0
- do iCell=1,nCells
- do k=2,maxLevelCell(iCell)
- drhoTopOfCell(k,iCell) = rho(k-1,iCell) - rhoDisplaced(k-1,iCell)
- end do
- end do
-
- ! interpolate drhoTopOfCell to drhoTopOfEdge
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeTop(iEdge)
- drhoTopOfEdge(k,iEdge) = &
- (drhoTopOfCell(k,cell1) + &
- drhoTopOfCell(k,cell2))/2
- end do
- end do
-
- ! du2TopOfEdge(k) = $u_{k-1}-u_k$
- du2TopOfEdge=0.0
- do iEdge=1,nEdges
- do k=2,maxLevelEdgeTop(iEdge)
- du2TopOfEdge(k,iEdge) = (u(k-1,iEdge) - u(k,iEdge))**2
- end do
- end do
-
- ! interpolate du2TopOfEdge to du2TopOfCell
- du2TopOfCell(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
- du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- end do
- end do
- do iCell = 1,nCells
- do k = 2,maxLevelCell(iCell)
- du2TopOfCell(k,iCell) = du2TopOfCell(k,iCell) / areaCell(iCell)
- end do
- end do
-
- ! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
- ! coef = -g/rho_0/2
- RiTopOfEdge = 0.0
- coef = -gravity/1000.0/2.0
- do iEdge = 1,nEdges
- do k = 2,maxLevelEdgeTop(iEdge)
- RiTopOfEdge(k,iEdge) = coef*drhoTopOfEdge(k,iEdge) &
- *(h_edge(k-1,iEdge)+h_edge(k,iEdge)) &
- / (du2TopOfEdge(k,iEdge) + 1e-20)
- end do
- end do
-
- ! compute RiTopOfCell using drhoTopOfCell and du2TopOfCell
- ! coef = -g/rho_0/2
- RiTopOfCell = 0.0
- coef = -gravity/1000.0/2.0
- do iCell = 1,nCells
- do k = 2,maxLevelCell(iCell)
- RiTopOfCell(k,iCell) = coef*drhoTopOfCell(k,iCell) &
- *(h(k-1,iCell)+h(k,iCell)) &
- / (du2TopOfCell(k,iCell) + 1e-20)
- end do
- end do
-
- deallocate(drhoTopOfCell, drhoTopOfEdge, &
- du2TopOfCell, du2TopOfEdge)
- endif
-
- !
- ! Compute vertical viscosity
- !
- if (config_vert_visc_type.eq.'const') then
-
- vertViscTopOfEdge = config_vert_visc
-
- elseif (config_vert_visc_type.eq.'tanh') then
-
- if (config_vert_grid_type.ne.'zlevel') then
- write(0,*) 'Abort: config_vert_visc_type.eq.tanh may only', &
- ' use config_vert_grid_type of zlevel at this time'
- call dmpar_abort(dminfo)
- endif
-
- do k=1,nVertLevels+1
- vertViscTopOfEdge(k,:) = -(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
- *tanh(-(zTopZLevel(k)-config_ZMid_tanh) &
- /config_zWidth_tanh) &
- + (config_max_visc_tanh+config_min_visc_tanh)/2
- end do
-
- elseif (config_vert_visc_type.eq.'rich') then
-
- vertViscTopOfEdge = 0.0
- do iEdge = 1,nEdges
- do k = 2,maxLevelEdgeTop(iEdge)
- ! mrp 110324 efficiency note: this if is inside iEdge and k loops.
- ! Perhaps there is a more efficient way to do this.
- if (RiTopOfEdge(k,iEdge)>0.0) then
- vertViscTopOfEdge(k,iEdge) = config_bkrd_vert_visc &
- + config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
- ! maltrud do limiting of coefficient--should not be necessary
- ! also probably better logic could be found
- if (vertViscTopOfEdge(k,iEdge) > config_convective_visc) then
- if( config_implicit_vertical_mix) then
- vertViscTopOfEdge(k,iEdge) = config_convective_visc
- else
- vertViscTopOfEdge(k,iEdge) = &
- ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
- end if
- end if
- else
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- if (config_implicit_vertical_mix) then
- ! for Ri<0 and implicit mix, use convective diffusion
- vertViscTopOfEdge(k,iEdge) = config_convective_visc
- else
- ! for Ri<0 and explicit vertical mix,
- ! use maximum diffusion allowed by CFL criterion
- ! mrp 110324 efficiency note: for z-level, could use fixed
- ! grid array hMeanTopZLevel and compute maxdiff on startup.
- vertViscTopOfEdge(k,iEdge) = &
- ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
- end if
- end if
- end do
- end do
-
- else
-
- write(0,*) 'Abort: unrecognized config_vert_visc_type'
- call dmpar_abort(dminfo)
-
- endif
-
- !
- ! Compute vertical tracer diffusion
- !
- if (config_vert_diff_type.eq.'const') then
-
- vertDiffTopOfCell = config_vert_diff
-
- elseif (config_vert_diff_type.eq.'tanh') then
-
- if (config_vert_grid_type.ne.'zlevel') then
- write(0,*) 'Abort: config_vert_diff_type.eq.tanh may only', &
- ' use config_vert_grid_type of zlevel at this time'
- call dmpar_abort(dminfo)
- endif
-
- do k=1,nVertLevels+1
- vertDiffTopOfCell(k,:) = -(config_max_diff_tanh-config_min_diff_tanh)/2.0 &
- *tanh(-(zTopZLevel(k)-config_ZMid_tanh) &
- /config_zWidth_tanh) &
- + (config_max_diff_tanh+config_min_diff_tanh)/2
- end do
-
- elseif (config_vert_diff_type.eq.'rich') then
-
- vertDiffTopOfCell = 0.0
- coef = -gravity/1000.0/2.0
- do iCell = 1,nCells
- do k = 2,maxLevelCell(iCell)
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- ! Perhaps there is a more efficient way to do this.
- if (RiTopOfCell(k,iCell)>0.0) then
- vertDiffTopOfCell(k,iCell) = config_bkrd_vert_diff &
- + (config_bkrd_vert_visc &
- + config_rich_mix / (1.0 + 5.0*RiTopOfCell(k,iCell))**2) &
- / (1.0 + 5.0*RiTopOfCell(k,iCell))
- ! maltrud do limiting of coefficient--should not be necessary
- ! also probably better logic could be found
- if (vertDiffTopOfCell(k,iCell) > config_convective_diff) then
- if (config_implicit_vertical_mix) then
- vertDiffTopOfCell(k,iCell) = config_convective_diff
- else
- vertDiffTopOfCell(k,iCell) = &
- ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
- end if
- end if
- else
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- if (config_implicit_vertical_mix) then
- ! for Ri<0 and implicit mix, use convective diffusion
- vertDiffTopOfCell(k,iCell) = config_convective_diff
- else
- ! for Ri<0 and explicit vertical mix,
- ! use maximum diffusion allowed by CFL criterion
- ! mrp 110324 efficiency note: for z-level, could use fixed
- ! grid array hMeanTopZLevel and compute maxdiff on startup.
- vertDiffTopOfCell(k,iCell) = &
- ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
- end if
- end if
- end do
- end do
-
- else
-
- write(0,*) 'Abort: unrecognized config_vert_diff_type'
- call dmpar_abort(dminfo)
-
- endif
-
- call timer_stop("compute_vertical_mix_coefficients")
- end subroutine compute_vertical_mix_coefficients!}}}
-
subroutine enforce_boundaryEdge(tend, grid)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
@@ -3251,421 +3084,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.
-! A is an nxn matrix, with:
-! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
-! b diagonal, filled from 1:n
-! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
-!
-! Input: a,b,c,r,n
-!
-! Output: x
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- integer,intent(in) :: n
- real (KIND=RKIND), dimension(n), intent(in) :: a,b,c,r
- real (KIND=RKIND), dimension(n), intent(out) :: x
- real (KIND=RKIND), dimension(n) :: bTemp,rTemp
- real (KIND=RKIND) :: m
- integer i
-
- call timer_start("tridiagonal_solve")
-
- ! Use work variables for b and r
- bTemp(1) = b(1)
- rTemp(1) = r(1)
-
- ! First pass: set the coefficients
- do i = 2,n
- m = a(i-1)/bTemp(i-1)
- bTemp(i) = b(i) - m*c(i-1)
- rTemp(i) = r(i) - m*rTemp(i-1)
- end do
-
- x(n) = rTemp(n)/bTemp(n)
- ! Second pass: back-substition
- do i = n-1, 1, -1
- x(i) = (rTemp(i) - c(i)*x(i+1))/bTemp(i)
- end do
-
- call timer_stop("tridiagonal_solve")
-
-end subroutine tridiagonal_solve!}}}
-
-subroutine tridiagonal_solve_mult(a,b,c,r,x,n,nDim,nSystems)!{{{
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Solve the matrix equation Ax=r for x, where A is tridiagonal.
-! A is an nxn matrix, with:
-! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
-! b diagonal, filled from 1:n
-! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
-!
-! Input: a,b,c,r,n
-!
-! Output: x
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- integer,intent(in) :: n, nDim, nSystems
- real (KIND=RKIND), dimension(n), intent(in) :: a,b,c
- real (KIND=RKIND), dimension(nSystems,nDim), intent(in) :: r
- real (KIND=RKIND), dimension(nSystems,nDim), intent(out) :: x
- real (KIND=RKIND), dimension(n) :: bTemp
- real (KIND=RKIND), dimension(nSystems,n) :: rTemp
- real (KIND=RKIND) :: m
- integer i,j
-
- call timer_start("tridiagonal_solve_mult")
-
- ! Use work variables for b and r
- bTemp(1) = b(1)
- do j = 1,nSystems
- rTemp(j,1) = r(j,1)
- end do
-
- ! First pass: set the coefficients
- do i = 2,n
- m = a(i-1)/bTemp(i-1)
- bTemp(i) = b(i) - m*c(i-1)
- do j = 1,nSystems
- rTemp(j,i) = r(j,i) - m*rTemp(j,i-1)
- end do
- end do
-
- do j = 1,nSystems
- x(j,n) = rTemp(j,n)/bTemp(n)
- end do
- ! Second pass: back-substition
- do i = n-1, 1, -1
- do j = 1,nSystems
- x(j,i) = (rTemp(j,i) - c(i)*x(j,i+1))/bTemp(i)
- end do
- end do
-
- call timer_stop("tridiagonal_solve_mult")
-
-end subroutine tridiagonal_solve_mult!}}}
-
end module time_integration
-!:vim foldmethod=marker
+! vim: foldmethod=marker
</font>
</pre>