<p><b>mpetersen@lanl.gov</b> 2010-03-22 10:58:28 -0600 (Mon, 22 Mar 2010)</p><p>Converted density to a 3D array.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- trunk/mpas/src/core_ocean/module_test_cases.F 2010-03-19 03:47:18 UTC (rev 149)
+++ trunk/mpas/src/core_ocean/module_test_cases.F 2010-03-22 16:58:28 UTC (rev 150)
@@ -9,13 +9,13 @@
subroutine setup_sw_test_case(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Configure grid metadata and model state for the shallow water test case
! specified in the namelist
!
! Output: block - a subset (not necessarily proper) of the model domain to be
! initialized
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -23,8 +23,7 @@
integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src
- real (kind=RKIND), dimension(:), pointer :: rho, alpha
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND) :: delta_rho
integer :: nCells, nEdges, nVertices, nVertLevels
@@ -32,60 +31,42 @@
write(0,*) 'Using initial conditions supplied in input file'
else if (config_test_case == 1) then
- write(0,*) 'Setting up shallow water test case 1'
- write(0,*) ' -- Advection of Cosine Bell over the Pole'
+ write(0,*) ' Setting up shallow water test case 1:'
+ write(0,*) ' Advection of Cosine Bell over the Pole'
block_ptr => domain % blocklist
do while (associated(block_ptr))
call sw_test_case_1(block_ptr % mesh, block_ptr % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, block_ptr % time_levs(i) % state)
- end do
-
block_ptr => block_ptr % next
end do
else if (config_test_case == 2) then
- write(0,*) 'Setting up shallow water test case 2'
- write(0,*) ' -- Setup shallow water test case 2: '// &
+ write(0,*) ' Setup shallow water test case 2: '// &
'Global Steady State Nonlinear Zonal Geostrophic Flow'
block_ptr => domain % blocklist
do while (associated(block_ptr))
call sw_test_case_2(block_ptr % mesh, block_ptr % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, block_ptr % time_levs(i) % state)
- end do
-
block_ptr => block_ptr % next
end do
else if (config_test_case == 5) then
- write(0,*) 'Setting up shallow water test case 5'
- write(0,*) ' -- Setup shallow water test case 5:'// &
+ write(0,*) ' Setup shallow water test case 5:'// &
' Zonal Flow over an Isolated Mountain'
block_ptr => domain % blocklist
do while (associated(block_ptr))
call sw_test_case_5(block_ptr % mesh, block_ptr % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, block_ptr % time_levs(i) % state)
- end do
-
block_ptr => block_ptr % next
end do
else if (config_test_case == 6) then
- write(0,*) 'Setting up shallow water test case 6'
- write(0,*) ' -- Rossby-Haurwitz Wave'
+ write(0,*) ' Set up shallow water test case 6:'
+ write(0,*) ' Rossby-Haurwitz Wave'
block_ptr => domain % blocklist
do while (associated(block_ptr))
call sw_test_case_6(block_ptr % mesh, block_ptr % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, block_ptr % time_levs(i) % state)
- end do
-
block_ptr => block_ptr % next
end do
@@ -105,10 +86,9 @@
do while (associated(block_ptr))
h => block_ptr % time_levs(1) % state % h % array
u => block_ptr % time_levs(1) % state % u % array
+ rho => block_ptr % time_levs(1) % state % rho % array
u_src => block_ptr % mesh % u_src % array
- rho => block_ptr % mesh % rho % array
- alpha => block_ptr % mesh % alpha % array
nCells = block_ptr % mesh % nCells
nEdges = block_ptr % mesh % nEdges
@@ -125,12 +105,14 @@
end do
! define the density for multiple layers
- delta_rho=0.1
+ delta_rho=0.0
do iLevel = 1,nVertLevels
- rho(iLevel) = delta_rho*(iLevel-1)
+ rho(iLevel,1) = delta_rho*(iLevel-1)
enddo
- rho = rho - sum(rho)/nVertLevels + 1000.0
- alpha = 1/rho
+ rho(:,1) = rho(:,1) - sum(rho(:,1))/nVertLevels + 1000.0
+ do iLevel = 1,nVertLevels
+ rho(iLevel,:) = rho(iLevel,1)
+ enddo
do i=2,nTimeLevs
call copy_state(block_ptr % time_levs(1) % state, &
@@ -144,7 +126,7 @@
' min h max h ',&
' min u_src max u_src '
do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel), &
+ print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
minval(u(iLevel,:)), maxval(u(iLevel,:)), &
minval(h(iLevel,:)), maxval(h(iLevel,:)), &
minval(u_src(iLevel,:)), maxval(u_src(iLevel,:))
@@ -158,13 +140,13 @@
subroutine sw_test_case_1(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup shallow water test case 1: Advection of Cosine Bell over the Pole
!
! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
! Approximations to the Shallow Water Equations in Spherical
! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -233,14 +215,14 @@
subroutine sw_test_case_2(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup shallow water test case 2: Global Steady State Nonlinear Zonal
! Geostrophic Flow
!
! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
! Approximations to the Shallow Water Equations in Spherical
! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -325,13 +307,13 @@
subroutine sw_test_case_5(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup shallow water test case 5: Zonal Flow over an Isolated Mountain
!
! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
! Approximations to the Shallow Water Equations in Spherical
! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -447,13 +429,13 @@
subroutine sw_test_case_6(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup shallow water test case 6: Rossby-Haurwitz Wave
!
! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
! Approximations to the Shallow Water Equations in Spherical
! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -519,10 +501,10 @@
real function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -538,9 +520,9 @@
real function AA(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! A, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -557,9 +539,9 @@
real function BB(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! B, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -575,9 +557,9 @@
real function CC(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! C, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
Modified: trunk/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2010-03-19 03:47:18 UTC (rev 149)
+++ trunk/mpas/src/core_ocean/module_time_integration.F 2010-03-22 16:58:28 UTC (rev 150)
@@ -13,14 +13,14 @@
subroutine timestep(domain, dt)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step
!
! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
! plus grid meta-data
! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
! model state advanced forward in time by dt seconds
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -55,7 +55,7 @@
subroutine rk4(domain, dt)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
! 4th order Runge-Kutta
!
@@ -63,7 +63,7 @@
! plus grid meta-data
! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
! model state advanced forward in time by dt seconds
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -237,14 +237,14 @@
subroutine compute_tend(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
! Input: s - current model state
! grid - grid metadata
!
! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -426,13 +426,13 @@
subroutine compute_scalar_tend(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
! grid - grid metadata
!
! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -471,13 +471,13 @@
subroutine compute_solve_diagnostics(dt, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
!
! Input: grid - grid metadata
!
! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -495,8 +495,8 @@
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
!mrp 100112:
real (kind=RKIND), dimension(:,:), pointer :: &
- zmid, zbot, pmid, pbot, MontPot
- real (kind=RKIND), dimension(:), pointer :: zSurface, rho, alpha
+ zmid, zbot, pmid, pbot, MontPot, rho
+ real (kind=RKIND), dimension(:), pointer :: zSurface
real (kind=RKIND) :: delta_p
character :: c1*6
!mrp 100112 end
@@ -523,6 +523,7 @@
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
!mrp 100112:
+ rho => s % rho % array
zmid => s % zmid % array
zbot => s % zbot % array
pmid => s % pmid % array
@@ -548,10 +549,6 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
- !mrp 100112:
- rho => grid % rho % array
- alpha => grid % alpha % array
- !mrp 100112 end
nCells = grid % nCells
nEdges = grid % nEdges
@@ -811,17 +808,17 @@
zSurface(iCell) = zbot(1,iCell) + h(1,iCell)
! assume pressure at the surface is zero for now.
- pmid(1,iCell) = 0.5*rho(1)*gravity* h(1,iCell) ! + psurf(iCell)
- pbot(1,iCell) = rho(1)*gravity* h(1,iCell) ! + psurf(iCell)
+ pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell) ! + psurf(iCell)
+ pbot(1,iCell) = rho(1,iCell)*gravity* h(1,iCell) ! + psurf(iCell)
MontPot(1,iCell) = gravity * zSurface(iCell)
do k=2,nVertLevels
- delta_p = rho(k)*gravity* h(k,iCell)
+ delta_p = rho(k,iCell)*gravity* h(k,iCell)
pmid(k,iCell) = pbot(k-1,iCell) + 0.5*delta_p
pbot(k,iCell) = pbot(k-1,iCell) + delta_p
- ! from delta M = p delta alpha
+ ! from delta M = p delta / rho
MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pbot(k-1,iCell)*(alpha(k) - alpha(k-1))
+ + pbot(k-1,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
end do
end do
!mrp 100112 end
Modified: trunk/mpas/src/core_ocean/module_vector_reconstruction.F
===================================================================
--- trunk/mpas/src/core_ocean/module_vector_reconstruction.F 2010-03-19 03:47:18 UTC (rev 149)
+++ trunk/mpas/src/core_ocean/module_vector_reconstruction.F 2010-03-22 16:58:28 UTC (rev 150)
@@ -11,13 +11,13 @@
contains
subroutine init_reconstruct(grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: reconstruct vector field at cell centers based on radial basis functions
!
! Input: grid meta data and vector component data residing at cell edges
!
! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
@@ -167,13 +167,13 @@
subroutine reconstruct(s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: reconstruct vector field at cell centers based on radial basis functions
!
! Input: grid meta data and vector component data residing at cell edges
!
! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
</font>
</pre>