<p><b>mpetersen@lanl.gov</b> 2011-08-18 16:46:58 -0600 (Thu, 18 Aug 2011)</p><p>Higdon split with a few bug corrections, and lots of added commenting. This version includes -1/gamma u term on bcl and btr mom equations, and a del2 operator on btr mom eqn.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-08-18 22:26:56 UTC (rev 947)
+++ branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-08-18 22:46:58 UTC (rev 948)
@@ -1,26 +1,23 @@
&sw_model
config_test_case = 0
config_time_integration = 'higdon_unsplit'
+ config_rk_filter_btr_mode = .false.
config_dt = 10.0
config_ntimesteps =100
config_output_interval = 100
config_stats_interval = 100
/
-
&io
/
-
&restart
config_restart_interval = 10000000
config_do_restart = .false.
config_restart_time = 1036800.0
/
-
&grid
config_vert_grid_type = 'zlevel'
config_rho0 = 1000
/
-
&timestep_higdon
config_n_ts_iter = 2
config_n_bcl_iter_beg = 1
@@ -30,12 +27,18 @@
config_n_btr_cor_iter = 2
config_compute_tr_midstage = .true.
config_u_correction = .true.
+ config_filter_btr_mode false
+ config_btr_mom_decay = .false.
+ config_btr_mom_decay_time = 3600.0
+ config_btr_mom_eddy_visc2 = 0.0
/
&hmix
config_h_mom_eddy_visc2 = 1.0e5
config_h_mom_eddy_visc4 = 0.0
config_h_tracer_eddy_diff2 = 1.0e4
config_h_tracer_eddy_diff4 = 0.0
+ config_mom_decay = .false.
+ config_mom_decay_time = 3600.0
/
&vmix
config_vert_visc_type = 'const'
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-08-18 22:26:56 UTC (rev 947)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-08-18 22:46:58 UTC (rev 948)
@@ -3,6 +3,7 @@
#
namelist integer sw_model config_test_case 5
namelist character sw_model config_time_integration RK4
+namelist logical sw_model config_rk_filter_btr_mode false
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
@@ -25,12 +26,17 @@
namelist logical timestep_higdon config_compute_tr_midstage true
namelist logical timestep_higdon config_u_correction true
namelist logical timestep_higdon config_filter_btr_mode false
+namelist logical timestep_higdon config_btr_mom_decay false
+namelist real timestep_higdon config_btr_mom_decay_time 3600.0
+namelist real timestep_higdon config_btr_mom_eddy_visc2 0.0
namelist logical sw_model config_h_ScaleWithMesh false
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
namelist real hmix config_h_tracer_eddy_diff4 0.0
namelist real hmix config_apvm_upwinding 0.5
+namelist logical hmix config_mom_decay false
+namelist real hmix config_mom_decay_time 3600.0
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
namelist logical vmix config_implicit_vertical_mix .true.
@@ -193,6 +199,10 @@
var persistent real FBtr ( nEdges Time ) 1 o FBtr state - -
var persistent real GBtrForcing ( nEdges Time ) 1 o GBtrForcing state - -
var persistent real uBcl ( nVertLevels nEdges Time ) 2 o uBcl state - -
+var persistent real circulationBtr ( nVertices Time ) 1 - circulationBtr state - -
+var persistent real divergenceBtr ( nCells Time ) 1 - divergenceBtr state - -
+var persistent real vorticityBtr ( nVertices Time ) 1 - vorticityBtr state - -
+var persistent real u_diffusionBtr ( nEdges Time ) 1 - u_diffusionBtr state - -
# Diagnostic fields: only written to output
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
@@ -236,4 +246,3 @@
var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 - vertViscTopOfEdge diagnostics - -
var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells Time ) 1 - vertDiffTopOfCell diagnostics - -
-
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-08-18 22:26:56 UTC (rev 947)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-08-18 22:46:58 UTC (rev 948)
@@ -80,10 +80,10 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+
call compute_mesh_scaling(mesh)
call rbfInterp_initialize(mesh)
-! mrp temp
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
@@ -119,13 +119,22 @@
if (.not. config_do_restart) then
block % state % time_levs(1) % state % xtime % scalar = 0.0
+
+! mrp 110808 add, so that variables are copied to * variables for Higdon
+ do i=2,nTimeLevs
+ call copy_state(block % state % time_levs(i) % state, &
+ block % state % time_levs(1) % state)
+ end do
+! mrp 110808 add end
+
+
else
do i=2,nTimeLevs
call copy_state(block % state % time_levs(i) % state, &
block % state % time_levs(1) % state)
end do
endif
-
+
end subroutine mpas_init_block
@@ -148,8 +157,8 @@
! Eventually, dt should be domain specific
dt = config_dt
ntimesteps = config_ntimesteps
-
- call write_output_frame(output_obj, output_frame, domain)
+ ! mrp 110808 turn off initial output frame
+! call write_output_frame(output_obj, output_frame, domain)
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
@@ -364,6 +373,7 @@
block % state % time_levs(1) % state % ssh % array(cell1) &
+ block % state % time_levs(1) % state % ssh % array(cell2) )
+ ! uBtr = sum(u)/sum(h) on each column
uhSum = (sshEdge + block % mesh % hZLevel % array(1)) &
* block % state % time_levs(1) % state % u % array(1,iEdge)
hSum = sshEdge + block % mesh % hZLevel % array(1)
@@ -377,29 +387,30 @@
enddo
block % state % time_levs(1) % state % uBtr % array(iEdge) = uhSum/hsum
+ ! uBcl(k,iEdge) = u(k,iEdge) - uBtr(iEdge)
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
= block % state % time_levs(1) % state % u % array(k,iEdge) &
- block % state % time_levs(1) % state % uBtr % array(iEdge)
enddo
+ ! uBcl=0, u=0 on land cells
do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1, block % mesh % nVertLevels
block % state % time_levs(1) % state % uBcl % array(k,iEdge) = 0.0
block % state % time_levs(1) % state % u % array(k,iEdge) = 0.0
enddo
enddo
- if (config_filter_btr_mode) then
- ! filter uBtr out of initial condition
- block % state % time_levs(1) % state % u % array(:,:) &
- = block % state % time_levs(1) % state % uBcl % array(:,:)
+ if (config_filter_btr_mode) then
+ ! filter uBtr out of initial condition
+ block % state % time_levs(1) % state % u % array(:,:) &
+ = block % state % time_levs(1) % state % uBcl % array(:,:)
- block % state % time_levs(1) % state % uBtr % array(:) = 0.0
- endif
+ block % state % time_levs(1) % state % uBtr % array(:) = 0.0
+ endif
endif
-
!print *, '11 u ',minval(domain % blocklist % state % time_levs(1) % state % u % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
! maxval(domain % blocklist % state % time_levs(1) % state % u % array(:,1:domain % blocklist % mesh % nEdgesSolve))
!print *, '11 uBtr ',minval(domain % blocklist % state % time_levs(1) % state % uBtr % array(1:domain % blocklist % mesh % nEdgesSolve)), &
@@ -409,24 +420,23 @@
! mrp temp testing - is uBcl vert sum zero?
- do iEdge=1,block % mesh % nEdges
- uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * block % state % time_levs(1) % state % uBcl % array(1,iEdge)
- hSum = sshEdge + block % mesh % hZLevel % array(1)
+! do iEdge=1,block % mesh % nEdges
+! uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * block % state % time_levs(1) % state % uBcl % array(1,iEdge)
+! hSum = sshEdge + block % mesh % hZLevel % array(1)
- do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + block % mesh % hZLevel % array(k) * block % state % time_levs(1) % state % uBcl % array(k,iEdge)
- hSum = hSum + block % mesh % hZLevel % array(k)
- enddo
- block % state % time_levs(1) % state % FBtr % array(iEdge) = uhSum/hSum
+! do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+! uhSum = uhSum + block % mesh % hZLevel % array(k) * block % state % time_levs(1) % state % uBcl % array(k,iEdge)
+! hSum = hSum + block % mesh % hZLevel % array(k)
+! enddo
+! block % state % time_levs(1) % state % FBtr % array(iEdge) = uhSum/hSum
- enddo ! iEdge
+! enddo ! iEdge
!print *, 'uBcl vert sum IC',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve)), &
! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve))
! mrp temp testing - is uBcl vert sum zero? end
-
block => block % next
end do
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-08-18 22:26:56 UTC (rev 947)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-08-18 22:46:58 UTC (rev 948)
@@ -81,6 +81,7 @@
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
integer :: nCells, nEdges, nVertLevels, num_tracers
+ real (kind=RKIND) :: coef
real (kind=RKIND), dimension(:,:), pointer :: &
u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell, ke_edge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -162,6 +163,13 @@
end if
call compute_tend_h(block % tend, provis, block % diagnostics, block % mesh)
call compute_tend_u(block % tend, provis, block % diagnostics, block % mesh)
+
+ ! mrp 110718 filter btr mode out of u_tend
+ ! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
+ if (config_rk_filter_btr_mode) then
+ call filter_btr_mode_tend_u(block % tend, provis, block % diagnostics, block % mesh)
+ endif
+
call compute_scalar_tend(block % tend, provis, block % diagnostics, block % mesh)
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
@@ -191,6 +199,13 @@
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % u % array(:,:)
+
+ ! mrp 110718 filter btr mode out of u
+ if (config_rk_filter_btr_mode) then
+ ! call filter_btr_mode_u(provis, block % mesh)
+ !block % tend % h % array(:,:) = 0.0 ! I should not need this
+ endif
+
provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
@@ -206,6 +221,7 @@
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
+
call compute_solve_diagnostics(dt, provis, block % mesh)
block => block % next
@@ -220,6 +236,12 @@
do while (associated(block))
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
+ rk_weights(rk_step) * block % tend % u % array(:,:)
+ ! mrp 110718 filter btr mode out of u
+ if (config_rk_filter_btr_mode) then
+ ! call filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
+ !block % tend % h % array(:,:) = 0.0 ! I should not need this
+ endif
+
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
+ rk_weights(rk_step) * block % tend % h % array(:,:)
@@ -310,6 +332,12 @@
end if
end do
+ ! mrp 110718 filter btr mode out of u
+ if (config_rk_filter_btr_mode) then
+ call filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
+ !block % tend % h % array(:,:) = 0.0 ! I should not need this
+ endif
+
!
! Implicit vertical solve for tracers
!
@@ -340,13 +368,33 @@
deallocate(A,C,uTemp,tracersTemp)
end if
+ ! mrp 110725 momentum decay term
+ if (config_mom_decay) then
+
+ !
+ ! Implicit solve for momentum decay
+ !
+ ! Add term to RHS of momentum equation: -1/gamma u
+ !
+ ! This changes the solve to:
+ ! u^{n+1} = u_provis^{n+1}/(1+dt/gamma)
+ !
+ coef = 1.0/(1.0 + dt/config_mom_decay_time)
+ do iEdge=1,block % mesh % nEdges
+ do k=1,maxLevelEdgeTop(iEdge)
+ u(k,iEdge) = coef*u(k,iEdge)
+ end do
+ end do
+
+ end if
+
+
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
-! mrp temp
call reconstruct(block % state % time_levs(2) % state, block % mesh)
block => block % next
@@ -374,11 +422,12 @@
integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
- n_bcl_iter(config_n_ts_iter)
+ n_bcl_iter(config_n_ts_iter), &
+ vertex1, vertex2, iVertex
type (block_type), pointer :: block
- real (kind=RKIND) :: uhSum, hSum, grav_rho0Inv, sshEdge, flux, &
- uPerp, uCorr, tracerTemp
+ real (kind=RKIND) :: uhSum, hSum, sshEdge, flux, &
+ uPerp, uCorr, tracerTemp, coef
integer :: num_tracers, ucorr_coef
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -390,21 +439,19 @@
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
- grav_rho0Inv = gravity/config_rho0
-
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Prep variables before first iteration
!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
block => domain % blocklist
do while (associated(block))
-!mrp 110512 need to add w, p, and ssh here?
-
do iEdge=1,block % mesh % nEdges
! The baroclinic velocity needs be recomputed at the beginning of a
! timestep because the implicit vertical mixing is conducted on the
- ! total u. We keep uBcl from the previous timestep.
+ ! total u. We keep uBtr from the previous timestep.
block % state % time_levs(1) % state % uBcl % array(:,iEdge) &
= block % state % time_levs(1) % state % u % array(:,iEdge) &
- block % state % time_levs(1) % state % uBtr % array(iEdge)
@@ -417,12 +464,10 @@
enddo ! iEdge
+ ! Initialize * variables that are used compute baroclinic tendencies below.
block % state % time_levs(2) % state % ssh % array(:) &
= block % state % time_levs(1) % state % ssh % array(:)
-! Are the following two needed?
-! block % state % time_levs(2) % state % ssh_edge % array(:,:) &
-! = block % state % time_levs(1) % state % ssh_edge % array(:,:)
block % state % time_levs(2) % state % h_edge % array(:,:) &
= block % state % time_levs(1) % state % h_edge % array(:,:)
@@ -433,6 +478,7 @@
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
= block % state % time_levs(1) % state % tracers % array(:,k,iCell)
end do
+
end do
block => block % next
@@ -468,9 +514,11 @@
block => block % next
end do
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stage 1: Baroclinic velocity (3D) prediction, explicit with long timestep
!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! compute velocity tendencies, T(u*,w*,p*)
@@ -480,7 +528,6 @@
call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
end if
call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
-! mrp soon: for split version: add g/rho grad ssh to tend_u
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
end do
@@ -510,21 +557,17 @@
uTemp = 0.0 ! could put this after with uTemp(maxleveledgetop+1:nvertlevels)=0
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- !{\bf G}_k = -f{\bf u}_{k,n+1}^{'\;\perp} + {\bf T}({\bf u}_k^*,w_k^*,p_k^*) +\frac{g}{\rho_0}</font>
<font color="gray">abla \ssh^*
-! This soon, with coriolis:
-! G(k) = tend_u(k,iEdge) -fEdge(iEdge)*uBclPerp(k,iEdge) - g/rho0 grad (h_1^*)
+ ! uBclNew = uBclOld + dt*(-f*uBclPerp + T(u*,w*,p*) + g*grad(SSH*) )
! Here uNew is a work variable containing -fEdge(iEdge)*uBclPerp(k,iEdge)
-! mrp temp new one:
uTemp(k) &
= block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ dt * (block % tend % u % array (k,iEdge) &
- + block % state % time_levs(2) % state % u % array (k,iEdge) &
- + split*grav_rho0Inv &
+ + block % state % time_levs(2) % state % u % array (k,iEdge) & ! this is f*uBcl^{perp}
+ + split*gravity &
*( block % state % time_levs(2) % state % ssh % array(cell2) &
- block % state % time_levs(2) % state % ssh % array(cell1) ) &
/block % mesh % dcEdge % array(iEdge) )
-
enddo
! Compute GBtrForcing, the vertically averaged forcing
@@ -553,27 +596,8 @@
+ uTemp(k) - dt * block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
enddo
-! mrp temp testing - is uBcl vert sum zero?
-! uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * block % state % time_levs(1) % state % uBcl % array(1,iEdge)
-! hSum = sshEdge + block % mesh % hZLevel % array(1)
-! do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
-! uhSum = uhSum + block % mesh % hZLevel % array(k) * block % state % time_levs(1) % state % uBcl % array(k,iEdge)
-! hSum = hSum + block % mesh % hZLevel % array(k)
-! enddo
-! block % state % time_levs(1) % state % FBtr % array(iEdge) = uhSum/hSum
-! mrp temp testing end
-
enddo ! iEdge
-! mrp temp testing - is uBcl vert sum zero?
-!print *, '12 uBcl ',minval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
-! maxval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve))
-
-!print *, 'uBcl vert sum ',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve)), &
-! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve))
-! mrp temp testing - is uBcl vert sum zero? end
-
-
deallocate(uTemp)
block => block % next
@@ -594,9 +618,11 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stage 2: Barotropic velocity (2D) prediction, explicitly subcycled
!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
oldBtrSubcycleTime = 1
newBtrSubcycleTime = 2
@@ -606,6 +632,7 @@
block => domain % blocklist
do while (associated(block))
+ ! For Higdon unsplit, simply set uBtrNew=0, uBtrSubcycle=0, and uNew=uBclNew
block % state % time_levs(2) % state % uBtr % array(:) = 0.0
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:) = 0.0
@@ -658,34 +685,106 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do j=1,config_n_btr_subcycles
- ! Compute thickness flux and new SSH
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: initial solve for velecity
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ uPerpTime = oldBtrSubcycleTime
+
block => domain % blocklist
do while (associated(block))
-!! mrp del allocate(tend_ssh(block % mesh % nCells))
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ ! Compute -f*uPerp
+ uPerp = 0.0
+ do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
+ eoe = block % mesh % edgesOnEdge % array(i,iEdge)
+ uPerp = uPerp + block % mesh % weightsOnEdge % array(i,iEdge) &
+ * block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ * block % mesh % fEdge % array(eoe)
+ end do
+
+ ! mrp 110606 efficiency note: could make this a 1D integer factor instead of an if statement.
+ if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) = 0.0
+ else
+
+ ! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + dt/config_n_btr_subcycles *( &
+ uPerp &
+ - gravity &
+ *( block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ - block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
+ /block % mesh % dcEdge % array(iEdge) &
+ + block % state % time_levs(1) % state % GBtrForcing % array(iEdge) )
+
+ endif
+
+ end do
+
+ ! Implicit solve for barotropic momentum decay
+ if ( config_btr_mom_decay) then
+ !
+ ! Add term to RHS of momentum equation: -1/gamma u
+ !
+ ! This changes the solve to:
+ ! u^{n+1} = u_provis^{n+1}/(1+dt/gamma)
+ !
+ coef = 1.0/(1.0 + dt/config_n_btr_subcycles/config_btr_mom_decay_time)
+ do iEdge=1,block % mesh % nEdges
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ * coef
+ end do
+
+ endif
+
+
+ block => block % next
+ end do ! block
+
+
+ ! boundary update on uBtrNew
+ block => domain % blocklist
+ do while (associated(block))
+
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:), &
+ block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ block => block % next
+ end do ! block
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Compute thickness flux and new SSH
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ block => domain % blocklist
+ do while (associated(block))
+
block % tend % ssh % array(:) = 0.0
+ ! Flux = 1/2*(uBtrNew+uBtrOld)*H where H it total thickness
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
-! FBtr(iEdge) = FBtr(iEdge) + flux
-! tend_ssh(cell1) = tend_ssh(cell1) - flux
-! tend_ssh(cell2) = tend_ssh(cell2) + flux
-
sshEdge = 0.5 &
*( block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
-! \zeta_{n+j/J}^{edge} = Interp(\zeta_{n+j/J})
-! {\bf F}_j = \ubtr_{n+j/J}
-! \left(\zeta_{n+j/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)
- flux = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ flux = 0.5*( block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ +block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* block % mesh % dvEdge % array(iEdge) &
* (sshEdge + hSum)
+
block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux
block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux
@@ -694,11 +793,9 @@
+ flux
end do
+ ! SSHnew = SSHold + dt/J*(-div(Flux))
do iCell=1,block % mesh % nCells
-! \zeta_{n+(j+1)/J} = \zeta_{n+j/J} + \frac{\Delta t}{J} \left( -
-! </font>
<font color="gray">abla \cdot {\bf F}_j \right)
-! sshSubcycleNew(iCell) = sshSubcycleOld(iCell) + dt/J*tend_ssh(iCell)/areaCell(iCell)
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ dt/config_n_btr_subcycles &
@@ -706,11 +803,10 @@
end do
-! mrp del deallocate(tend_ssh)
block => block % next
end do ! block
-! boundary update on \zeta_{n+(j+1)/J}
+ ! boundary update on SSNnew
block => domain % blocklist
do while (associated(block))
@@ -729,39 +825,106 @@
do iCell=1,block % mesh % nCells
-
-! sshNew(iCell) = sshNew(iCell) + sshSubcycleNew(iCell)
+ ! Accumulate SSH in running sum over the subcycles.
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(2) % state % ssh % array(iCell) &
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
end do
-! mrp del deallocate(tend_ssh)
+ block => block % next
+ end do ! block
+! mrp 110801
+! compute btr_divergence and btr_vorticity for del2(u_btr)
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % u_diffusionBtr % array(:) = 0.0
+ if ( config_btr_mom_eddy_visc2 > 0.0 ) then
+ !
+ ! Compute circulation and relative vorticity at each vertex
+ !
+ block % state % time_levs(1) % state % circulationBtr % array(:) = 0.0
+ do iEdge=1,block % mesh % nEdges
+ vertex1 = block % mesh % verticesOnEdge % array(1,iEdge)
+ vertex2 = block % mesh % verticesOnEdge % array(2,iEdge)
+ block % state % time_levs(1) % state % circulationBtr % array(vertex1) &
+ = block % state % time_levs(1) % state % circulationBtr % array(vertex1) &
+ - block % mesh % dcEdge % array (iEdge) &
+ *block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+
+ block % state % time_levs(1) % state % circulationBtr % array(vertex2) &
+ = block % state % time_levs(1) % state % circulationBtr % array(vertex2) &
+ + block % mesh % dcEdge % array (iEdge) &
+ *block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+ end do
+ do iVertex=1,block % mesh % nVertices
+ block % state % time_levs(1) % state % vorticityBtr % array(iVertex) &
+ = block % state % time_levs(1) % state % circulationBtr % array(iVertex) / block % mesh % areaTriangle % array (iVertex)
+ end do
+
+ !
+ ! Compute the divergence at each cell center
+ !
+ block % state % time_levs(1) % state % divergenceBtr % array(:) = 0.0
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+ block % state % time_levs(1) % state % divergenceBtr % array (cell1) &
+ = block % state % time_levs(1) % state % divergenceBtr % array (cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ *block % mesh % dvEdge % array(iEdge)
+
+ block % state % time_levs(1) % state % divergenceBtr % array (cell2) &
+ = block % state % time_levs(1) % state % divergenceBtr % array (cell2) &
+ - block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ *block % mesh % dvEdge % array(iEdge)
+ end do
+ do iCell = 1,block % mesh % nCells
+ block % state % time_levs(1) % state % divergenceBtr % array(iCell) &
+ = block % state % time_levs(1) % state % divergenceBtr % array(iCell) &
+ /block % mesh % areaCell % array(iCell)
+ enddo
+
+ !
+ ! Compute Btr diffusion
+ !
+ do iEdge=1,block % mesh % nEdgesSolve
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+ vertex1 = block % mesh % verticesOnEdge % array(1,iEdge)
+ vertex2 = block % mesh % verticesOnEdge % array(2,iEdge)
+
+ ! Here -( vorticityBtr(vertex2) - vorticityBtr(vertex1) ) / dvEdge % array (iEdge)
+ ! is - </font>
<font color="blue">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
+ ! + k \times </font>
<font color="gray">abla vorticity pointing from cell1 to cell2.
+
+ block % state % time_levs(1) % state % u_diffusionBtr % array(iEdge) = block % mesh % meshScalingDel2 % array (iEdge) * config_btr_mom_eddy_visc2 * &
+ (( block % state % time_levs(1) % state % divergenceBtr % array(cell2) - block % state % time_levs(1) % state % divergenceBtr % array(cell1) ) / block % mesh % dcEdge % array (iEdge) &
+ -( block % state % time_levs(1) % state % vorticityBtr % array(vertex2) - block % state % time_levs(1) % state % vorticityBtr % array(vertex1) ) / block % mesh % dvEdge % array (iEdge))
+
+ end do
+ end if
block => block % next
end do ! block
- ! compute new barotropic velocity
-! do iEdge=1,nEdges
-! cell1 = cellsOnEdge(1,iEdge)
-! cell2 = cellsOnEdge(2,iEdge)
-! grad_ssh = - gravity*rho0Inv*( sshSubcycleNew(cell2) &
-! - sshSubcycleNew(cell1) )/dcEdge(iEdge)
-! uBtrNew(iEdge) = uBtrNew(iEdge) + uBtrSubcycleNew(iEdge)
-! end do
+
+
+
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Final solve for velocity. Iterate for Coriolis term.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
do BtrCorIter=1,config_n_btr_cor_iter
- if (BtrCorIter==1) then
- uPerpTime = oldBtrSubcycleTime
- else
- uPerpTime = newBtrSubcycleTime
- endif
+ uPerpTime = newBtrSubcycleTime
- block => domain % blocklist
- do while (associated(block))
+ block => domain % blocklist
+ do while (associated(block))
do iEdge=1,block % mesh % nEdges
@@ -777,51 +940,61 @@
* block % mesh % fEdge % array(eoe)
end do
- ! timestep in uBtrSubcycle:
-! \ubtr_{n+(j+1)/J} = \ubtr_{n+j/J} + \frac{\Delta t}{J} \left(
-! - f\ubtr_{n+j/J}^{\perp}
-! - \frac{g}{\rho_0}</font>
<font color="gray">abla \zeta_{n+(j+1)/J} + {\overline {\bf G}}
-! \right),
-! uBtrSubCycleNew(iEdge) = uBtrSubcycleOld(iEdge) + dt*(fEdge(iEdge)*uBtrPerp(iEdge)
-! - grad_ssh + GBtrForcing(iEdge))
-
! mrp 110606 efficiency note: could make this a 1D integer factor instead of an if statement.
if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) = 0.0
else
+ ! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ dt/config_n_btr_subcycles *( &
uPerp &
- - grav_rho0Inv &
+ - gravity &
*( block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
/block % mesh % dcEdge % array(iEdge) &
- + block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+ + block % state % time_levs(1) % state % GBtrForcing % array(iEdge) &
+ + block % state % time_levs(1) % state % u_diffusionBtr % array(iEdge))
+ ! added del2 diffusion to btr solve
endif
end do
+ ! Implicit solve for barotropic momentum decay
+ if ( config_btr_mom_decay) then
+ ! Add term to RHS of momentum equation: -1/gamma u
+ !
+ ! This changes the solve to:
+ ! u^{n+1} = u_provis^{n+1}/(1+dt/gamma)
+ !
+ coef = 1.0/(1.0 + dt/config_n_btr_subcycles/config_btr_mom_decay_time)
+ do iEdge=1,block % mesh % nEdges
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ * coef
+ end do
+ endif
+
block => block % next
end do ! block
-! boundary update on \ubtr_{n+(j+1)/J}
+ ! boundary update on uBtrNew
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field1dReal(domain % dminfo, &
- block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:), &
- block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:), &
+ block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
block => block % next
end do ! block
- end do !do BtrCorIter=1,2
+ end do !do BtrCorIter=1,config_n_btr_cor_iter
! uBtrNew = uBtrNew + uBtrSubcycleNEW
@@ -887,31 +1060,40 @@
end do ! block
- ! Check that \zeta_{n} + \Delta t \left( - </font>
<font color="gray">abla \cdot {\overline {\bf F}} \right)
+ ! Check that you can compute SSH using the total sum or the individual increments
+ ! over the barotropic subcycles.
+ ! efficiency: This next block of code is really a check for debugging, and can
+ ! be removed later.
block => domain % blocklist
do while (associated(block))
allocate(uTemp(block % mesh % nVertLevels))
-! allocate(tend_ssh(block % mesh % nCells))
-
+ ! Accumulate fluxes in the tend % ssh variable
block % tend % ssh % array(:) = 0.0
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - block % state % time_levs(1) % state % FBtr % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + block % state % time_levs(1) % state % FBtr % array(iEdge)
+ block % tend % ssh % array(cell1) &
+ = block % tend % ssh % array(cell1) &
+ - block % state % time_levs(1) % state % FBtr % array(iEdge)
+ block % tend % ssh % array(cell2) &
+ = block % tend % ssh % array(cell2) &
+ + block % state % time_levs(1) % state % FBtr % array(iEdge)
+
end do
do iCell=1,block % mesh % nCells
-
+
+ ! SSHnew = SSHold + dt*(-div(Flux))
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell) &
+ dt &
* block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
end do
+
! Now can compare sshSubcycleNEW (big step using summed fluxes) with
! sshSubcycleOLD (individual steps to get there)
!print *, 'ssh, by substeps',minval(block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve)), &
@@ -919,10 +1101,11 @@
!print *, 'ssh, by 1 step ',minval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve)), &
! maxval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve))
- ! Correction velocity, u^{corr}:
-!u^{corr} = \left( {\overline {\bf F}}
-! - \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) u_k^* \right)
-!\left/ \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) \right.
+ ! Correction velocity uCorr = (Flux - Sum(h u*))/H
+ ! or, for the full latex version:
+ !u^{corr} = \left( {\overline {\bf F}}
+ ! - \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) u_k^* \right)
+ !\left/ \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) \right.
if (config_u_correction) then
ucorr_coef = 1
@@ -930,14 +1113,6 @@
ucorr_coef = 0
endif
-
-!print *, 'uBtr ',minval(block % state % time_levs(2) % state % uBtr % array(1:block % mesh % nEdges)), &
-! maxval(block % state % time_levs(2) % state % uBtr % array(1:block % mesh % nEdges))
-!print *, 'ssh ',minval(block % state % time_levs(2) % state % ssh % array(1:block % mesh % nEdges)), &
-! maxval(block % state % time_levs(2) % state % ssh % array(1:block % mesh % nEdges))
-!print *, 'FBtr ',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdges)), &
-! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdges))
-
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -962,8 +1137,6 @@
uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) &
/block % mesh % dvEdge % array(iEdge) &
- uhSum)/hSum)
-! mrp temp for printing
-!block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = uCorr
! put u^{tr}, the velocity for tracer transport, in uNew
! mrp 060611 not sure if boundary enforcement is needed here.
@@ -973,7 +1146,7 @@
block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
endif
- ! Put new SSH values in h array, for the compute_scalar_tend call below.
+ ! Put new sshEdge values in h_edge array, for the compute_scalar_tend call below.
block % state % time_levs(2) % state % h_edge % array(1,iEdge) &
= sshEdge + block % mesh % hZLevel % array(1)
@@ -984,9 +1157,6 @@
end do ! iEdge
-!print *, 'uCorr ',minval(block % state % time_levs(1) % state % GBtrForcing % array(1:block % mesh % nEdgesSolve)), &
-! maxval(block % state % time_levs(1) % state % GBtrForcing % array(1:block % mesh % nEdgesSolve))
-
! Put new SSH values in h array, for the compute_scalar_tend call below.
do iCell=1,block % mesh % nCells
block % state % time_levs(2) % state % h % array(1,iCell) &
@@ -1001,9 +1171,6 @@
enddo
enddo ! iCell
-
-
-! deallocate(tend_ssh)
deallocate(uTemp)
block => block % next
@@ -1012,11 +1179,11 @@
endif ! higdon_split
-
-
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stage 3: Tracer, density, pressure, vertical velocity prediction
!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
block => domain % blocklist
do while (associated(block))
@@ -1032,7 +1199,7 @@
block => block % next
end do
-! --- update halos for prognostic variables
+ ! --- update halos for prognostic variables
block => domain % blocklist
do while (associated(block))
@@ -1131,17 +1298,12 @@
do k=1,block % mesh % maxLevelCell % array(iCell)
do i=1,block % state % time_levs(1) % state % num_tracers
! This is Phi at n+1
-! mrp temp orig
block % state % time_levs(2) % state % tracers % array(i,k,iCell) &
= ( block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
* block % state % time_levs(1) % state % h % array(k,iCell) &
+ dt * block % tend % tracers % array(i,k,iCell) &
) / hNew(k)
-! mrp replace with this for constant tracers:
-! block % state % time_levs(2) % state % tracers % array(i,k,iCell) &
-! = block % state % time_levs(1) % state % tracers % array(i,k,iCell)
-
enddo
end do
end do
@@ -1163,48 +1325,23 @@
block => domain % blocklist
do while (associated(block))
-! nCells = block % mesh % nCells
-! nEdges = block % mesh % nEdges
-! nVertLevels = block % mesh % nVertLevels
-
-! printing:
-!print *, 'uold ',minval(block % state % time_levs(1) % state % u % array(:,1:nEdges)),&
-! maxval(block % state % time_levs(1) % state % u % array(:,1:nEdges))
-!print *, 'hold ',minval(block % state % time_levs(1) % state % h % array(:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % h % array(:,1:nCells))
-!print *, 'hold1 ',minval(block % state % time_levs(1) % state % h % array(1,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % h % array(1,1:nCells))
-!print *, 'Told ',minval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells))
-!print *, 'Sold ',minval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells))
-!print *, 'unew ',minval(block % state % time_levs(2) % state % u % array(:,1:nEdges)),&
-! maxval(block % state % time_levs(2) % state % u % array(:,1:nEdges))
-!print *, 'hnew ',minval(block % state % time_levs(2) % state % h % array(:,1:nCells)),&
-! maxval(block % state % time_levs(2) % state % h % array(:,1:nCells))
-!print *, 'hnew1 ',minval(block % state % time_levs(2) % state % h % array(1,1:nCells)),&
-! maxval(block % state % time_levs(2) % state % h % array(1,1:nCells))
-!print *, 'Tnew ',minval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells)),&
-! maxval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells))
-!print *, 'Snew ',minval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells)),&
-! maxval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells))
-!print *, 'Tr1Old1',minval(block % state % time_levs(1) % state % tracers % array(3,1,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(3,1,1:nCells))
-!print *, 'Tr1New1',minval(block % state % time_levs(2) % state % tracers % array(3,1,1:nCells)),&
-! maxval(block % state % time_levs(2) % state % tracers % array(3,1,1:nCells))
-! block % state % time_levs(2) % state % tracers % array(3,1,1:nCells) = domain % dminfo % my_proc_id
-! printing end
-
-
! Recompute final u to go on to next step.
! u_{n+1} = uBtr_{n+1} + uBcl_{n+1}
! Right now uBclNew is at time n+1/2, so back compute to get uBcl at time n+1
! using uBcl_{n+1/2} = 1/2*(uBcl_n + u_Bcl_{n+1})
! so the following lines are
! u_{n+1} = uBtr_{n+1} + 2*uBcl_{n+1/2} - uBcl_n
+ ! note that uBcl is recomputed at the beginning of the next timestep due to Imp Vert mixing,
+ ! so uBcl does not have to be recomputed here.
do iEdge=1,block % mesh % nEdges
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+
+ ! uBtrNew = uBtrSubcycleNew (old here is because counter already flipped)
+ ! This line is not needed if u is resplit at the beginning of the timestep.
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+
block % state % time_levs(2) % state % u % array(k,iEdge) &
= block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+2*block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
@@ -1218,7 +1355,6 @@
enddo ! iEdges
-
! do iEdge=1,block % mesh % nEdges
! block % state % time_levs(2) % state % u % array(:,iEdge) &
! = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
@@ -1246,6 +1382,12 @@
end do ! iCell
end if ! higdon_split
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Implicit vertical mixing, done after timestep is complete
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
u => block % state % time_levs(2) % state % u % array
tracers => block % state % time_levs(2) % state % tracers % array
h => block % state % time_levs(2) % state % h % array
@@ -1256,14 +1398,6 @@
maxLevelCell => block % mesh % maxLevelCell % array
maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
-
-! dividing by h above for higdon.
-! do iCell=1,nCells
-! do k=1,maxLevelCell(iCell)
-! tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
-! end do
-! end do
-
if (config_implicit_vertical_mix) then
allocate(A(block % mesh % nVertLevels),C(block % mesh % nVertLevels),uTemp(block % mesh % nVertLevels), &
tracersTemp(num_tracers,block % mesh % nVertLevels))
@@ -1337,6 +1471,26 @@
deallocate(A,C,uTemp,tracersTemp)
end if
+ ! mrp 110725 adding momentum decay term
+ if (config_mom_decay) then
+
+ !
+ ! Implicit solve for momentum decay
+ !
+ ! Add term to RHS of momentum equation: -1/gamma u
+ !
+ ! This changes the solve to:
+ ! u^{n+1} = u_provis^{n+1}/(1+dt/gamma)
+ !
+ coef = 1.0/(1.0 + dt/config_mom_decay_time)
+ do iEdge=1,block % mesh % nEdges
+ do k=1,maxLevelEdgeTop(iEdge)
+ u(k,iEdge) = coef*u(k,iEdge)
+ end do
+ end do
+
+ end if
+
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
@@ -1625,10 +1779,11 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- + q &
- - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + q &
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+
end do
end do
@@ -1651,11 +1806,11 @@
/ (zMidZLevel(k-1) - zMidZLevel(k))
end do
w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0
-
! Average w*du/dz from vertical interface to vertical middle of cell
do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
enddo
enddo
deallocate(w_dudzTopEdge)
@@ -1679,10 +1834,12 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( pressure(k,cell2) &
- - pressure(k,cell1) )/dcEdge(iEdge)
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - rho0Inv*( pressure(k,cell2) &
+ - pressure(k,cell1) )/dcEdge(iEdge)
end do
+
enddo
endif
@@ -1704,8 +1861,13 @@
! 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)
+! mrp 110816 orig
+! u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+! -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+! mrp 110816 remove vorticity
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) ! &
+! -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+! mrp 110816 end
u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
@@ -1738,9 +1900,15 @@
do k=1,maxLevelEdgeTop(iEdge)
+! mrp 110816 orig
+! delsq_u(k,iEdge) = &
+! ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+! -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+! mrp 110816 remove vort
delsq_u(k,iEdge) = &
- ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+ ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) ! &
+! -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+! mrp 110816 end
end do
end do
@@ -1792,11 +1960,21 @@
vertex2 = verticesOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
+ delsq_u(k,iEdge) = &
+ ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+! mrp 110816 orig
+! u_diffusion = ( delsq_divergence(k,cell2) &
+! - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+! -( delsq_vorticity(k,vertex2) &
+! - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+! mrp 110816 remove vort
u_diffusion = ( delsq_divergence(k,cell2) &
- - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+ - delsq_divergence(k,cell1) ) / dcEdge(iEdge) ! &
+! -( delsq_vorticity(k,vertex2) &
+! - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+! mrp 110816 end
u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
@@ -1873,10 +2051,245 @@
end subroutine compute_tend_u
- ! Put f*uBcl^{perp} in u as a work variable
+ subroutine filter_btr_mode_tend_u(tend, s, d, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Filter and remove barotropic mode from the tendencies
+ !
+ ! 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(in) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND) :: vertSum, uhSum, hSum, sshEdge
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
+ tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+! mrp 110516 cleanup fvertex fedge not used in this subroutine
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ tend_u => tend % u % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ u_src => grid % u_src % array
+
+ do iEdge=1,grid % nEdges
+
+ ! I am using hZLevel here. This assumes that SSH is zero everywhere already,
+ ! which should be the case if the barotropic mode is filtered.
+ ! The more general case is to use sshEdge or h_edge.
+ uhSum = (grid % hZLevel % array(1)) * tend_u(1,iEdge)
+ hSum = grid % hZLevel % array(1)
+
+ do k=2,grid % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum + grid % hZLevel % array(k) * tend_u(k,iEdge)
+ hSum = hSum + grid % hZLevel % array(k)
+ enddo
+
+ vertSum = uhSum/hSum
+
+ do k=1,grid % maxLevelEdgeTop % array(iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
+ enddo
+
+ enddo ! iEdge
+
+ end subroutine filter_btr_mode_tend_u
+
+ subroutine filter_btr_mode_u(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Filter and remove barotropic mode.
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND) :: vertSum, uhSum, hSum, sshEdge
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
+ tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+! mrp 110516 cleanup fvertex fedge not used in this subroutine
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ u_src => grid % u_src % array
+
+ do iEdge=1,grid % nEdges
+
+ ! I am using hZLevel here. This assumes that SSH is zero everywhere already,
+ ! which should be the case if the barotropic mode is filtered.
+ ! The more general case is to use sshedge or h_edge.
+ uhSum = (grid % hZLevel % array(1)) * u(1,iEdge)
+ hSum = grid % hZLevel % array(1)
+
+ do k=2,grid % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum + grid % hZLevel % array(k) * u(k,iEdge)
+ hSum = hSum + grid % hZLevel % array(k)
+ enddo
+
+ vertSum = uhSum/hSum
+ do k=1,grid % maxLevelEdgeTop % array(iEdge)
+ u(k,iEdge) = u(k,iEdge) - vertSum
+ enddo
+
+ enddo ! iEdge
+
+ end subroutine filter_btr_mode_u
+
subroutine compute_f_uperp(s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute height and normal wind tendencies, as well as diagnostic variables
+ ! Put f*uBcl^{perp} in u as a work variable
!
! Input: s - current model state
! grid - grid metadata
@@ -2574,18 +2987,6 @@
endif
10 format(2i8,10e20.10)
-
- ! print some diagnostics - for debugging
-! print *, 'after vertical mixing',&
-! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
-! do iTracer=1,num_tracers
-! do k = 1,nVertLevels
-! print '(2i5,20es12.4)', iTracer,k, &
-! minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))
-! enddo
-! enddo
-
-
end subroutine compute_scalar_tend
</font>
</pre>