<p><b>ringler@lanl.gov</b> 2013-03-21 12:48:18 -0600 (Thu, 21 Mar 2013)</p><p><br>
adding Euler forward time stepping to dissipation terms in theta equation<br>
</p><hr noshade><pre><font color="gray">Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F        2013-03-21 17:49:04 UTC (rev 2650)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F        2013-03-21 18:48:18 UTC (rev 2651)
@@ -932,6 +932,10 @@
tend_theta(:,:) = 0.
+! tdr
+ write(6,*) ' zeroing euler_tend_theta'
+ euler_tend_theta(:,:) = 0.0
+! tdr
!
! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
@@ -948,6 +952,20 @@
flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
tend_theta(k,cell1) = tend_theta(k,cell1) + flux
tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+
+ ! tdr
+ ! (factor of 3.0 is to compensate for dt=dt/3 on the first RK step)
+ if(doEuler) then
+ if(rk_step.eq.1) then
+ euler_tend_theta(k,cell1) = euler_tend_theta(k,cell1) + 3.0*flux
+ euler_tend_theta(k,cell2) = euler_tend_theta(k,cell2) - 3.0*flux
+ endif
+ else
+ tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ endif
+ ! tdr
+
end do
end do
@@ -986,8 +1004,19 @@
theta_turb_flux = meshScalingDel4(iEdge)*h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ ! tdr
+ ! (factor of 3.0 is to compensate for dt=dt/3 on the first RK step)
+ if(doEuler) then
+ if(rk_step.eq.1) then
+ euler_tend_theta(k,cell1) = euler_tend_theta(k,cell1) - 3.0*flux
+ euler_tend_theta(k,cell2) = euler_tend_theta(k,cell2) + 3.0*flux
+ endif
+ else
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ endif
+ ! tdr
+
end do
end do
@@ -996,7 +1025,40 @@
end if
+ !
+ ! vertical mixing for theta - 2nd order
+ !
+ if ( v_theta_eddy_visc2 > 0.0 ) then
+ do iCell = 1, grid % nCellsSolve
+ do k=2,nVertLevels-1
+ z1 = geopotential(k-1,iCell)/gravity
+ z2 = geopotential(k ,iCell)/gravity
+ z3 = geopotential(k+1,iCell)/gravity
+ z4 = geopotential(k+2,iCell)/gravity
+
+ zm = 0.5*(z1+z2)
+ z0 = 0.5*(z2+z3)
+ zp = 0.5*(z3+z4)
+
+ ! tdr
+ ! (factor of 3.0 is to compensate for dt=dt/3 on the first RK step)
+ if(doEuler) then
+ if(rk_step.eq.1) then
+ euler_tend_theta(k,iCell) = euler_tend_theta(k,iCell) + 3.0*v_theta_eddy_visc2*prandtl*h(k,iCell)*( &
+ (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
+ -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
+                 endif
+ else
+ tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*h(k,iCell)*( &
+ (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
+ -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
+ endif
+ end do
+ end do
+
+ end if
+
!
! horizontal advection for theta
!
@@ -1121,35 +1183,11 @@
!
- ! vertical mixing for theta - 2nd order
- !
- if ( v_theta_eddy_visc2 > 0.0 ) then
-
- do iCell = 1, grid % nCellsSolve
- do k=2,nVertLevels-1
- z1 = geopotential(k-1,iCell)/gravity
- z2 = geopotential(k ,iCell)/gravity
- z3 = geopotential(k+1,iCell)/gravity
- z4 = geopotential(k+2,iCell)/gravity
-
- zm = 0.5*(z1+z2)
- z0 = 0.5*(z2+z3)
- zp = 0.5*(z3+z4)
-
- tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*h(k,iCell)*( &
- (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
- -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
- end do
- end do
-
- end if
-
- !
! Add theta tendency from physics
!
do iCell=1,grid % nCells
do k=1,grid % nVertLevels
- tend_theta(k,iCell) = tend_theta(k,iCell) + theta_phys_tend(k,iCell)
+ tend_theta(k,iCell) = tend_theta(k,iCell) + theta_phys_tend(k,iCell) + euler_tend_theta(k,iCell)
end do
end do
</font>
</pre>