<p><b>ringler@lanl.gov</b> 2011-08-17 11:42:21 -0600 (Wed, 17 Aug 2011)</p><p><br>
updated forcing functions with normalization to provide constant pot. enstrophy flux into global domain.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/namelist.input.pvsw
===================================================================
--- branches/pv_based_swm/mpas/namelist.input.pvsw 2011-08-16 20:57:45 UTC (rev 943)
+++ branches/pv_based_swm/mpas/namelist.input.pvsw 2011-08-17 17:42:21 UTC (rev 944)
@@ -4,23 +4,27 @@
config_alpha = 0.0
config_test_case = 0
config_time_integration = 'RK4'
- config_dt = 2400
- config_ntimesteps = 36
- config_output_interval = 1
- config_stats_interval = 0
+ config_dt = 300
+ config_ntimesteps = 100000
+ config_output_interval = 1000
+ config_stats_interval = 120
config_h_mom_eddy_visc2 = 0.0
config_h_mom_eddy_visc4 = 0.0
- config_h_tracer_eddy_diff2 = 1260.0
+ config_h_tracer_eddy_diff2 = 88.0
config_h_tracer_eddy_diff4 = 0.0
config_thickness_adv_order = 2
config_tracer_adv_order = 4
config_positive_definite = .false.
config_monotonic = .false.
config_wind_stress = .false.
- config_bottom_drag = .false.
+ config_bottom_drag = .true.
config_bottom_drag_term1 = .true.
config_bottom_drag_term2 = .true.
- config_dragCoefficient = 2.0e-5
+ config_dragCoefficient = 5.0e-6
+ config_forcing = .true.
+ config_Lx = 10080000.0000000
+ config_Ly = 10080535.7000509
+ config_enstrophy_injection_rate = 1.75e-18
/
&io
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-08-16 20:57:45 UTC (rev 943)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-08-17 17:42:21 UTC (rev 944)
@@ -16,6 +16,7 @@
use test_cases
use poisson_solver
use alpha_solver
+ use time_integration
implicit none
@@ -23,11 +24,44 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
+ real (kind=RKIND) :: gVor, localArea, globalgVor, globalArea
write(6,*) 'calling sw_test_case'
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ ! -------------
+ ! needs to be moved out of this top level
+ !
+ ! find amplitude of forcing based on user-specified enstrophy injection rate
+ if(config_forcing) then
+ write(6,*) ' init: calling compute_injection_rate '
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_injection_rate(block % state % time_levs(1) % state, block % mesh, gVor, localArea)
+ block => block % next
+ enddo
+ write(6,*) ' init: calling dmpar_sum_real'
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_sum_real(domain % dminfo, gVor, globalgVor)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
+ block => block % next
+ end do
+ write(6,*) ' init: filling pvForcingAmp'
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % pvForcingAmp % scalar = config_enstrophy_injection_rate / (globalgVor / globalArea) ! units of s^-2
+ write(6,10) 'pvForcingAmp :',block % state % time_levs(1) % state % xtime % scalar, block % state % time_levs(1) % state % pvForcingAmp % scalar, globalgVor
+ 10 format(a15,3e20.10)
+ block => block % next
+ enddo
+ endif
+ !
+ ! needs to be moved out of this top level
+ !
+ ! -------------
+
!
! Initialize core
!
@@ -71,7 +105,7 @@
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
-
+
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
call rbfInterp_initialize(mesh)
@@ -127,7 +161,7 @@
call shift_time_levels_state(domain % blocklist % state)
call timer_start("model write")
- if (mod(itimestep, config_output_interval) == 0) then
+ if (mod(itimestep, config_output_interval) == 0 .or. itimestep == 1) then
call write_output_frame(output_obj, output_frame, domain)
end if
if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
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-08-16 20:57:45 UTC (rev 943)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-08-17 17:42:21 UTC (rev 944)
@@ -70,9 +70,36 @@
integer :: rk_step
+ real (kind=RKIND) :: gVor, localArea, globalgVor, globalArea
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
call timer_start("rk4 init")
+
+ ! find amplitude of forcing based on user-specified enstrophy injection rate
+ if(config_forcing) then
+ ! write(6,*) ' calling compute_injection_rate '
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_injection_rate(block % state % time_levs(1) % state, block % mesh, gVor, localArea)
+ block => block % next
+ enddo
+ ! write(6,*) ' calling dmpar_sum_real'
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_sum_real(domain % dminfo, gVor, globalgVor)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
+ block => block % next
+ end do
+ ! write(6,*) ' filling pvForcingAmp'
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % pvForcingAmp % scalar = config_enstrophy_injection_rate / (globalgVor / globalArea) ! units of s^-2
+ write(6,10) 'pvForcingAmp :',block % state % time_levs(1) % state % xtime % scalar, block % state % time_levs(1) % state % pvForcingAmp % scalar, globalgVor
+ 10 format(a15,3e20.10)
+ block => block % next
+ enddo
+ endif
+
block => domain % blocklist
call allocate_state(provis, &
block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
@@ -867,6 +894,65 @@
end subroutine compute_drag
+ subroutine compute_injection_rate(s, grid, gVor, localArea)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! compute \int (G * vor) dArea
+ !
+ ! 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), intent(out) :: gVor, localArea
+
+ real (kind=RKIND), pointer, dimension(:) :: areaCell, fCell
+ real (kind=RKIND), pointer, dimension(:) :: xCell, yCell, zCell
+ real (kind=RKIND), pointer, dimension(:,:) :: h, pv
+ real (kind=RKIND) :: r1, kwave, lwave, pii, phix, phiy, nuy, wx, rtime
+ integer :: iCell, k
+
+ write(6,*) ' inside compute_injection_rate '
+
+ pii = 4.0*atan(1.0)
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ areaCell => grid % areaCell % array
+ fCell => grid % fCell % array
+
+ h => s % h % array
+ pv => s % pv % array
+
+ rtime = s % xtime % scalar
+
+ kwave = 2.0*pii*4.0 / config_Lx
+ lwave = 2.0*pii*4.0 / config_Ly
+
+ nuy = 2e-7 ! 1/s (approx 1 / 231 days)
+ wx = 1.2e-6 ! 1/s
+
+ phiy = pii * sin (2 * pii * nuy * rtime)
+ phix = pii * sin (wx * rtime)
+
+ gVor = 0.0
+ localArea = 0.0
+ do iCell = 1, grid % nCellsSolve
+ do k=1,grid % nVertLevels
+ r1 = cos(lwave*yCell(iCell)+phiy) - cos(kwave*xCell(iCell)+phix)
+ gVor = gVor + r1 * (pv(k,iCell)*h(k,iCell) - fCell(iCell)) * areaCell(iCell)
+ localArea = localArea + areaCell(iCell)
+ enddo
+ enddo
+
+ write(6,*) ' exiting compute_injection_rate ', gVor
+
+ end subroutine compute_injection_rate
+
+
subroutine compute_forcing(s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! forcing on RHS of PV equation
@@ -883,7 +969,7 @@
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
+ real (kind=RKIND) :: wx, wy, kwave, lwave, phasex, phasey, rtime, pii, amp, nuy, phiy, phix
integer :: iCell, k
pii = 4.0*atan(1.0)
@@ -898,25 +984,22 @@
rtime = s % xtime % scalar
amp = s % pvForcingAmp % scalar
+ write(6,*) 'compute forcing : ',rtime, amp
- amp = 1.0e-8
-
if(config_forcing) then
- wx = 2.0 * pii / (10.0*86400.0)
- wy = 2.0 * 3.0 / (10.0*86400.0)
-
- 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
+ nuy = 2e-7 ! 1/s (approx 1 / 231 days)
+ wx = 1.2e-6 ! 1/s
+
+ phiy = pii * sin (2 * pii * nuy * rtime)
+ phix = pii * sin (wx * rtime)
+
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)
+ pv_forcing(k,iCell) = amp * (cos(lwave*yCell(iCell)+phiy) - cos(kwave*xCell(iCell)+phix)) / h(k,iCell)
enddo
enddo
</font>
</pre>