<p><b>laura@ucar.edu</b> 2011-04-07 16:27:39 -0600 (Thu, 07 Apr 2011)</p><p>corrections to subroutines advance_scalars and advance_scalars_mono: commented out lines that set scalars_tend to zero; corrected the claculation of scalar_new at the bottom of subroutine subroutine advance_scalars_mono<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-04-07 19:52:10 UTC (rev 787)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-04-07 22:27:39 UTC (rev 788)
@@ -9,6 +9,7 @@
#ifdef DO_PHYSICS
use module_driver_microphysics
use module_physics_todynamics
+ use module_physics_utilities
#endif
contains
@@ -459,6 +460,8 @@
end do
! if(debug) then
+ 101 format(' min, max scalar',i4,2(1x,e17.10))
+ write(0,*)
block => domain % blocklist
do while (associated(block))
scalar_min = 0.
@@ -490,7 +493,7 @@
scalar_max = max(scalar_max, block % state % time_levs(2) % state % scalars % array(iScalar,k,iCell))
enddo
enddo
- write(0,*) ' min, max scalar ',iScalar, scalar_min, scalar_max
+ write(0,101) iScalar,scalar_min,scalar_max
end do
block => block % next
@@ -1225,6 +1228,11 @@
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, scalar_turb_flux, z1,z2,z3,z4,zm,z0,zp
+#ifdef DO_PHYSICS
+ character(len=80):: errmess
+ real(kind=RKIND):: scalar_max,scalar_min
+#endif
+
real (kind=RKIND) :: flux3, flux4
real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
@@ -1274,7 +1282,19 @@
qv_init => grid % qv_init % array
zgrid => grid % zgrid % array
- scalar_tend = 0. ! testing purposes - we have no sources or sinks
+!#ifndef DO_PHYSICS
+! scalar_tend = 0. ! testing purposes - we have no sources or sinks
+!#else
+!testing:
+ write(0,*)
+ write(0,*) 'subroutine advance_scalars:'
+ write(0,*) 'max tend qv=',maxval(scalar_tend(tend%index_qv,:,:))
+ write(0,*) 'max tend qc=',maxval(scalar_tend(tend%index_qc,:,:))
+ write(0,*) 'max tend qr=',maxval(scalar_tend(tend%index_qr,:,:))
+ write(0,*) 'max tend qi=',maxval(scalar_tend(tend%index_qi,:,:))
+ write(0,*) 'max tend qs=',maxval(scalar_tend(tend%index_qs,:,:))
+ write(0,*)
+!#endif
!
! Runge Kutta integration, so we compute fluxes from scalar_new values, update starts form scalar_old
@@ -1380,6 +1400,9 @@
! horizontal mixing for scalars - we could combine this with transport...
+!ldf (04-05-2011): do not do horizontal or vertical mixing if advection of scalars is set to be monotic.
+ if(.not. config_monotonic) then
+
if ( h_theta_eddy_visc2 > 0.0 ) then
do iEdge=1,grid % nEdges
@@ -1465,6 +1488,8 @@
end do
end if
+
+ end if !config_monotonic
!
! vertical flux divergence
@@ -1517,6 +1542,17 @@
end do
end do
+#ifdef DO_PHYSICS
+ 101 format(' min, max scalar',i4,2(1x,e17.10))
+ !testing:
+ do iScalar=1,s_old%num_scalars
+ scalar_max = maxval(scalar_new(iScalar,:,:))
+ scalar_min = minval(scalar_new(iScalar,:,:))
+ write(0,101) iScalar,scalar_min,scalar_max
+ enddo
+#endif
+
+
end subroutine advance_scalars
@@ -1565,6 +1601,13 @@
real (kind=RKIND) :: flux3, flux4
real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+#ifdef DO_PHYSICS
+ character(len=80):: errmess
+ integer:: imax,imin
+ integer:: kmax,kmin
+ real(kind=RKIND):: scalar_max,scalar_min
+#endif
+
flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
@@ -1599,7 +1642,18 @@
nVertLevels = grid % nVertLevels
- scalar_tend = 0. ! testing purposes - we have no sources or sinks
+!#ifndef DO_PHYSICS
+! scalar_tend = 0. ! testing purposes - we have no sources or sinks
+!#else
+ write(0,*)
+ write(0,*) 'subroutine advance_scalars_mono:'
+ write(0,*) 'max tend qv=',maxval(scalar_tend(tend%index_qv,:,:))
+ write(0,*) 'max tend qc=',maxval(scalar_tend(tend%index_qc,:,:))
+ write(0,*) 'max tend qr=',maxval(scalar_tend(tend%index_qr,:,:))
+ write(0,*) 'max tend qi=',maxval(scalar_tend(tend%index_qi,:,:))
+ write(0,*) 'max tend qs=',maxval(scalar_tend(tend%index_qs,:,:))
+ write(0,*)
+!#endif
!
! Runge Kutta integration, so we compute fluxes from scalar_new values, update starts from scalar_old
@@ -1875,19 +1929,29 @@
end do
! decouple from mass
- if (k > 1) then
- do iCell=1,grid % nCells
- do iScalar=1,s_old % num_scalars
- s_update(iScalar,iCell,km1) = s_update(iScalar,iCell,km1) / h_new(k-1,iCell)
- end do
- end do
+! if (k > 1) then
+! do iCell=1,grid % nCells
+! do iScalar=1,s_old % num_scalars
+! s_update(iScalar,iCell,km1) = s_update(iScalar,iCell,km1) / h_new(k-1,iCell)
+! end do
+! end do
- do iCell=1,grid % nCells
- do iScalar=1,s_old % num_scalars
- scalar_new(iScalar,k-1,iCell) = s_update(iScalar,iCell,km1)
- end do
- end do
- end if
+! do iCell=1,grid % nCells
+! do iScalar=1,s_old % num_scalars
+! scalar_new(iScalar,k-1,iCell) = s_update(iScalar,iCell,km1)
+! end do
+! end do
+! end if
+!ldf begin:
+ if(k > 1) then
+ do iCell = 1,grid%nCells
+ do iScalar = 1,s_old%num_scalars
+ scalar_new(iScalar,k-1,iCell) = (s_update(iScalar,iCell,km1)+dt*scalar_tend(iScalar,k-1,iCell)) &
+ / h_new(k-1,iCell)
+ enddo
+ enddo
+ endif
+!ldf end.
ktmp = km1
km1 = km0
@@ -1895,12 +1959,32 @@
end do
- do iCell=1,grid % nCells
- do iScalar=1,s_old % num_scalars
- scalar_new(iScalar,grid % nVertLevels,iCell) = s_update(iScalar,iCell,km1) / h_new(grid%nVertLevels,iCell)
- end do
- end do
+! do iCell=1,grid % nCells
+! do iScalar=1,s_old % num_scalars
+! scalar_new(iScalar,grid % nVertLevels,iCell) = s_update(iScalar,iCell,km1) / h_new(grid%nVertLevels,iCell)
+! end do
+! end do
+!ldf begin:
+ do iCell = 1,grid%nCells
+ do iScalar = 1,s_old%num_scalars
+ scalar_new(iScalar,grid%nVertLevels,iCell) = (s_update(iScalar,iCell,km1)+dt*scalar_tend(iScalar,grid%nVertLevels,iCell)) &
+ / h_new(grid%nVertLevels,iCell)
+ enddo
+ enddo
+!ldf end.
+#ifdef DO_PHYSICS
+ 101 format(' min, max scalar ',i4,2(1x,e17.10))
+ 102 format(i6,i4,8(1x,e17.10))
+ !testing:
+ do iScalar=1,s_old%num_scalars
+ scalar_max = maxval(scalar_new(iScalar,:,:))
+ scalar_min = minval(scalar_new(iScalar,:,:))
+ write(0,101) iScalar,scalar_min,scalar_max
+ where(scalar_new(iScalar,:,:) .lt. 0.) scalar_new(iScalar,:,:) = 0.
+ enddo
+#endif
+
end subroutine advance_scalars_mono
!----
</font>
</pre>