<p><b>duda</b> 2009-11-23 10:31:54 -0700 (Mon, 23 Nov 2009)</p><p>Add horizontal divergence as a diagnostic field from repository trunk.<br>
Add del2 mixing for velocity from repository trunk.<br>
Add namelist options to control eddy viscosities for theta and <br>
velocity in the horizontal and vertical.<br>
<br>
M src/module_time_integration.F<br>
M src/module_constants.F<br>
M Registry/Registry<br>
M namelist.input<br>
</p><hr noshade><pre><font color="gray">Modified: branches/hyd_model/Registry/Registry
===================================================================
--- branches/hyd_model/Registry/Registry 2009-11-20 00:50:42 UTC (rev 78)
+++ branches/hyd_model/Registry/Registry 2009-11-23 17:31:54 UTC (rev 79)
@@ -6,6 +6,10 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
+namelist real sw_model config_h_mom_eddy_visc2 0.0
+namelist real sw_model config_v_mom_eddy_visc2 0.0
+namelist real sw_model config_h_theta_eddy_visc2 0.0
+namelist real sw_model config_v_theta_eddy_visc2 0.0
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
@@ -98,6 +102,7 @@
# Diagnostic fields: only written to output
var real v ( nVertLevels nEdges Time ) o v
+var real divergence ( nVertLevels nCells Time ) o divergence
var real vorticity ( nVertLevels nVertices Time ) o vorticity
var real pv_edge ( nVertLevels nEdges Time ) o pv_edge
var real h_edge ( nVertLevels nEdges Time ) o h_edge
Modified: branches/hyd_model/namelist.input
===================================================================
--- branches/hyd_model/namelist.input 2009-11-20 00:50:42 UTC (rev 78)
+++ branches/hyd_model/namelist.input 2009-11-23 17:31:54 UTC (rev 79)
@@ -5,6 +5,10 @@
config_ntimesteps = 240
config_output_interval = 24
config_number_of_sub_steps = 4
+ config_h_mom_eddy_visc2 = 0.0
+ config_v_mom_eddy_visc2 = 0.0
+ config_h_theta_eddy_visc2 = 0.0
+ config_v_theta_eddy_visc2 = 0.0
/
&restart
Modified: branches/hyd_model/src/module_constants.F
===================================================================
--- branches/hyd_model/src/module_constants.F 2009-11-20 00:50:42 UTC (rev 78)
+++ branches/hyd_model/src/module_constants.F 2009-11-23 17:31:54 UTC (rev 79)
@@ -8,7 +8,7 @@
real (kind=RKIND), parameter :: cp = 1003.
real (kind=RKIND), parameter :: cv = 716. ! cp - rgas
real (kind=RKIND), parameter :: cvpm = -.71385842 ! -cv/cp
- real (kind=RKIND), parameter :: prandtl = 3.
+ real (kind=RKIND), parameter :: prandtl = 1.0
contains
Modified: branches/hyd_model/src/module_time_integration.F
===================================================================
--- branches/hyd_model/src/module_time_integration.F 2009-11-20 00:50:42 UTC (rev 78)
+++ branches/hyd_model/src/module_time_integration.F 2009-11-23 17:31:54 UTC (rev 79)
@@ -491,12 +491,13 @@
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
integer :: nCells, nEdges, nVertices, nVertLevels
-! real (kind=RKIND), pointer :: vnu, hnu
- real (kind=RKIND) :: vnu, hnu
+ real (kind=RKIND) :: h_mom_eddy_visc2, v_mom_eddy_visc2
+ real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2
+ real (kind=RKIND) :: u_diffusion
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, p_s
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, geopotential, theta, ww, h_diabatic, &
- tend_theta
+ circulation, divergence, vorticity, ke, pv_edge, geopotential, theta, ww, &
+ h_diabatic, tend_theta
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -511,6 +512,8 @@
u => s % u % array
h_edge => s % h_edge % array
circulation => s % circulation % array
+ divergence => s % divergence % array
+ vorticity => s % vorticity % array
ke => s % ke % array
pv_edge => s % pv_edge % array
geopotential => s % geopotential % array
@@ -542,11 +545,12 @@
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
-! vnu => grid % vnu
-! hnu => grid % hnu
- hnu = 0.
- vnu = 0.
+ h_mom_eddy_visc2 = config_h_mom_eddy_visc2
+ v_mom_eddy_visc2 = config_v_mom_eddy_visc2
+ h_theta_eddy_visc2 = config_h_theta_eddy_visc2
+ v_theta_eddy_visc2 = config_v_theta_eddy_visc2
+
!
! Compute u (normal) velocity tendency for each edge (cell face)
!
@@ -557,6 +561,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
do k=1,nVertLevels
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
@@ -601,58 +606,86 @@
end do
#endif
+
!
+ ! horizontal mixing for u
+ !
+ if ( h_mom_eddy_visc2 > 0.0 ) then
+
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
+ ! only valid for h_mom_eddy_visc2 == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = h_mom_eddy_visc2 * u_diffusion
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
+
+ end if
+
+
+ !
! vertical advection for u
!
- do iEdge=1,grid % nEdges
+ do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
wdun(1) = 0.
- if (cell1 > 0 .and. cell2 > 0) then
do k=2,nVertLevels
- wdun(k) = &
- (ww(k,cell1)/(h(k,cell1)+h(k-1,cell1)) +ww(k,cell2)/(h(k,cell2)+h(k-1,cell2)))* &
- rdnu(k)*(u(k,iEdge)-u(k-1,iEdge))
+ wdun(k) = &
+ (ww(k,cell1)/(h(k,cell1)+h(k-1,cell1)) +ww(k,cell2)/(h(k,cell2)+h(k-1,cell2)))* &
+ rdnu(k)*(u(k,iEdge)-u(k-1,iEdge))
end do
- end if
wdun(nVertLevels+1) = 0.
do k=1,nVertLevels
- tend_u(k,iEdge) = tend_u(k,iEdge) - 0.5*(wdun(k+1)+wdun(k))
+ tend_u(k,iEdge) = tend_u(k,iEdge) - 0.5*(wdun(k+1)+wdun(k))
end do
end do
-! --- any mixing for u needs to be added here.
!
- ! ---- vertical mixing - 2nd order
+ ! vertical mixing for u - 2nd order
!
+ if ( v_mom_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdges
+ do iEdge=1,grid % nEdgesSolve
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=2,nVertLevels-1
+
+ z1 = 0.5*(geopotential(k-1,cell1)+geopotential(k-1,cell2))/gravity
+ z2 = 0.5*(geopotential(k ,cell1)+geopotential(k ,cell2))/gravity
+ z3 = 0.5*(geopotential(k+1,cell1)+geopotential(k+1,cell2))/gravity
+ z4 = 0.5*(geopotential(k+2,cell1)+geopotential(k+2,cell2))/gravity
+
+ zm = 0.5*(z1+z2)
+ z0 = 0.5*(z2+z3)
+ zp = 0.5*(z3+z4)
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) + v_mom_eddy_visc2*( &
+ (u(k+1,iEdge)-u(k ,iEdge))/(zp-z0) &
+ -(u(k ,iEdge)-u(k-1,iEdge))/(z0-zm) )/(0.5*(zp-zm))
+ end do
+ end do
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ end if
- if (cell1 > 0 .and. cell2 > 0) then
- do k=2,nVertLevels-1
- z1 = 0.5*(geopotential(k-1,cell1)+geopotential(k-1,cell2))/gravity
- z2 = 0.5*(geopotential(k ,cell1)+geopotential(k ,cell2))/gravity
- z3 = 0.5*(geopotential(k+1,cell1)+geopotential(k+1,cell2))/gravity
- z4 = 0.5*(geopotential(k+2,cell1)+geopotential(k+2,cell2))/gravity
-
- zm = 0.5*(z1+z2)
- z0 = 0.5*(z2+z3)
- zp = 0.5*(z3+z4)
-
- tend_u(k,iEdge) = tend_u(k,iEdge) + vnu*( &
- (u(k+1,iEdge)-u(k ,iEdge))/(zp-z0) &
- -(u(k ,iEdge)-u(k-1,iEdge))/(z0-zm) )/(0.5*(zp-zm))
- end do
- end if
- end do
-
!----------- rhs for theta
tend_theta(:,:) = 0.
@@ -660,89 +693,99 @@
! delsq_theta(:,:) = 0.
!
- ! horizontal mixing - we could combine this with advection directly (i.e. as a turbulent flux),
+ ! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
! but here we can also code in hyperdiffusion if we wish (2nd order at present)
!
+ if ( h_theta_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdges
+ do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
if (cell1 > 0 .and. cell2 > 0) then
+
do k=1,grid % nVertLevels
- theta_turb_flux = hnu*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) + flux
- tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
+ tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) - flux
-! delsq_theta(k,cell1) = delsq_theta(k,iCell) + dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
-! delsq_theta(k,cell2) = delsq_theta(k,iCell) - dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+! delsq_theta(k,cell1) = delsq_theta(k,iCell) + dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+! delsq_theta(k,cell2) = delsq_theta(k,iCell) - dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ end do
- end do
end if
- end do
+ end do
+ end if
+
+
!
- ! horizontal advection
+ ! horizontal advection for theta
!
-
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- theta_edge = 0.5 * (theta(k,cell1) + theta(k,cell2))
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
-! if(u(k,iEdge) .gt. 0.) then
-! theta_edge = theta_edge - (1./6.)*delsq_theta(k,cell1)
-! else
-! theta_edge = theta_edge - (1./6.)*delsq_theta(k,cell2)
-! end if
+ do k=1,grid % nVertLevels
+ theta_edge = 0.5 * (theta(k,cell1) + theta(k,cell2))
- flux = dvEdge (iEdge) * h_edge(k,iEdge) * u(k,iEdge)* theta_edge
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
- end do
- end if
+! if(u(k,iEdge) .gt. 0.) then
+! theta_edge = theta_edge - (1./6.)*delsq_theta(k,cell1)
+! else
+! theta_edge = theta_edge - (1./6.)*delsq_theta(k,cell2)
+! end if
+
+ flux = dvEdge (iEdge) * h_edge(k,iEdge) * u(k,iEdge)* theta_edge
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
+
+ end if
end do
+
!
- ! vertical advection plus diabatic term
- ! Note: we are also dividing through by the cell area after the horizontal flux divergence
+ ! vertical advection plus diabatic term
+ ! Note: we are also dividing through by the cell area after the horizontal flux divergence
!
-
do iCell = 1, nCells
- wdtn(1) = 0.
- do k=2,nVertLevels
- wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
- end do
- wdtn(nVertLevels+1) = 0.
- do k=1,nVertLevels
- tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdnw(k)*(wdtn(k+1)-wdtn(k))
-!! tend_theta(k,iCell) = tend_theta(k) + h(k,iCell)*h_diabatic(k,iCell)
- end do
+ wdtn(1) = 0.
+ do k=2,nVertLevels
+ wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
+ end do
+ wdtn(nVertLevels+1) = 0.
+ do k=1,nVertLevels
+ tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdnw(k)*(wdtn(k+1)-wdtn(k))
+!! tend_theta(k,iCell) = tend_theta(k) + h(k,iCell)*h_diabatic(k,iCell)
+ end do
end do
+
!
- ! ---- vertical mixing - 2nd order
+ ! vertical mixing for theta - 2nd order
!
+ if ( v_theta_eddy_visc2 > 0.0 ) then
- do iCell = 1, nCells
- do k=2,nVertLevels-1
- z1 = geopotential(k-1,iCell)/gravity
- z2 = geopotential(k ,iCell)/gravity
- z3 = geopotential(k+1,iCell)/gravity
- z4 = geopotential(k+2,iCell)/gravity
+ do iCell = 1, grid % nCellsSolve
+ do k=2,nVertLevels-1
+ z1 = geopotential(k-1,iCell)/gravity
+ z2 = geopotential(k ,iCell)/gravity
+ z3 = geopotential(k+1,iCell)/gravity
+ z4 = geopotential(k+2,iCell)/gravity
+
+ zm = 0.5*(z1+z2)
+ z0 = 0.5*(z2+z3)
+ zp = 0.5*(z3+z4)
+
+ tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*h(k,iCell)*( &
+ (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
+ -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
+ end do
+ end do
- zm = 0.5*(z1+z2)
- z0 = 0.5*(z2+z3)
- zp = 0.5*(z3+z4)
+ end if
- tend_theta(k,iCell) = tend_theta(k,iCell) + vnu*prandtl*h(k,iCell)*( &
- (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
- -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
- end do
- end do
-
end subroutine compute_dyn_tend
!---------------------------------------------------------------------------------------------------------
@@ -1119,12 +1162,13 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, r
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt
+ circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, &
+ divergence
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -1138,6 +1182,7 @@
tend_u => s % u % array
circulation => s % circulation % array
vorticity => s % vorticity % array
+ divergence => s % divergence % array
ke => s % ke % array
pv_edge => s % pv_edge % array
pv_vertex => s % pv_vertex % array
@@ -1203,7 +1248,34 @@
end do
end do
+
!
+ ! Compute the divergence at each cell center
+ !
+ divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0) then
+ do k=1,nVertLevels
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ enddo
+ endif
+ if(cell2 > 0) then
+ do k=1,nVertLevels
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ enddo
+ end if
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
+ enddo
+
+
+ !
! Compute kinetic energy in each cell
!
ke(:,:) = 0.0
</font>
</pre>