<p><b>duda</b> 2011-08-12 11:36:59 -0600 (Fri, 12 Aug 2011)</p><p>BRANCH COMMIT<br>
<br>
Add curvature terms (from Sang-Hun).<br>
<br>
<br>
M src/core_nhyd_atmos/module_time_integration.F<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-08-12 17:35:36 UTC (rev 935)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-08-12 17:36:59 UTC (rev 936)
@@ -1965,6 +1965,12 @@
! logical, parameter :: debug = .true.
logical, parameter :: debug = .false.
+ !SHP-curvature
+ logical, parameter :: curvature = .true.
+ real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+ real (kind=RKIND) :: r_earth
+ real (kind=RKIND), dimension(:,:), pointer :: ur_cell, vr_cell
+
real (kind=RKIND), parameter :: c_s = 0.25
real (kind=RKIND), dimension( grid % nVertLevels ) :: d_diag, d_off_diag
real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
@@ -1985,6 +1991,11 @@
!-----------
+ !SHP-curvature
+ r_earth = a
+ ur_cell => diag % uReconstructZonal % array
+ vr_cell => diag % uReconstructMeridional % array
+
coef_3rd_order = config_coef_3rd_order
rho => s % rho % array
@@ -2430,7 +2441,26 @@
end do
end do
+ !SHP-curvature
+ if (curvature) then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 2.*omega_e*cos(grid % angleEdge % array(iEdge))*cos(grid % latEdge % array(iEdge)) &
+ *.25*(rw(k,cell1)+rw(k+1,cell1)+rw(k,cell2)+rw(k+1,cell2)) &
+ - u(k,iEdge)*.25*(rw(k+1,cell1)+rw(k,cell1)+rw(k,cell2)+rw(k+1,cell2))/r_earth
+ ! + .5*(rho(k,cell1)+rho(k,cell2))*v(k,iEdge)*tan(grid % latEdge % array(iEdge)) &
+ ! *(u(k,iEdge)*cos(grid % angleEdge % array(iEdge)) &
+ ! -v(k,iEdge)*sin(grid % angleEdge % array(iEdge)))/r_earth
+ end do
+ end do
+ end if
+
+
!----------- rhs for w
tend_w(:,:) = 0.
@@ -2735,6 +2765,22 @@
deallocate(qtot)
+ !SHP-curvature
+ if (curvature) then
+
+ do iCell = 1, grid % nCellsSolve
+ do k=2,nVertLevels
+ tend_w(k,iCell) = tend_w(k,iCell) &
+ + rho(k,iCell)*( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
+ +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth &
+ + 2.*omega_e*cos(grid % latCell % array(iCell))*rho(k,iCell) &
+ *(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))
+
+ end do
+ end do
+
+ end if
+
!----------- rhs for theta
tend_theta(:,:) = 0.
</font>
</pre>