<p><b>ringler@lanl.gov</b> 2011-07-27 16:14:06 -0600 (Wed, 27 Jul 2011)</p><p><br>
first cut at forcing module for 2D turblence studies.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-07-25 11:47:14 UTC (rev 930)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-07-27 22:14:06 UTC (rev 931)
@@ -24,6 +24,11 @@
namelist real sw_model config_alpha 0.0
namelist logical sw_model config_alpha_closure false
namelist real sw_model config_dragCoefficient 0.001
+namelist real sw_model config_enstrophy_injection_rate 0.0
+namelist real sw_model config_enstrophy_forcing_amp 0.0
+namelist real sw_model config_Lx 1.0
+namelist real sw_model config_Ly 1.0
+namelist logical sw_model config_forcing false
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-07-25 11:47:14 UTC (rev 930)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-07-27 22:14:06 UTC (rev 931)
@@ -79,6 +79,8 @@
call reconstructGradKE(block % state % time_levs(1) % state, mesh)
call compute_drag(block % state % time_levs(1) % state, mesh)
+ call compute_forcing(block % state % time_levs(1) % state, mesh)
+
if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-07-25 11:47:14 UTC (rev 930)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-07-27 22:14:06 UTC (rev 931)
@@ -7,7 +7,7 @@
implicit none
- real (KIND=RKIND), parameter :: threshold = 1.0e-8
+ real (KIND=RKIND), parameter :: threshold = 1.0e-9
integer, parameter :: nIterPerGS = 2, nCheckPoisson = 25, nQuitPoisson = 5000
private
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-07-25 11:47:14 UTC (rev 930)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-07-27 22:14:06 UTC (rev 931)
@@ -238,6 +238,7 @@
call reconstruct(provis, block % mesh)
call reconstructGradKE(provis, block % mesh)
call compute_drag(provis, block % mesh)
+ call compute_forcing(provis, block % mesh)
block => block % next
end do
@@ -327,6 +328,7 @@
call reconstruct(block % state % time_levs(2) % state, block % mesh)
call reconstructGradKE(block % state % time_levs(2) % state, block % mesh)
call compute_drag(block % state % time_levs(2) % state, block % mesh)
+ call compute_forcing(provis, block % mesh)
block => block % next
end do
@@ -457,7 +459,8 @@
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
- real (kind=RKIND), dimension(:,:), pointer :: tend_pv_advection, tend_pv_dissipation, tend_pv_forcing, tend_pv_drag, pv_drag
+ real (kind=RKIND), dimension(:,:), pointer :: tend_pv_advection, tend_pv_dissipation, tend_pv_forcing, tend_pv_drag
+ real (kind=RKIND), dimension(:,:), pointer :: pv_drag, pv_forcing, pv_dissipation
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -466,6 +469,8 @@
u => s % u % array
h_edge => s % h_edge % array
pv_drag => s % pv_drag % array
+ pv_forcing => s % pv_forcing % array
+ pv_dissipation => s % pv_dissipation % array
tracer_tend => tend % tracers % array
tend_pv_advection => tend % pv_advection % array
@@ -730,10 +735,14 @@
enddo
enddo
+ ! add forcing here
+ iTracer=1
+ do iCell = 1,grid%nCells
+ do k=1,grid % nVertLevels
+ tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + pv_forcing(k,iCell)
+ enddo
+ enddo
- ! compute forcing here
-
-
! store PV forcing term
do iCell = 1,grid%nCells
do k=1,grid % nVertLevels
@@ -857,6 +866,75 @@
end subroutine compute_drag
+ subroutine compute_forcing(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! forcing on RHS of PV equation
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), pointer, dimension(:) :: xCell, yCell, zCell
+ real (kind=RKIND), pointer, dimension(:,:) :: pv_forcing, h
+ real (kind=RKIND) :: wx, wy, kwave, lwave, phasex, phasey, rtime, pii, amp
+ integer :: iCell, k
+
+ pii = 4.0*atan(1.0)
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+
+ h => s % h % array
+ pv_forcing => s % pv_forcing % array
+ pv_forcing(:,:) = 0.0
+
+ if(config_forcing) then
+
+ write(6,*) ' inside compute forcing'
+
+! following lines need to be removed
+
+ rtime = 0.0
+ amp = 1.0e-12
+
+! above lines need to be removed
+
+ wx = 2.0*pii / (10.0*86400.0)
+ wy = wx * 3.0 / pii
+
+ phasex = pii * sin(wx*rtime)
+ phasey = pii * sin(wy*rtime)
+
+ kwave = 2.0*pii*4.0 / config_Lx
+ lwave = 2.0*pii*4.0 / config_Ly
+
+ write(6,*) kwave, lwave
+ write(6,*) phasex, phasey
+ write(6,*) wx, wy
+ write(6,*) config_enstrophy_forcing_amp, amp
+ write(6,*) pii
+ write(6,*)
+
+ do iCell = 1,grid%nCells
+ do k=1,grid % nVertLevels
+ pv_forcing(k,iCell) = - amp * cos( kwave*xCell(iCell) + phasex)
+ pv_forcing(k,iCell) = pv_forcing(k,iCell) + amp * cos( lwave*yCell(iCell) + phasey)
+ pv_forcing(k,iCell) = pv_forcing(k,iCell) / h(k,iCell)
+ enddo
+ enddo
+
+ endif
+
+ end subroutine compute_forcing
+
+
subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
</font>
</pre>