<p><b>mpetersen@lanl.gov</b> 2012-02-14 10:25:12 -0700 (Tue, 14 Feb 2012)</p><p>Adding variables and stub tendency routines for ztilde.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ztilde/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/Registry 2012-02-13 23:43:23 UTC (rev 1504)
+++ branches/ocean_projects/ztilde/src/core_ocean/Registry 2012-02-14 17:25:12 UTC (rev 1505)
@@ -22,6 +22,9 @@
namelist character grid config_vert_grid_type isopycnal
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
+namelist real grid config_div_low_freq_time 86400.0
+namelist real grid config_h_restore_time 86400.0
+namelist real grid config_h_high_freq_diff 1.0e-5
namelist integer split_explicit_ts config_n_ts_iter 2
namelist integer split_explicit_ts config_n_bcl_iter_beg 2
namelist integer split_explicit_ts config_n_bcl_iter_mid 2
@@ -191,6 +194,8 @@
var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
var persistent real salinity ( nVertLevels nCells Time ) 2 iro salinity state tracers dynamics
var persistent real tracer1 ( nVertLevels nCells Time ) 2 iro tracer1 state tracers testing
+var persistent real hHighFreq ( nVertLevels nCells Time ) 2 r hHighFreq state - -
+var persistent real divergenceLowFreq ( nVertLevels nCells Time ) 2 r divergenceLowFreq state - -
% Tendency variables: neither read nor written to any files
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
@@ -199,6 +204,8 @@
var persistent real tend_temperature ( nVertLevels nCells Time ) 1 - temperature tend tracers dynamics
var persistent real tend_salinity ( nVertLevels nCells Time ) 1 - salinity tend tracers dynamics
var persistent real tend_tracer1 ( nVertLevels nCells Time ) 1 - tracer1 tend tracers testing
+var persistent real tend_hHighFreq ( nVertLevels nCells Time ) 1 - hHighFreq tend - -
+var persistent real tend_divergenceLowFreq ( nVertLevels nCells Time ) 1 - divergenceLowFreq tend - -
% state variables for Split Explicit timesplitting
var persistent real uBtr ( nEdges Time ) 2 - uBtr state - -
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-13 23:43:23 UTC (rev 1504)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-14 17:25:12 UTC (rev 1505)
@@ -106,7 +106,8 @@
config_vert_grid_type.ne.'zlevel'.and. &
config_vert_grid_type.ne.'zstar1'.and. &
config_vert_grid_type.ne.'zstar'.and. &
- config_vert_grid_type.ne.'zstarWeights') then
+ config_vert_grid_type.ne.'zstarWeights'.and. &
+ config_vert_grid_type.ne.'ztilde') then
print *, ' Incorrect choice of config_vert_grid_type.'
call mpas_dmpar_abort(dminfo)
endif
@@ -232,7 +233,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(1) % state, mesh)
call mpas_timer_stop("diagnostic solve", initDiagSolveTimer)
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh)
+ call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh, block % tend)
call ocn_compute_mesh_scaling(mesh)
@@ -287,12 +288,18 @@
end do
! mrp 110808 add end
+ else ! Very first run, not from restart.
- else
- do i=2,nTimeLevs
- call mpas_copy_state(block % state % time_levs(i) % state, &
- block % state % time_levs(1) % state)
- end do
+ do i=2,nTimeLevs
+ call mpas_copy_state(block % state % time_levs(i) % state, &
+ block % state % time_levs(1) % state)
+ end do
+
+ ! for ztilde coordinates, initialize
+ if (config_vert_grid_type.eq.'ztilde') then
+ block % state % time_levs(1) % state % hHighFreq % array(:,:) = 0.0
+ block % state % time_levs(1) % state % divergenceLowFreq % array(:,:) = 0.0
+ endif
endif
end subroutine mpas_init_block!}}}
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_tendency.F 2012-02-13 23:43:23 UTC (rev 1504)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_tendency.F 2012-02-14 17:25:12 UTC (rev 1505)
@@ -397,6 +397,117 @@
!***********************************************************************
!
+! routine ocn_tend_hHighFreq
+!
+!> \brief Compute high frequency thickness tendency
+!> \author Mark Petersen
+!> \date 14 February 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine compute high frequency thickness tendency for z-tilde vertical grid
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_tend_hHighFreq(tend, s, d, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute high frequency thickness tendency for z-tilde vertical grid
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(inout) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:), pointer :: tend_hHighFreq
+
+ integer :: err
+
+ call mpas_timer_start("ocn_tend_hHighFreq")
+
+ tend_hHighFreq => tend % hHighFreq % array
+
+ !
+ ! height tendency: start accumulating tendency terms
+ !
+ tend_hHighFreq = 0.0
+
+ !
+ ! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
+ !
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ ! for explanation of divergence operator.
+ !
+! call mpas_timer_start("hadv", .false., thickHadvTimer)
+! call ocn_thick_hadv_tend(grid, uTransport, h_edge, tend_h, err)
+! call mpas_timer_stop("hadv", thickHadvTimer)
+
+ !
+ ! height tendency: vertical advection term -d/dz(hw)
+ !
+! call mpas_timer_start("vadv", .false., thickVadvTimer)
+! call ocn_thick_vadv_tend(grid, wtop, tend_h, err)
+! call mpas_timer_stop("vadv", thickVadvTimer)
+
+ call mpas_timer_stop("ocn_tend_hHighFreq")
+
+ end subroutine ocn_tend_hHighFreq!}}}
+
+!***********************************************************************
+!
+! routine ocn_tend_divergenceLowFreq
+!
+!> \brief Compute low frequency divergence tendency
+!> \author Mark Petersen
+!> \date 14 February 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine compute low frequency divergence tendency for z-tilde vertical grid
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_tend_divergenceLowFreq(tend, s, d, grid)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute height and normal wind tendencies, as well as diagnostic variables
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(inout) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:), pointer :: tend_divergenceLowFreq
+
+ integer :: err
+
+ call mpas_timer_start("ocn_tend_divergenceLowFreq")
+
+ tend_divergenceLowFreq => tend % divergenceLowFreq % array
+
+ !
+ ! height tendency: start accumulating tendency terms
+ !
+ tend_divergenceLowFreq = 0.0
+
+ call mpas_timer_stop("ocn_tend_divergenceLowFreq")
+
+ end subroutine ocn_tend_divergenceLowFreq!}}}
+
+!***********************************************************************
+!
! routine ocn_diagnostic_solve
!
!> \brief Computes diagnostic variables
@@ -888,7 +999,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_wtop(s1,s2, grid)!{{{
+ subroutine ocn_wtop(s1,s2, grid, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
!
@@ -902,6 +1013,7 @@
type (state_type), intent(inout) :: s1
type (state_type), intent(inout) :: s2
type (mesh_type), intent(in) :: grid
+ type (tend_type), intent(in) :: tend
! mrp 110512 could clean this out, remove pointers?
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
@@ -910,9 +1022,8 @@
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge, tend_hHighFreq
real (kind=RKIND), dimension(:,:), allocatable:: div_hu
real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col, h_weights
@@ -928,6 +1039,8 @@
u => s2 % u % array
wTop => s2 % wTop % array
+ tend_hHighFreq => tend % hHighFreq % array
+
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
maxLevelCell => grid % maxLevelCell % array
@@ -1011,12 +1124,14 @@
! mode, div_hu_btr, among all the layers. Distribute in proportion
! to the layer thickness.
+ h_weights = 1.0
+
do iCell=1,nCells
hSum = 0.0
do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + h(k,iCell)
+ h_tend_col(k) = - h_weights(k)*h(k,iCell)*div_hu_btr(iCell)
+ hSum = hSum + h_weights(k)*h(k,iCell)
end do
h_tend_col = h_tend_col / hSum
@@ -1057,6 +1172,32 @@
end do
end do
+ elseif (config_vert_grid_type.eq.'ztilde') then
+
+ ! This is identical to zstar, but also includes a high frequency
+ ! thickness tendency.
+
+ h_weights = 1.0
+
+ do iCell=1,nCells
+
+ hSum = 0.0
+ do k=1,maxLevelCell(iCell)
+ h_tend_col(k) = - h_weights(k)*h(k,iCell)*div_hu_btr(iCell) &
+ + tend_hHighFreq(k,iCell)
+ hSum = hSum + h_weights(k)*h(k,iCell)
+ end do
+ h_tend_col = h_tend_col / hSum
+
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
+ end do
+ end do
+
endif
deallocate(div_hu, div_hu_btr, h_tend_col, h_weights)
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-13 23:43:23 UTC (rev 1504)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-14 17:25:12 UTC (rev 1505)
@@ -111,6 +111,10 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ if (config_vert_grid_type.eq.'ztilde') then
+ block % state % time_levs(2) % state % hHighFreq% array(:,:) = block % state % time_levs(1) % state % hHighFreq% array(:,:)
+ block % state % time_levs(2) % state % divergenceLowFreq % array(:,:) = block % state % time_levs(1) % state % divergenceLowFreq % array(:,:)
+ endif
do iCell=1,block % mesh % nCells ! couple tracers to h
do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
@@ -163,12 +167,18 @@
! --- compute tendencies
+
call mpas_timer_start("RK4-tendency computations")
block => domain % blocklist
do while (associated(block))
+ if (config_vert_grid_type.eq.'ztilde') then
+ !call ocn_tend_hHighFreq(block % tend, provis, block % diagnostics, block % mesh)
+ !call ocn_tend_divergenceLowFreq(block % tend, provis, block % diagnostics, block % mesh)
+ endif
+
! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(provis, provis, block % mesh)
+ call ocn_wtop(provis, provis, block % mesh, block % tend)
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, provis, block % diagnostics, err)
@@ -201,6 +211,17 @@
call mpas_dmpar_exch_halo_field3d_real(domain % dminfo, block % tend % tracers % array(:,:,:), &
block % tend % num_tracers, block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+ if (config_vert_grid_type.eq.'ztilde') then
+
+ call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % hHighFreq % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % divergenceLowFreq % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ endif
+
block => block % next
end do
call mpas_timer_stop("RK4-pronostic halo update")
@@ -227,6 +248,17 @@
end do
end do
+
+ if (config_vert_grid_type.eq.'ztilde') then
+
+ provis % hHighFreq % array(:,:) = block % state % time_levs(1) % state % hHighFreq % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % hHighFreq % array(:,:)
+
+ provis % divergenceLowFreq % array(:,:) = block % state % time_levs(1) % state % divergenceLowFreq % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % divergenceLowFreq % array(:,:)
+
+ end if
+
if (config_test_case == 1) then ! For case 1, wind field should be fixed
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
@@ -259,6 +291,16 @@
end do
end do
+ if (config_vert_grid_type.eq.'ztilde') then
+ block % state % time_levs(2) % state % hHighFreq % array(:,:) &
+ = block % state % time_levs(2) % state % hHighFreq % array(:,:) &
+ + rk_weights(rk_step) * block % tend % hHighFreq % array(:,:)
+
+ block % state % time_levs(2) % state % divergenceLowFreq % array(:,:) &
+ = block % state % time_levs(2) % state % divergenceLowFreq % array(:,:) &
+ + rk_weights(rk_step) * block % tend % divergenceLowFreq % array(:,:)
+ end if
+
block => block % next
end do
call mpas_timer_stop("RK4-RK4 accumulate update")
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-13 23:43:23 UTC (rev 1504)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-14 17:25:12 UTC (rev 1505)
@@ -718,7 +718,7 @@
block => domain % blocklist
do while (associated(block))
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
+ call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh,block % tend )
call ocn_tend_h (block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
</font>
</pre>