<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 &gt; 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)*(  &amp;
+                                             (theta(k+1,iCell)-theta(k  ,iCell))/(zp-z0)                 &amp;
+                                            -(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)*(  &amp;
+                                          (theta(k+1,iCell)-theta(k  ,iCell))/(zp-z0)                 &amp;
+                                         -(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 &gt; 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)*(  &amp;
-                                        (theta(k+1,iCell)-theta(k  ,iCell))/(zp-z0)                 &amp;
-                                       -(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>