<p><b>ringler@lanl.gov</b> 2011-06-16 12:54:45 -0600 (Thu, 16 Jun 2011)</p><p><br>
completed addition and testing of bottom drag<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/namelist.input.pvsw
===================================================================
--- branches/pv_based_swm/mpas/namelist.input.pvsw 2011-06-10 22:25:01 UTC (rev 899)
+++ branches/pv_based_swm/mpas/namelist.input.pvsw 2011-06-16 18:54:45 UTC (rev 900)
@@ -5,19 +5,22 @@
config_test_case = 0
config_time_integration = 'RK4'
config_dt = 2400
- config_ntimesteps = 3600
- config_output_interval = 36
+ config_ntimesteps = 36
+ config_output_interval = 1
config_stats_interval = 0
config_h_mom_eddy_visc2 = 0.0
config_h_mom_eddy_visc4 = 0.0
- config_h_tracer_eddy_diff2 = 0.0
+ config_h_tracer_eddy_diff2 = 1260.0
config_h_tracer_eddy_diff4 = 0.0
config_thickness_adv_order = 2
- config_tracer_adv_order = 2
+ config_tracer_adv_order = 4
config_positive_definite = .false.
config_monotonic = .false.
config_wind_stress = .false.
config_bottom_drag = .false.
+ config_bottom_drag_term1 = .true.
+ config_bottom_drag_term2 = .true.
+ config_dragCoefficient = 2.0e-5
/
&io
Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-06-10 22:25:01 UTC (rev 899)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-06-16 18:54:45 UTC (rev 900)
@@ -18,6 +18,8 @@
namelist logical sw_model config_monotonic false
namelist logical sw_model config_wind_stress false
namelist logical sw_model config_bottom_drag false
+namelist logical sw_model config_bottom_drag_term1 false
+namelist logical sw_model config_bottom_drag_term2 false
namelist real sw_model config_apvm_upwinding 0.5
namelist real sw_model config_alpha 0.0
namelist logical sw_model config_alpha_closure false
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-06-10 22:25:01 UTC (rev 899)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-06-16 18:54:45 UTC (rev 900)
@@ -72,11 +72,13 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
- call compute_drag(block % state % time_levs(1) % state, block % mesh)
-
+
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
+
+ call reconstructGradKE(block % state % time_levs(1) % state, mesh)
+ call compute_drag(block % state % time_levs(1) % state, mesh)
if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-06-10 22:25:01 UTC (rev 899)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-06-16 18:54:45 UTC (rev 900)
@@ -7,7 +7,7 @@
implicit none
- real (KIND=RKIND), parameter :: threshold = 1.0e-10
+ real (KIND=RKIND), parameter :: threshold = 1.0e-8
integer, parameter :: nIterPerGS = 2, nCheckPoisson = 25, nQuitPoisson = 5000
private
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-06-10 22:25:01 UTC (rev 899)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-06-16 18:54:45 UTC (rev 900)
@@ -101,13 +101,13 @@
block % state % time_levs(2) % state % pv_forcing % array(:,:) = 0.0
block % state % time_levs(2) % state % pv_drag % array(:,:) = 0.0
- block % state % time_levs(2) % state % uReconstructX % array(:,:) = block % state % time_levs(1) % state % uReconstructX % array(:,:)
- block % state % time_levs(2) % state % uReconstructY % array(:,:) = block % state % time_levs(1) % state % uReconstructY % array(:,:)
- block % state % time_levs(2) % state % uReconstructZ % array(:,:) = block % state % time_levs(1) % state % uReconstructZ % array(:,:)
+ ! block % state % time_levs(2) % state % uReconstructX % array(:,:) = block % state % time_levs(1) % state % uReconstructX % array(:,:)
+ ! block % state % time_levs(2) % state % uReconstructY % array(:,:) = block % state % time_levs(1) % state % uReconstructY % array(:,:)
+ ! block % state % time_levs(2) % state % uReconstructZ % array(:,:) = block % state % time_levs(1) % state % uReconstructZ % array(:,:)
- block % state % time_levs(2) % state % gradKeReconstructX % array(:,:) = block % state % time_levs(1) % state % gradKeReconstructX % array(:,:)
- block % state % time_levs(2) % state % gradKeReconstructY % array(:,:) = block % state % time_levs(1) % state % gradKeReconstructY % array(:,:)
- block % state % time_levs(2) % state % gradKeReconstructZ % array(:,:) = block % state % time_levs(1) % state % gradKeReconstructZ % array(:,:)
+ ! block % state % time_levs(2) % state % gradKeReconstructX % array(:,:) = block % state % time_levs(1) % state % gradKeReconstructX % array(:,:)
+ ! block % state % time_levs(2) % state % gradKeReconstructY % array(:,:) = block % state % time_levs(1) % state % gradKeReconstructY % array(:,:)
+ ! block % state % time_levs(2) % state % gradKeReconstructZ % array(:,:) = block % state % time_levs(1) % state % gradKeReconstructZ % array(:,:)
do iCell=1,block % mesh % nCells ! couple tracers to h
@@ -239,9 +239,6 @@
call reconstructGradKE(provis, block % mesh)
call compute_drag(provis, block % mesh)
- write(6,*) ' uR min', minval(provis % uReconstructX % array(:,:))
- write(6,*) ' uR max ', maxval(provis % uReconstructX % array(:,:))
-
block => block % next
end do
end if
@@ -330,7 +327,6 @@
call reconstruct(block % state % time_levs(2) % state, block % mesh)
call reconstructGradKE(block % state % time_levs(2) % state, block % mesh)
call compute_drag(block % state % time_levs(2) % state, block % mesh)
- write(6,*) ' uR minc', minval(block % state % time_levs(2) % state % uReconstructX % array(:,:))
block => block % next
end do
@@ -469,6 +465,14 @@
u => s % u % array
h_edge => s % h_edge % array
+ pv_drag => s % pv_drag % array
+
+ tracer_tend => tend % tracers % array
+ tend_pv_advection => tend % pv_advection % array
+ tend_pv_dissipation => tend % pv_dissipation % array
+ tend_pv_forcing => tend % pv_forcing % array
+ tend_pv_drag => tend % pv_drag % array
+
dcEdge => grid % dcEdge % array
deriv_two => grid % deriv_two % array
dvEdge => grid % dvEdge % array
@@ -477,12 +481,6 @@
boundaryCell=> grid % boundaryCell % array
boundaryEdge=> grid % boundaryEdge % array
areaCell => grid % areaCell % array
- tracer_tend => tend % tracers % array
- pv_drag => s % pv_drag % array
- tend_pv_advection => tend % pv_advection % array
- tend_pv_dissipation => tend % pv_dissipation % array
- tend_pv_forcing => tend % pv_forcing % array
- tend_pv_drag => tend % pv_drag % array
coef_3rd_order = 0.
if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
@@ -801,27 +799,49 @@
pv_drag(:,:) = 0.0
- ! do iCell = 1,grid%nCells
- ! do k=1,grid % nVertLevels
- ! pv_drag(k,iCell) = -config_dragCoefficient * sqrt(2.0*ke(k,iCell)) * vor(k,iCell) / h(k,iCell)
- ! enddo
- ! enddo
+ if(config_bottom_drag) then
+ if(config_bottom_drag_term1) then
+
+ ! term1: -C_d * |u| * relative PV
+ ! the drag has to be spread throughout the column, so divide by thickness
+ ! units of drag have to be equivalent to d/dt (h*q) or 1/s * m/s * 1/(ms) == 1/s^2
do iCell = 1,grid%nCells
do k=1,grid % nVertLevels
+ pv_drag(k,iCell) = -config_dragCoefficient * sqrt(2.0*ke(k,iCell)) * vor(k,iCell) / h(k,iCell)
+ enddo
+ enddo
+
+ endif
+
+ if(config_bottom_drag_term2) then
+
+ ! term2: -C_d * grad |u| cross u
+ ! the drag has to be spread throughout the column, so divide by thickness
+ ! units of drag have to be equivalent to d/dt (h*q) or 1/s * m/s * 1/(ms) == 1/s^2
+ do iCell = 1,grid%nCells
+ do k=1,grid % nVertLevels
p(1) = gx(k,iCell); p(2) = gy(k,iCell); p(3) = gz(k,iCell)
q(1) = ux(k,iCell); q(2) = uy(k,iCell); q(3) = uz(k,iCell)
call cross_product_in_R3(p,q,t)
- r = xCell(iCell)*t(1) + yCell(iCell)*t(2) + zCell(iCell)*t(3)
+ if(grid % on_a_sphere) then
+ r = xCell(iCell)*t(1) + yCell(iCell)*t(2) + zCell(iCell)*t(3)
+ else
+ r = t(3)
+ endif
pv_drag(k,iCell) = pv_drag(k,iCell) - config_dragCoefficient * r / h(k,iCell)
- write(6,*) 'p ',p(:)
- write(6,*) 'q ',q(:)
- write(6,*) 't ',t(:)
- write(6,*) 'r ',r
- write(6,*) 'drag ',config_dragCoefficient * r / h(k,iCell)
+ ! write(6,*) 'p ',p(:)
+ ! write(6,*) 'q ',q(:)
+ ! write(6,*) 't ',t(:)
+ ! write(6,*) 'r ',r
+ ! write(6,*) 'drag ',vor(k,iCell),-config_dragCoefficient * r / h(k,iCell)
enddo
enddo
+ endif
+
+ endif
+
contains
subroutine cross_product_in_R3(p_1,p_2,p_out)
@@ -957,13 +977,13 @@
!
! Compute kinetic energy at every edge
- !
+ ! (actually this is the gradient of |u|)
gradKE(:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,nVertLevels
- gradKE(k,iEdge) = (ke(k,cell2) - ke(k,cell1)) / dcEdge(iEdge)
+ gradKE(k,iEdge) = (sqrt(2.0*ke(k,cell2)) - sqrt(2.0*ke(k,cell1)) ) / dcEdge(iEdge)
enddo
enddo
Modified: branches/pv_based_swm/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/pv_based_swm/mpas/src/operators/module_vector_reconstruction.F 2011-06-10 22:25:01 UTC (rev 899)
+++ branches/pv_based_swm/mpas/src/operators/module_vector_reconstruction.F 2011-06-16 18:54:45 UTC (rev 900)
@@ -238,7 +238,7 @@
edgesOnCell => grid % edgesOnCell % array
nEdgesOnCell=> grid % nEdgesOnCell % array
nCellsSolve = grid % nCellsSolve
- u => state % u % array
+ u => state % gradKE % array
uReconstructX => state % gradKeReconstructX % array
uReconstructY => state % gradKeReconstructY % array
uReconstructZ => state % gradKeReconstructZ % array
</font>
</pre>