<p><b>qchen3@fsu.edu</b> 2013-04-12 10:07:24 -0600 (Fri, 12 Apr 2013)</p><p><br>
BRANCH COMMIT<br>
<br>
Code compiles and runs. uBolus and the Redi diffusion are computed but NOT applied <br>
in the simulation. Their correctness have not been verified.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/gm_implement/namelist.input.ocean
===================================================================
--- branches/ocean_projects/gm_implement/namelist.input.ocean 2013-04-11 21:37:04 UTC (rev 2741)
+++ branches/ocean_projects/gm_implement/namelist.input.ocean 2013-04-12 16:07:24 UTC (rev 2742)
@@ -20,7 +20,7 @@
/
&time_integration
config_dt = 500.0
- config_time_integrator = 'split_explicit'
+ config_time_integrator = 'RK4'
/
&grid
config_num_halos = 3
@@ -61,7 +61,10 @@
config_Leith_visc2_max = 2.5e3
/
&standard_GM
- config_h_kappa = 10.0
+ config_use_GM = .true.
+ config_use_Redi_diff = .true.
+ config_h_kappa = 600.0
+ config_Redi_kappa = 0.0
config_h_kappa_q = 0.0
/
&Rayleigh_damping
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F 2013-04-11 21:37:04 UTC (rev 2741)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F 2013-04-12 16:07:24 UTC (rev 2742)
@@ -122,6 +122,7 @@
end do
end do
+
! Compute vertical derivative of density (rhoz) at cell center and layer interface
do iCell = 1, nCells
do k = 2, maxLevelCell(iCell)
@@ -135,27 +136,34 @@
rhoz(maxLevelCell(iCell)+1,iCell) = 0.0
end do
+!WRITE(*,*) 'DEBUG MARK #1'
+
! Interpolate grad_rho and grad_zMid to layer interface
do iEdge = 1, nEdges
- do k = 2, maxLevelEdgeTop(iEdge)
- h1 = h_edge(k-1,iEdge)
- h2 = h_edge(k,iEdge)
+ ! The interpolation can only be carried out on non-boundary edges
+ if (maxLevelEdgeTop(iEdge) .GE. 1) then
+ do k = 2, maxLevelEdgeTop(iEdge)
+ h1 = h_edge(k-1,iEdge)
+ h2 = h_edge(k,iEdge)
- ! Using second-order interpolation below
- grad_rho_interface(k,iEdge) = (h2 * grad_rho(k-1,iEdge) + h1 * grad_rho(k,iEdge)) / (h1 + h2)
- grad_zMid_interface(k,iEdge) = (h2 * grad_zMid(k-1,iEdge) + h1 * grad_zMid(k,iEdge)) / (h1 + h2)
+ ! Using second-order interpolation below
+ grad_rho_interface(k,iEdge) = (h2 * grad_rho(k-1,iEdge) + h1 * grad_rho(k,iEdge)) / (h1 + h2)
+ grad_zMid_interface(k,iEdge) = (h2 * grad_zMid(k-1,iEdge) + h1 * grad_zMid(k,iEdge)) / (h1 + h2)
- end do
+ end do
- ! Approximation of values on the top and bottom interfaces through the idea of having ghost cells above
- ! the top and below the bottom layers of the same depths and density.
- grad_rho_interface(1,iEdge) = grad_rho(1,iEdge)
- grad_rho_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_rho(maxLevelEdgeTop(iEdge),iEdge)
- grad_zMid_interface(1,iEdge) = grad_zMid(1,iEdge)
- grad_zMid_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_zMid(maxLevelEdgeTop(iEdge),iEdge)
+ ! Approximation of values on the top and bottom interfaces through the idea of having ghost cells above
+ ! the top and below the bottom layers of the same depths and density.
+ grad_rho_interface(1,iEdge) = grad_rho(1,iEdge)
+ grad_rho_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_rho(maxLevelEdgeTop(iEdge),iEdge)
+ grad_zMid_interface(1,iEdge) = grad_zMid(1,iEdge)
+ grad_zMid_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_zMid(maxLevelEdgeTop(iEdge),iEdge)
+ end if
end do
+!WRITE(*,*) 'DEBUG MARK #2'
+
! Interpolate rhoz to cell edge and layer interface
do iEdge = 1, nEdges
do k = 1, maxLevelEdgeTop(iEdge)+1
@@ -175,6 +183,8 @@
! Compute the density gradient along constant z surfaces.
grad_rho_constZ(:,:) = grad_rho_interface(:,:) - rhoz_edge(:,:) * grad_zMid_interface(:,:)
+!WRITE(*,*) 'DEBUG MARK #3'
+
! Compute slopeRelative at cell edge and layer interface
! QC Note: We may need to overwrite the invalid or unrealistically large values in the variable.
slopeRelative(:,:) = -(1.0 - config_diapycnal_diff) * grad_rho_interface(:,:) / rhoz_edge(:,:)
@@ -199,6 +209,8 @@
end do
end do
+!WRITE(*,*) 'DEBUG MARK #4'
+
!
!Compute the stream function
!
@@ -209,36 +221,39 @@
cell2 = cellsOnEdge(2,iEdge)
! Construct the tridiagonal matrix
- ! First row
- k = 2
- BruntVaisalaFreqTopEdge = 0.5 * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
- tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
- tridiagC(k-1) = 2.*c/h_edge(k,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
- rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
-
- ! Second to next to the last rows
- do k = 3, maxLevelEdgeTop(iEdge)-1
+ if (maxLevelEdgeTop(iEdge) .GE. 3) then
+ ! First row
+ k = 2
BruntVaisalaFreqTopEdge = 0.5 * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
- tridiagA(k-2) = 2.*c/h_edge(k-1,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
tridiagC(k-1) = 2.*c/h_edge(k,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
- end do
- ! Last row
- k = maxLevelEdgeTop(iEdge)
- tridiagA(k-2) = 2.*c/h_edge(k-1,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
- tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
- rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
+ ! Second to next to the last rows
+ do k = 3, maxLevelEdgeTop(iEdge)-1
+ BruntVaisalaFreqTopEdge = 0.5 * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
+ tridiagA(k-2) = 2.*c/h_edge(k-1,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
+ tridiagC(k-1) = 2.*c/h_edge(k,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
+ end do
- ! Total number of rows
- N = maxLevelEdgeTop(iEdge) - 1
+ ! Last row
+ k = maxLevelEdgeTop(iEdge)
+ tridiagA(k-2) = 2.*c/h_edge(k-1,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
+ rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
+
+ ! Total number of rows
+ N = maxLevelEdgeTop(iEdge) - 1
- ! Call the tridiagonal solver
- call tridiagonal_solve(tridiagA, tridiagB, tridiagC, rightHandSide, gmStreamFunc(2:maxLevelEdgeTop(iEdge),iEdge), N)
-
+ ! Call the tridiagonal solver
+ call tridiagonal_solve(tridiagA, tridiagB, tridiagC, rightHandSide, gmStreamFunc(2:maxLevelEdgeTop(iEdge),iEdge), N)
+ end if
end do
+!WRITE(*,*) 'DEBUG MARK #5'
+
! Compute uBolusGM from the stream function
do iEdge = 1, nEdges
do k = 1, maxLevelEdgeTop(iEdge)
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-04-11 21:37:04 UTC (rev 2741)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-04-12 16:07:24 UTC (rev 2742)
@@ -207,16 +207,17 @@
call ocn_diagnostic_solve(dt, block % provis, block % mesh)
- ! Compute uBolusGM and then tracer transport velocity
- if (config_h_kappa .GE. epsilon(0D0)) then
+ block % provis % uTransport % array(:,:) = block % provis % u % array(:,:)
+
+ ! Compute uBolusGM, slopeRelative and RediDiffVertCoef if respective flags are turned on
+ if (config_use_GM .OR. config_use_Redi_diff) then
call ocn_gm_compute_uBolus(block % provis, block % mesh, block % scratch)
- else
- block % provis % uBolusGM % array(:,:) = 0.0
end if
- if (config_use_GM) then
+ !if (config_use_GM) then
+ if ( .False. ) then ! For testing and debuging: compute uBolusGM but leave them inactive
block % provis % uTransport % array(:,:) &
- = block % provis % u % array(:,:) &
+ = block % provis % uTransport % array(:,:) &
+ block % provis % uBolusGM % array(:,:)
end if
@@ -283,16 +284,17 @@
! implicit vmix routine that follows.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ block % state % time_levs(2) % state % uTransport % array(:,:) = block % state % time_levs(2) % state % u % array(:,:)
+
! Compute uBolusGM and the tracer transport velocity
- if (config_h_kappa .GE. epsilon(0D0)) then
+ if (config_use_GM .OR. config_use_Redi_diff) then
call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh, block % scratch)
- else
- block % state % time_levs(2) % state % uBolusGM % array(:,:) = 0.0
end if
- if (config_use_GM) then
+ !if (config_use_GM) then
+ if (.false.) then ! For testing only: compute uBolusGM et al but leave them inactive
block % state % time_levs(2) % state % uTransport % array(:,:) &
- = block % state % time_levs(2) % state % u % array(:,:) &
+ = block % state % time_levs(2) % state % uTransport % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
end if
@@ -323,16 +325,17 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ block % state % time_levs(2) % state % uTransport % array(:,:) = block % state % time_levs(2) % state % u % array(:,:)
+
! Compute uBolusGM and the tracer transport velocity
- if (config_h_kappa .GE. epsilon(0D0)) then
+ if (config_use_GM .OR. config_use_Redi_diff) then
call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh, block % scratch)
- else
- block % state % time_levs(2) % state % uBolusGM % array(:,:) = 0.0
end if
- if (config_use_GM) then
+ !if (config_use_GM) then
+ if (.false.) then ! For testing: compute uBolus but leave it inactive
block % state % time_levs(2) % state % uTransport % array(:,:) &
- = block % state % time_levs(2) % state % u % array(:,:) &
+ = block % state % time_levs(2) % state % uTransport % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
end if
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-04-11 21:37:04 UTC (rev 2741)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-04-12 16:07:24 UTC (rev 2742)
@@ -217,13 +217,17 @@
! COMPUTE the extra terms arising due to mismatch between the constant coordiate surfaces and the
! isopycnal surfaces.
!
- ! Compute vertical derivative of tracers at cell center and layer interface
- do iTracer = 1, num_tracers
- do iCell = 1, nCells
- do k = 2, maxLevelCell(iCell)
- dhTracerdZ(k,iCell) = (h(k-1,iCell)*tracers(iTracer,k-1,iCell) - h(k,iCell)*tracers(iTracer,k,iCell)) / (zMid(k-1,iCell) - zMid(k,iCell))
- end do
+ !if (config_use_Redi_diff) then
+ if ( .false. ) then ! For testing: leave redi diff inactive
+
+ ! Compute vertical derivative of tracers at cell center and layer interface
+ do iTracer = 1, num_tracers
+ do iCell = 1, nCells
+ do k = 2, maxLevelCell(iCell)
+ dhTracerdZ(k,iCell) = (h(k-1,iCell)*tracers(iTracer,k-1,iCell) - h(k,iCell)*tracers(iTracer,k,iCell)) / (zMid(k-1,iCell) - zMid(k,iCell))
+ end do
+
! Approximation of dhTracerdZ on the top and bottom interfaces through the idea of having
! ghost cells above the top and below the bottom layers of the same depths and tracer density.
! Essentially, this enforces the boundary condition (d tracer)/dz = 0 at the top and bottom.
@@ -302,6 +306,7 @@
end do
end do
+ end if ! config_use_Redi_diff
call mpas_deallocate_scratch_field(scratch % grad_hTracer)
call mpas_deallocate_scratch_field(scratch % grad_hTracer_interface)
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_vmix_coefs_redi.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_vmix_coefs_redi.F 2013-04-11 21:37:04 UTC (rev 2741)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_vmix_coefs_redi.F 2013-04-12 16:07:24 UTC (rev 2742)
@@ -121,7 +121,10 @@
vertDiffTopOfCell => d % vertDiffTopOfCell % array
vertRediDiff => s % k33 % array
- call ocn_tracer_vmix_coefs_redi(grid, vertDiffTopOfCell, vertRediDiff, err)
+ !if (config_use_Redi_diff) then
+ if (.False.) then
+ call ocn_tracer_vmix_coefs_redi(grid, vertDiffTopOfCell, vertRediDiff, err)
+ end if
!--------------------------------------------------------------------
</font>
</pre>