<p><b>mpetersen@lanl.gov</b> 2011-02-07 16:04:29 -0700 (Mon, 07 Feb 2011)</p><p>BRANCH COMMIT: A few more updates. Add config_vert_tracer_adv_order to Registry and namelist.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/ocean_projects/vert_adv_mrp
___________________________________________________________________
Added: svn:mergeinfo
+ /trunk/mpas:703-714
Modified: branches/ocean_projects/vert_adv_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/vert_adv_mrp/namelist.input.ocean 2011-02-07 22:23:18 UTC (rev 715)
+++ branches/ocean_projects/vert_adv_mrp/namelist.input.ocean 2011-02-07 23:04:29 UTC (rev 716)
@@ -45,7 +45,8 @@
config_eos_type = 'linear'
/
&advection
- config_vert_tracer_adv = 'upwind'
+ config_vert_tracer_adv = 'spline'
+ config_vert_tracer_adv_order = 3
config_tracer_adv_order = 2
config_thickness_adv_order = 2
config_positive_definite = .false.
Modified: branches/ocean_projects/vert_adv_mrp/namelist.input.sw
===================================================================
--- branches/ocean_projects/vert_adv_mrp/namelist.input.sw 2011-02-07 22:23:18 UTC (rev 715)
+++ branches/ocean_projects/vert_adv_mrp/namelist.input.sw 2011-02-07 23:04:29 UTC (rev 716)
@@ -13,6 +13,8 @@
config_tracer_adv_order = 2
config_positive_definite = .false.
config_monotonic = .false.
+ config_wind_stress = .false.
+ config_bottom_drag = .false.
/
&io
Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/Registry 2011-02-07 22:23:18 UTC (rev 715)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/Registry 2011-02-07 23:04:29 UTC (rev 716)
@@ -32,7 +32,8 @@
namelist real vmix config_vmixTanhZMid -100
namelist real vmix config_vmixTanhZWidth 100
namelist character eos config_eos_type linear
-namelist character advection config_vert_tracer_adv 'centered'
+namelist character advection config_vert_tracer_adv stencil
+namelist integer advection config_vert_tracer_adv_order 4
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
namelist logical advection config_positive_definite false
Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-07 22:23:18 UTC (rev 715)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-07 23:04:29 UTC (rev 716)
@@ -897,8 +897,11 @@
allocate(tracerTop(num_tracers,nVertLevels+1))
tracerTop(:,1)=0.0
- if (config_vert_tracer_adv.eq.'centered') then
+ if (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.2) then
+ ! Compute tracerTop using centered stencil, a simple average.
+
do iCell=1,grid % nCellsSolve
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
@@ -917,54 +920,12 @@
end do
end do
- elseif (config_vert_tracer_adv.eq.'upwind') then
+ elseif (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.3) then
- do iCell=1,grid % nCellsSolve
- do k=2,maxLevelCell(iCell)
- if (wTop(k,iCell)>0.0) then
- upwindCell = k
- else
- upwindCell = k-1
- endif
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
- end do
- end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
- end do
+ ! Compute tracerTop using 3rd order stencil. This is the same
+ ! as 4th order, but includes upwinding.
- elseif (config_vert_tracer_adv.eq.'linear') then
-
- do iCell=1,grid % nCellsSolve
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! Linear interpolation of tracer value at Top interface.
- ! Note hRatio on the k side is multiplied by tracer at k-1
- ! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
- tracerTop(iTracer,k) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
- end do
- end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
- end do
-
- elseif (config_vert_tracer_adv.eq.'flux3') then
-
! Hardwire flux3Coeff at 1.0 for now. Could add this to the
! namelist, if desired.
flux3Coef = 1.0
@@ -1002,8 +963,11 @@
end do
end do
- elseif (config_vert_tracer_adv.eq.'flux4') then
+ elseif (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.4) then
+ ! Compute tracerTop using 4rd order stencil [-1 7 7 -1]
+
do iCell=1,grid % nCellsSolve
k=2
do iTracer=1,num_tracers
@@ -1037,23 +1001,54 @@
end do
end do
- elseif (config_vert_tracer_adv.eq.'CubicSpline') then
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.2) then
+ ! Compute tracerTop using linear interpolation.
+
+ do iCell=1,grid % nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! Note hRatio on the k side is multiplied by tracer at k-1
+ ! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
+ tracerTop(iTracer,k) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ end do
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ end do
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.3) then
+
+ ! Compute tracerTop using cubic spline interpolation.
+
allocate(tracer2ndDer(nVertLevels))
allocate(tracersIn(nVertLevels),tracersOut(nVertLevels), &
- posZMidZLevel(nVertLevels), posZTopZLevel(nVertLevels))
+ posZMidZLevel(nVertLevels), posZTopZLevel(nVertLevels-1))
+ ! For the ocean, zlevel coordinates are negative and decreasing,
+ ! but spline functions assume increasing, so flip to positive.
+
+ posZMidZLevel = -zMidZLevel(1:nVertLevels)
+ posZTopZLevel = -zTopZLevel(2:nVertLevels)
+
do iCell=1,grid % nCellsSolve
! mrp 110201 efficiency note: push tracer loop down
! into spline subroutines to improve efficiency
do iTracer=1,num_tracers
! Place data in arrays to avoid creating new temporary arrays for every
- ! subroutine call. For the ocean, zlevel coordinates are negative and
- ! decreasing, but spline functions assume increasing, so flip to positive
+ ! subroutine call.
tracersIn(1:maxLevelCell(iCell))=tracers(iTracer,1:maxLevelCell(iCell),iCell)
- posZMidZLevel(1:maxLevelCell(iCell)) = -zMidZLevel(1:maxLevelCell(iCell))
- posZTopZLevel(1:maxLevelCell(iCell)-1) = -zTopZLevel(2:maxLevelCell(iCell))
call CubicSplineCoefficients(posZMidZLevel, &
tracersIn, &
@@ -1079,8 +1074,18 @@
deallocate(tracer2ndDer)
deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
- endif
+ else
+
+ print *, 'Abort: Incorrect combination of ', &
+ 'config_vert_tracer_adv and config_vert_tracer_adv_order.'
+ print *, 'Use:'
+ print *, 'config_vert_tracer_adv=''stencil'' and_vert_tracer_adv_order=2,3,4 or'
+ print *, 'config_vert_tracer_adv=''spline'' and_vert_tracer_adv_order=2,3'
+ call dmpar_abort(dminfo)
+
+ endif ! vertical tracer advection method
+
deallocate(tracerTop)
endif ! ZLevel
Modified: branches/ocean_projects/vert_adv_mrp/src/core_sw/Registry
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_sw/Registry 2011-02-07 22:23:18 UTC (rev 715)
+++ branches/ocean_projects/vert_adv_mrp/src/core_sw/Registry 2011-02-07 23:04:29 UTC (rev 716)
@@ -15,6 +15,8 @@
namelist integer sw_model config_tracer_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic false
+namelist logical sw_model config_wind_stress false
+namelist logical sw_model config_bottom_drag false
namelist real sw_model config_apvm_upwinding 0.5
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -112,6 +114,7 @@
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
+var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
Modified: branches/ocean_projects/vert_adv_mrp/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_sw/module_time_integration.F 2011-02-07 22:23:18 UTC (rev 715)
+++ branches/ocean_projects/vert_adv_mrp/src/core_sw/module_time_integration.F 2011-02-07 23:04:29 UTC (rev 716)
@@ -269,8 +269,10 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+ real (kind=RKIND) :: ke_edge
-
h => s % h % array
u => s % u % array
v => s % v % array
@@ -308,6 +310,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ u_src => grid % u_src % array
!
! Compute height tendency for each cell
@@ -497,6 +500,28 @@
end if
+ ! Compute u (velocity) tendency from wind stress (u_src)
+ if(config_wind_stress) then
+ do iEdge=1,grid % nEdges
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ end do
+ endif
+
+ if (config_bottom_drag) then
+ do iEdge=1,grid % nEdges
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
+ ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
+ + ke(1,cellsOnEdge(2,iEdge)))
+
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ - 1.0e-3*u(1,iEdge) &
+ *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
+ end do
+ endif
+
end subroutine compute_tend
</font>
</pre>