<p><b>ringler@lanl.gov</b> 2011-02-10 14:39:49 -0700 (Thu, 10 Feb 2011)</p><p><br>
add a branch that contains the GM closure.<br>
<br>
The closure is only applicable to isopycnal coordinates.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/GMClosure/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2011-02-09 22:20:51 UTC (rev 735)
+++ branches/ocean_projects/GMClosure/namelist.input.ocean 2011-02-10 21:39:49 UTC (rev 736)
@@ -1,10 +1,12 @@
&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 = 45.0
+ config_ntimesteps = 23040000
+ config_output_interval = 7680
+ config_stats_interval = 7680
+ config_uTPV = .false.
+ config_diffPV = .false.
/
&io
@@ -21,12 +23,14 @@
&grid
config_vert_grid_type = 'isopycnal'
- config_rho0 = 1015
+ config_rho0 = 1013
/
&hmix
+ config_h_kappa = 0.0
+ config_h_kappa_q = 0.0
config_h_mom_eddy_visc2 = 0.0
config_h_mom_eddy_visc4 = 5.0e8
- config_h_tracer_eddy_diff2 = 10.0
+ config_h_tracer_eddy_diff2 = 0.0
config_h_tracer_eddy_diff4 = 0.0
/
&vmix
Modified: branches/ocean_projects/GMClosure/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2011-02-09 22:20:51 UTC (rev 735)
+++ branches/ocean_projects/GMClosure/src/core_ocean/Registry 2011-02-10 21:39:49 UTC (rev 736)
@@ -1,21 +1,24 @@
#
# namelist type namelist_record name default_value
#
+namelist logical sw_model config_uTPV false
+namelist logical sw_model config_diffPV false
namelist integer sw_model config_test_case 5
namelist character sw_model config_time_integration RK4
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist integer sw_model config_stats_interval 100
-namelist character io config_input_name grid.nc
-namelist character io config_output_name output.nc
-namelist character io config_restart_name restart.nc
-namelist character io config_decomp_file_prefix graph.info.part.
+namelist character io config_input_name grid.nc
+namelist character io config_output_name output.nc
+namelist character io config_restart_name restart.nc
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
namelist character grid config_vert_grid_type isopycnal
namelist real grid config_rho0 1028
+namelist real hmix config_h_kappa 0.0
+namelist real hmix config_h_kappa_q 0.0
namelist real hmix config_h_mom_eddy_visc2 0.0
namelist real hmix config_h_mom_eddy_visc4 0.0
namelist real hmix config_h_tracer_eddy_diff2 0.0
@@ -155,6 +158,13 @@
var persistent real tend_tracer2 ( nVertLevels nCells Time ) 1 - tracer2 tend tracers testing
# Diagnostic fields: only written to output
+var persistent real uBolus ( nVertLevels nEdges Time ) 2 o uBolus state - -
+var persistent real uTransport ( nVertLevels nEdges Time ) 2 o uTransport state - -
+var persistent real uBolusX ( nVertLevels nCells Time ) 2 o uBolusX state - -
+var persistent real uBolusY ( nVertLevels nCells Time ) 2 o uBolusY state - -
+var persistent real uBolusZ ( nVertLevels nCells Time ) 2 o uBolusZ state - -
+var persistent real uBolusZonal ( nVertLevels nCells Time ) 2 o uBolusZonal state - -
+var persistent real uBolusMeridional ( nVertLevels nCells Time ) 2 o uBolusMeridional state - -
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
var persistent real vorticity ( nVertLevels nVertices Time ) 2 - vorticity state - -
Modified: branches/ocean_projects/GMClosure/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2011-02-09 22:20:51 UTC (rev 735)
+++ branches/ocean_projects/GMClosure/src/core_ocean/module_time_integration.F 2011-02-10 21:39:49 UTC (rev 736)
@@ -125,6 +125,12 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % uBolus % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % uTransport % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
if (config_h_mom_eddy_visc4 > 0.0) then
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
@@ -275,8 +281,8 @@
zMidZLevel, zTopZLevel
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
+ tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, pv_vertex, &
+ MontPot, wTop, divergence, uTransport
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
@@ -297,6 +303,7 @@
h => s % h % array
u => s % u % array
+ uTransport => s % uTransport % array
v => s % v % array
wTop => s % wTop % array
h_edge => s % h_edge % array
@@ -306,6 +313,7 @@
ke => s % ke % array
ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
+ pv_vertex => s % pv_vertex % array
MontPot => s % MontPot % array
pressure => s % pressure % array
@@ -366,7 +374,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,kmax
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) - flux
tend_h(k,cell2) = tend_h(k,cell2) + flux
end do
@@ -451,7 +459,7 @@
! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity )
! strictly only valid for config_h_mom_eddy_visc2 == constant
!
- if ( config_h_mom_eddy_visc2 > 0.0 ) then
+ if ( config_h_mom_eddy_visc2 > 0.0 .or. config_diffPV ) then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -464,10 +472,18 @@
! is - </font>
<font color="black">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
! + k \times </font>
<font color="gray">abla vorticity pointing from cell1 to cell2.
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
+ if(config_diffPV) then
+ u_diffusion = -( pv_vertex(k,vertex2) - pv_vertex(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = config_h_kappa_q * h_edge(k,iEdge) * u_diffusion
+
+ else
+
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
+ endif
+
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
end do
@@ -582,7 +598,11 @@
do j = 1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+ if(config_uTPV) then
+ q = q + weightsOnEdge(j,iEdge) * uTransport(k,eoe) * workpv * h_edge(k,eoe)
+ else
+ q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+ endif
end do
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ q &
@@ -698,7 +718,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, uTransport
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:,:), pointer :: boundaryEdge
@@ -719,6 +739,7 @@
u => s % u % array
+ uTransport => s % uTransport % array
h => s % h % array
boundaryCell=> grid % boundaryCell % array
wTop => s % wTop % array
@@ -769,7 +790,7 @@
do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
end do
@@ -811,13 +832,13 @@
!-- if u > 0:
if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ flux = dvEdge(iEdge) * uTransport(k,iEdge) * h_edge(k,iEdge) * ( &
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
!-- else u <= 0:
else
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ flux = dvEdge(iEdge) * uTransport(k,iEdge) * h_edge(k,iEdge) * ( &
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
@@ -863,7 +884,7 @@
endif
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ flux = dvEdge(iEdge) * uTransport(k,iEdge) * h_edge(k,iEdge) * ( &
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
@@ -1123,7 +1144,7 @@
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- rho, temperature, salinity
+ rho, temperature, salinity, uBolus, uTransport
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:,:), allocatable:: div_u
@@ -1143,6 +1164,8 @@
h => s % h % array
u => s % u % array
+ uBolus => s % uBolus % array
+ uTransport => s % uTransport % array
v => s % v % array
wTop => s % wTop % array
h_edge => s % h_edge % array
@@ -1309,6 +1332,7 @@
! used -1e34 so error clearly occurs if these values are used.
!
u(:,nEdges+1) = -1e34
+ uTransport(:,nEdges+1) = -1e34
h(:,nCells+1) = -1e34
tracers(s % index_temperature,:,nCells+1) = -1e34
tracers(s % index_salinity,:,nCells+1) = -1e34
@@ -1351,7 +1375,21 @@
enddo
!
+ !
+ ! Compute Bolus velocity as bolus == - /kappa (/grad h) / h
+ !
+ uBolus(:,:) = 0.0
+ do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ uBolus(k,iEdge) = -config_h_kappa * (h(k,cell2) - h(k,cell1))/dcEdge(iEdge)/h_edge(k,iEdge)
+ end do
+ end do
+
+ !
! Compute kinetic energy in each cell
+ ! NOTE: This includes the Bolus velocity, should reduce to RTKS in the limit of Bolus == 0
!
ke(:,:) = 0.0
do iCell=1,nCells
@@ -1359,6 +1397,7 @@
iEdge = edgesOnCell(i,iCell)
do k=1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ if(config_uTPV) ke(k,iCell) = ke(k,iCell) + 0.50 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)*uBolus(k,iEdge)
end do
end do
do k=1,maxLevelCell(iCell)
@@ -1367,6 +1406,12 @@
end do
!
+ !
+ ! Compute the transport velocity (this velocity transports thickness and PV)
+ !
+ uTransport(:,:) = u(:,:) + uBolus(:,:)
+
+ !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
</font>
</pre>