<p><b>maltrud@lanl.gov</b> 2011-05-25 14:48:44 -0600 (Wed, 25 May 2011)</p><p>added bottom drag as implicit when vertical mixing is implicit.<br>
<br>
added bottom drag coefficient to the vmix namelist.<br>
<br>
changed Richardson Number calculation to use nearest z-level as the stability<br>
reference instead of the surface.<br>
<br>
added limiting on the Richardson Number-based mixing coefficients<br>
to be <= the convective value. theoretically, this should not be<br>
necessary, but may be causing some buggy behaviour.<br>
<br>
modified equation_of_state to accept an extra flag which specifies whether<br>
potential density is calculated referenced to a relative or absolute model level.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry 2011-05-24 19:31:25 UTC (rev 851)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry 2011-05-25 20:48:44 UTC (rev 852)
@@ -26,6 +26,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/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-05-24 19:31:25 UTC (rev 851)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-05-25 20:48:44 UTC (rev 852)
@@ -78,7 +78,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
@@ -244,6 +244,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
@@ -288,7 +289,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)
@@ -731,9 +733,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
@@ -1394,7 +1398,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
@@ -1737,12 +1740,11 @@
! equation of state
!
! For an isopycnal model, density should remain constant.
+ ! For zlevel, calculate in-istu 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
!
@@ -1867,7 +1869,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 :: &
@@ -1885,6 +1887,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
@@ -1918,15 +1921,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
@@ -2027,6 +2032,14 @@
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
+ if (vertViscTopOfEdge(k,iEdge) > config_convective_visc &
+ .and. 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
else
! mrp 110324 efficiency note: this if is inside iCell and k loops.
if (config_implicit_vertical_mix) then
@@ -2086,6 +2099,14 @@
+ (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
+ if (vertDiffTopOfCell(k,iCell) > config_convective_diff &
+ .and. 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
else
! mrp 110324 efficiency note: this if is inside iCell and k loops.
if (config_implicit_vertical_mix) then
@@ -2156,7 +2177,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.
@@ -2176,6 +2197,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
@@ -2201,7 +2223,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:',&
@@ -2212,7 +2234,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.
@@ -2239,8 +2261,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
@@ -2342,6 +2364,8 @@
bup2s1t1 = 6.128773e-8, &
bup2s1t2 = 6.207323e-10
+ integer :: k_test, k_ref
+
tracers => s % tracers % array
nCells = grid % nCells
@@ -2372,26 +2396,42 @@
+ 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>