<p><b>mpetersen@lanl.gov</b> 2011-05-31 13:41:50 -0600 (Tue, 31 May 2011)</p><p>BRANCH COMMIT: Merged ocean_projects/imp_vert_mix_mrp branch into ocean_projects/timesplitting_mrp. Branch compiles and runs, KE matches to 6 digits after 10 steps compared to previous commit. The test was with config_bottom_drag_coeff=0.0, config_vert_visc_type = 'const', and config_vert_diff_type = 'const', so I expected an exact match.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/ocean_projects/timesplitting_mrp
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/ocean_projects/vert_adv_mrp:704-745
+ /branches/ocean_projects/imp_vert_mix_mrp:827-863
/branches/ocean_projects/vert_adv_mrp:704-745
Modified: branches/ocean_projects/timesplitting_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-05-31 17:14:20 UTC (rev 863)
+++ branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-05-31 19:41:50 UTC (rev 864)
@@ -40,6 +40,7 @@
config_implicit_vertical_mix = .false.
config_convective_visc = 1.0
config_convective_diff = 1.0
+ config_bottom_drag_coeff = 1.0e-3
/
&vmix_const
config_vert_visc = 2.5e-5
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-31 17:14:20 UTC (rev 863)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-31 19:41:50 UTC (rev 864)
@@ -33,6 +33,7 @@
namelist logical vmix config_implicit_vertical_mix .true.
namelist real vmix config_convective_visc 1.0
namelist real vmix config_convective_diff 1.0
+namelist real vmix config_bottom_drag_coeff 1.0e-3
namelist real vmix_const config_vert_visc 2.5e-4
namelist real vmix_const config_vert_diff 2.5e-5
namelist real vmix_rich config_bkrd_vert_visc 1.0e-4
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-31 17:14:20 UTC (rev 863)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-31 19:41:50 UTC (rev 864)
@@ -82,7 +82,7 @@
integer :: nCells, nEdges, nVertLevels, num_tracers
real (kind=RKIND), dimension(:,:), pointer :: &
- u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell
+ u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell, ke_edge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
@@ -249,6 +249,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
@@ -293,7 +294,8 @@
/ (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
/ h_edge(k,iEdge)
enddo
- A(maxLevelEdgeTop(iEdge)) = 0.0
+ 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)
@@ -1559,9 +1561,11 @@
! -c |u| u where c is unitless and 1.0e-3.
! see POP Reference guide, section 3.4.4.
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - 1.0e-3*u(k,iEdge) &
- *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+ if (.not.config_implicit_vertical_mix) then
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - config_bottom_drag_coeff*u(k,iEdge) &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+ end if
endif
@@ -2356,7 +2360,6 @@
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: r, h1, h2
-
h => s % h % array
u => s % u % array
v => s % v % array
@@ -2715,12 +2718,11 @@
! equation of state
!
! For an isopycnal model, density should remain constant.
+ ! For zlevel, calculate in-situ density
if (config_vert_grid_type.eq.'zlevel') then
- call equation_of_state(s,grid,-1)
-
- ! mrp 110324 In order to vissualize rhoDisplaced, include the following
- !call equation_of_state(s, grid, 1)
-
+ 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
!
@@ -2900,7 +2902,7 @@
real (kind=RKIND) :: coef
real (kind=RKIND), dimension(:,:), pointer :: &
vertViscTopOfEdge, vertDiffTopOfCell, &
- RiTopOfEdge, RiTopOfCell, rhoDisplaced, &
+ RiTopOfEdge, RiTopOfCell, rhoDisplaced, rho, &
kiteAreasOnVertex, h_edge, h, u
real (kind=RKIND), dimension(:), pointer :: &
@@ -2918,6 +2920,7 @@
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: r, h1, h2
+ rho => s % rho % array
rhoDisplaced => s % rhoDisplaced % array
u => s % u % array
h => s % h % array
@@ -2951,15 +2954,17 @@
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^*$,
+ ! compute density of parcel displaced to next deeper z-level,
! in state % rhoDisplaced
- call equation_of_state(s, grid, 1)
+!maltrud make sure rho is current--check this for redundancy
+ call equation_of_state(s, grid, 0, 'relative')
+ call equation_of_state(s, grid, 1, 'relative')
! 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)
+ drhoTopOfCell(k,iCell) = rho(k-1,iCell) - rhoDisplaced(k-1,iCell)
end do
end do
@@ -3060,6 +3065,16 @@
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
@@ -3119,6 +3134,16 @@
+ (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
@@ -3189,7 +3214,7 @@
end subroutine enforce_boundaryEdge
- subroutine equation_of_state(s, grid, k_displaced)
+ 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.
@@ -3209,6 +3234,7 @@
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
@@ -3234,7 +3260,7 @@
elseif (config_eos_type.eq.'jm') then
- call equation_of_state_jm(s, grid,k_displaced)
+ call equation_of_state_jm(s, grid, k_displaced, displacement_type)
else
print *, ' Incorrect choice of config_eos_type:',&
@@ -3245,7 +3271,7 @@
end subroutine equation_of_state
- subroutine equation_of_state_jm(s, grid,k_displaced)
+ 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.
@@ -3272,8 +3298,8 @@
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
@@ -3375,6 +3401,8 @@
bup2s1t1 = 6.128773e-8, &
bup2s1t2 = 6.207323e-10
+ integer :: k_test, k_ref
+
tracers => s % tracers % array
nCells = grid % nCells
@@ -3405,26 +3433,43 @@
+ 0.100766*depth + 2.28405e-7*depth**2
enddo
- ! 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.
- if (k_displaced.le.0) then
+ ! 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
- elseif (k_displaced.le.nVertLevels) then
+ else ! k_displaced /= 0
rhoPointer => s % rhoDisplaced % array
do k=1,nVertLevels
- p(k) = pRefEOS(k_displaced)
+ 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
- else
- write(0,*) 'Abort: In equation_of_state_jm', &
- ' k_displaced must not be larger than nVertLevels'
- call dmpar_abort(dminfo)
endif
do iCell=1,nCells
</font>
</pre>