<p><b>mpetersen@lanl.gov</b> 2011-03-14 15:57:46 -0600 (Mon, 14 Mar 2011)</p><p>BRANCH COMMIT: Added Richardson-number based vertical mixing (Pacanowski and Philander)<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/imp_vert_mix_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/namelist.input.ocean 2011-03-11 23:38:00 UTC (rev 757)
+++ branches/ocean_projects/imp_vert_mix_mrp/namelist.input.ocean 2011-03-14 21:57:46 UTC (rev 758)
@@ -1,10 +1,10 @@
&sw_model
config_test_case = 0
config_time_integration = 'RK4'
- config_dt = 90.0
- config_ntimesteps = 1920000
- config_output_interval = 19200
- config_stats_interval = 1920
+ config_dt = 120.0
+ config_ntimesteps = 100
+ config_output_interval = 100
+ config_stats_interval = 100
/
&io
@@ -14,35 +14,45 @@
/
&restart
- config_restart_interval = 115200
+ config_restart_interval = 10000000
config_do_restart = .false.
- config_restart_time = 31104000
+ config_restart_time = 1036800.0
/
&grid
- config_vert_grid_type = 'isopycnal'
- config_rho0 = 1015
+ config_vert_grid_type = 'zlevel'
+ config_rho0 = 1000
/
&hmix
- config_h_mom_eddy_visc2 = 0.0
- config_h_mom_eddy_visc4 = 5.0e8
- config_h_tracer_eddy_diff2 = 10.0
+ config_h_mom_eddy_visc2 = 1.0e5
+ config_h_mom_eddy_visc4 = 0.0
+ config_h_tracer_eddy_diff2 = 1.0e4
config_h_tracer_eddy_diff4 = 0.0
/
&vmix
- config_vert_visc_type = 'const'
- config_vert_diff_type = 'const'
- config_vmixTanhViscMax = 2.5e-1
- config_vmixTanhViscMin = 1.0e-4
- config_vmixTanhDiffMax = 2.5e-2
- config_vmixTanhDiffMin = 1.0e-5
- config_vmixTanhZMid = -100
- config_vmixTanhZWidth = 100
- config_vert_viscosity = 1.0e-4
- config_vert_diffusion = 1.0e-4
+ config_vert_visc_type = 'rich'
+ config_vert_diff_type = 'rich'
+ config_implicit_vertical_mix = .true.
/
+&vmix_const
+ config_vert_viscosity = 2.5e-5
+ config_vert_diffusion = 2.5e-5
+/
+&vmix_rich
+ config_bkrd_vert_visc = 1.0e-4
+ config_bkrd_vert_diff = 1.0e-5
+ config_rich_mix = 50.0
+/
+&vmix_tanh
+ config_max_visc_tanh = 2.5e-1
+ config_min_visc_tanh = 1.0e-4
+ config_max_diff_tanh = 2.5e-2
+ config_min_diff_tanh = 1.0e-5
+ config_zMid_tanh = -100
+ config_zWidth_tanh = 100
+/
&eos
- config_eos_type = 'linear'
+ config_eos_type = 'jm'
/
&advection
config_vert_tracer_adv = 'spline'
Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry 2011-03-11 23:38:00 UTC (rev 757)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry 2011-03-14 21:57:46 UTC (rev 758)
@@ -23,14 +23,18 @@
namelist real hmix config_apvm_upwinding 0.5
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
-namelist real vmix config_vert_viscosity 2.5e-4
-namelist real vmix config_vert_diffusion 2.5e-5
-namelist real vmix config_vmixTanhViscMax 2.5e-1
-namelist real vmix config_vmixTanhViscMin 1.0e-4
-namelist real vmix config_vmixTanhDiffMax 2.5e-2
-namelist real vmix config_vmixTanhDiffMin 1.0e-5
-namelist real vmix config_vmixTanhZMid -100
-namelist real vmix config_vmixTanhZWidth 100
+namelist logical vmix config_implicit_vertical_mix .true.
+namelist real vmix_const config_vert_viscosity 2.5e-4
+namelist real vmix_const config_vert_diffusion 2.5e-5
+namelist real vmix_rich config_bkrd_vert_visc 1.0e-4
+namelist real vmix_rich config_bkrd_vert_diff 1.0e-5
+namelist real vmix_rich config_rich_mix 50
+namelist real vmix_tanh config_max_visc_tanh 2.5e-1
+namelist real vmix_tanh config_min_visc_tanh 1.0e-4
+namelist real vmix_tanh config_max_diff_tanh 2.5e-2
+namelist real vmix_tanh config_min_diff_tanh 1.0e-5
+namelist real vmix_tanh config_zMid_tanh -100
+namelist real vmix_tanh config_zWidth_tanh 100
namelist character eos config_eos_type linear
namelist character advection config_vert_tracer_adv stencil
namelist integer advection config_vert_tracer_adv_order 4
@@ -177,6 +181,7 @@
var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 o pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
+var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 o rhoDisplaced state - -
# Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -193,3 +198,10 @@
var persistent real volumeEdgeGlobal ( Time ) 2 o volumeEdgeGlobal state - -
var persistent real CFLNumberGlobal ( Time ) 2 o CFLNumberGlobal state - -
+# Diagnostics fields, only one time level required
+var persistent real RiTopOfCell ( nVertLevelsP1 nCells ) 0 o RiTopOfCell diagnostics - -
+var persistent real RiTopOfEdge ( nVertLevelsP1 nEdges ) 0 o RiTopOfEdge diagnostics - -
+var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges ) 0 o vertViscTopOfEdge diagnostics - -
+var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells ) 0 o vertDiffTopOfCell diagnostics - -
+
+
Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-11 23:38:00 UTC (rev 757)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-14 21:57:46 UTC (rev 758)
@@ -143,8 +143,11 @@
block => domain % blocklist
do while (associated(block))
- call compute_tend(block % tend, provis, block % mesh)
- call compute_scalar_tend(block % tend, provis, block % mesh)
+ if (.not.config_implicit_vertical_mix) then
+ call compute_vertical_mix_coefficients(provis, block % diagnostics, block % mesh)
+ end if
+ call compute_tend(block % tend, provis, block % diagnostics, block % mesh)
+ call compute_scalar_tend(block % tend, provis, block % diagnostics, block % mesh)
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
end do
@@ -189,6 +192,7 @@
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
call compute_solve_diagnostics(dt, provis, block % mesh)
+
block => block % next
end do
end if
@@ -225,13 +229,33 @@
!
block => domain % blocklist
do while (associated(block))
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % maxLevelCell % array(iCell)
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- / block % state % time_levs(2) % state % h % array(k,iCell)
+
+ if (config_implicit_vertical_mix) then
+
+ call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+
+ do iCell=1,block % mesh % nCells
+ !N=maxLevelCell(iCell)
+
+ ! do k=1,maxLevelCell(iCell)
+ ! Compute A(k), C(k), R(k) for momentum
+ ! enddo
+ !call tridiagonal_solve(u(:,iCell),A(2:N),C(1:N),A(1:N-1),R)
+
+ ! do k=1,maxLevelCell(iCell)
+ ! Compute A(k), C(k), R(iTracer,k) for tracers
+ !enddo
+ !call tridiagonal_solve_mult(tracers(:,:,iCell),A(2:N),C(1:N),A(1:N-1),R)
end do
- end do
+ else
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % maxLevelCell % array(iCell)
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ / block % state % time_levs(2) % state % h % array(k,iCell)
+ end do
+ end do
+ end if
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
@@ -249,7 +273,7 @@
end subroutine rk4
- subroutine compute_tend(tend, s, grid)
+ subroutine compute_tend(tend, s, d, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -263,6 +287,7 @@
type (tend_type), intent(inout) :: tend
type (state_type), intent(in) :: s
+ type (diagnostics_type), intent(in) :: d
type (mesh_type), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
@@ -277,7 +302,7 @@
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- MontPot, wTop, divergence
+ MontPot, wTop, divergence, vertViscTopOfEdge
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
@@ -286,7 +311,7 @@
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge, vertViscTop
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
@@ -309,6 +334,7 @@
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
pressure => s % pressure % array
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -640,33 +666,12 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- allocate(vertViscTop(nVertLevels+1))
- if (config_vert_visc_type.eq.'const') then
- vertViscTop = config_vert_viscosity
- 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
- vertViscTop(k) = -(config_vmixTanhViscMax-config_vmixTanhViscMin)/2.0 &
- *tanh(-(zTopZLevel(k)-config_vmixTanhZMid) &
- /config_vmixTanhZWidth) &
- + (config_vmixTanhViscMax+config_vmixTanhViscMin)/2
- enddo
- else
- write(0,*) 'Abort: unrecognized config_vert_visc_type'
- call dmpar_abort(dminfo)
- endif
-
allocate(fluxVertTop(nVertLevels+1))
fluxVertTop(1) = 0.0
do iEdge=1,grid % nEdgesSolve
do k=2,maxLevelEdgeTop(iEdge)
- fluxVertTop(k) = vertViscTop(k) &
+ fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
* 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
enddo
@@ -679,12 +684,12 @@
enddo
end do
- deallocate(fluxVertTop, vertViscTop)
+ deallocate(fluxVertTop)
end subroutine compute_tend
- subroutine compute_scalar_tend(tend, s, grid)
+ subroutine compute_scalar_tend(tend, s, d, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -697,6 +702,7 @@
type (tend_type), intent(inout) :: tend
type (state_type), intent(in) :: s
+ type (diagnostics_type), intent(in) :: d
type (mesh_type), intent(in) :: grid
integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
@@ -706,7 +712,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,wTop, h_edge
+ u,h,wTop, h_edge, vertDiffTopOfCell
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:,:), pointer :: boundaryEdge
@@ -717,7 +723,7 @@
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
hRatioZLevelK, hRatioZLevelKm1
- real (kind=RKIND), dimension(:), allocatable:: vertDiffTop, tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
+ real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
posZTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
@@ -733,6 +739,7 @@
wTop => s % wTop % array
tracers => s % tracers % array
h_edge => s % h_edge % array
+ vertDiffTopOfCell => d % vertDiffTopOfCell % array
tend_tr => tend % tracers % array
@@ -1185,27 +1192,6 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
- allocate(vertDiffTop(nVertLevels+1))
- if (config_vert_diff_type.eq.'const') then
- vertDiffTop = config_vert_diffusion
- 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
- vertDiffTop(k) = -(config_vmixTanhDiffMax-config_vmixTanhDiffMin)/2.0 &
- *tanh(-(zTopZLevel(k)-config_vmixTanhZMid) &
- /config_vmixTanhZWidth) &
- + (config_vmixTanhDiffMax+config_vmixTanhDiffMin)/2
- enddo
- else
- write(0,*) 'Abort: unrecognized config_vert_diff_type'
- call dmpar_abort(dminfo)
- endif
-
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
do iCell=1,nCellsSolve
@@ -1213,7 +1199,7 @@
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
- fluxVertTop(iTracer,k) = vertDiffTop(k) &
+ fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
* (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
* 2 / (h(k-1,iCell) + h(k,iCell))
enddo
@@ -1230,7 +1216,7 @@
enddo
enddo ! iCell loop
- deallocate(fluxVertTop, vertDiffTop)
+ deallocate(fluxVertTop)
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
@@ -1735,7 +1721,296 @@
end subroutine compute_solve_diagnostics
+ 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
+
+
+! mrp clean out these after subroutine complete
+ 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
+
+ real (kind=RKIND), dimension(:,:), allocatable:: &
+ drhoTopOfCell, drhoTopOfEdge, &
+ du2TopOfCell, du2TopOfEdge
+ real (kind=RKIND) :: coef
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ vertViscTopOfEdge, vertDiffTopOfCell, &
+ RiTopOfEdge, RiTopOfCell, rhoDisplaced, &
+ kiteAreasOnVertex, h_edge, h, u
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ hZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ v, w, pressure,&
+ circulation, vorticity, ke, ke_edge, MontPot, wTop, &
+ pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
+ temperature, salinity
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND), dimension(:), allocatable:: pTop
+ real (kind=RKIND), dimension(:,:), allocatable:: div_u
+ character :: c1*6
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: r, h1, h2
+
+ rhoDisplaced => s % rhoDisplaced % array
+
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+ vertDiffTopOfCell => d % vertDiffTopOfCell % array
+ RiTopOfEdge => d % RiTopOfEdge % array
+ RiTopOfCell => d % RiTopOfCell % array
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ pv_vertex => s % pv_vertex % array
+ pv_cell => s % pv_cell % array
+ gradPVn => s % gradPVn % array
+ gradPVt => s % gradPVt % array
+ tracers => s % tracers % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+
+
+ zTopZLevel => grid % zTopZLevel % array
+
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ hZLevel => grid % hZLevel % array
+ deriv_two => grid % deriv_two % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexTop => grid % maxLevelVertexTop % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
+
+ boundaryEdge => grid % boundaryEdge % array
+ boundaryCell => grid % boundaryCell % array
+
+ !
+ ! 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 surface, $\rho^*$,
+ ! in state % rhoDisplaced
+ call equation_of_state(s, grid, 1)
+
+ ! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
+ drhoTopOfCell = 0.0
+ do iCell=1,nCells
+ do k=2,maxLevelCell(iCell)
+ drhoTopOfCell(k,iCell) = rhoDisplaced(k-1,iCell) - rhoDisplaced(k,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
+
+!test this interp
+
+ ! 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
+
+
+! what about boundary update?
+
+
+ deallocate(drhoTopOfCell, drhoTopOfEdge, &
+ du2TopOfCell, du2TopOfEdge)
+ endif
+
+ if (config_vert_visc_type.eq.'const') then
+
+ vertViscTopOfEdge = config_vert_viscosity
+
+ 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)
+ vertViscTopOfEdge(k,iEdge) = config_bkrd_vert_visc &
+ + config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
+ end do
+ end do
+
+ else
+
+ write(0,*) 'Abort: unrecognized config_vert_visc_type'
+ call dmpar_abort(dminfo)
+
+ endif
+
+ if (config_vert_diff_type.eq.'const') then
+
+ vertDiffTopOfCell = config_vert_diffusion
+
+ 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)
+ 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))**2
+ end do
+ end do
+
+ else
+
+ write(0,*) 'Abort: unrecognized config_vert_diff_type'
+ call dmpar_abort(dminfo)
+
+ endif
+
+
+ end subroutine compute_vertical_mix_coefficients
+
+
subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
@@ -1779,7 +2054,7 @@
end subroutine enforce_boundaryEdge
- subroutine equation_of_state(s, grid,k_displaced)
+ subroutine equation_of_state(s, grid, k_displaced)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This module contains routines necessary for computing the density
! from model temperature and salinity using an equation of state.
@@ -1787,9 +2062,9 @@
! 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
+ ! 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
@@ -1806,14 +2081,13 @@
integer :: nCells, iCell, k
type (dm_info) :: dminfo
- rho => s % rho % array
- tracers => s % tracers % array
+ if (config_eos_type.eq.'linear') then
- maxLevelCell => grid % maxLevelCell % array
- nCells = grid % nCells
+ rho => s % rho % array
+ tracers => s % tracers % array
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
- if (config_eos_type.eq.'linear') then
-
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
! Linear equation of state
@@ -1874,7 +2148,7 @@
real (kind=RKIND), dimension(:), pointer :: &
zMidZLevel, pRefEOS
real (kind=RKIND), dimension(:,:), pointer :: &
- rho
+ rhoPointer
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: maxLevelCell
@@ -1966,7 +2240,6 @@
bup2s1t1 = 6.128773e-8, &
bup2s1t2 = 6.207323e-10
- rho => s % rho % array
tracers => s % tracers % array
nCells = grid % nCells
@@ -2002,11 +2275,13 @@
! adiabatically displaced from its original level to level
! k_displaced.
if (k_displaced.le.0) then
+ rhoPointer => s % rho % array
do k=1,nVertLevels
p(k) = pRefEOS(k)
p2(k) = p(k)*p(k)
enddo
elseif (k_displaced.le.nVertLevels) then
+ rhoPointer => s % rhoDisplaced % array
do k=1,nVertLevels
p(k) = pRefEOS(k_displaced)
p2(k) = p(k)*p(k)
@@ -2058,7 +2333,7 @@
DENOMK = 1.0/(BULK_MOD - p(k))
- rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+ rhoPointer(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
end do
end do
</font>
</pre>