<p><b>mhoffman@lanl.gov</b> 2012-03-07 16:27:50 -0700 (Wed, 07 Mar 2012)</p><p>BRANCH COMMIT -- land_ice<br>
<br>
Cleanup of land ice core, including:<br>
* Removal of test cases module (test cases will be run from python I.C. setup scripts from now on)<br>
* Addition of needed state initialization to land_ice_mpas_core (that had previously been done in the test cases module)<br>
* Removal of unused parts of the time integration module<br>
* Symlinking namelist.input to the default land ice namelist<br>
* Removal of dimension nTracers from the Registry (not used)<br>
* Added plots of upper and lowerSurface to the dome viz script to allow validation that the needed state init occurred.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/namelist.input
===================================================================
--- branches/land_ice_projects/implement_core/namelist.input 2012-03-07 22:50:57 UTC (rev 1603)
+++ branches/land_ice_projects/implement_core/namelist.input 2012-03-07 23:27:50 UTC (rev 1604)
@@ -1 +1 @@
-link namelist.input.ocean
\ No newline at end of file
+link namelist.input.land_ice
\ No newline at end of file
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-03-07 22:50:57 UTC (rev 1603)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-03-07 23:27:50 UTC (rev 1604)
@@ -40,7 +40,6 @@
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nVertLevelsPlus2 nVertLevels+2
-dim nTracers nTracers
%
% var persistence type name_in_file ( dims ) time_levs iro- name_in_code struct super-array array_class
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-03-07 22:50:57 UTC (rev 1603)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-03-07 23:27:50 UTC (rev 1604)
@@ -34,8 +34,6 @@
err = 0
- if (.not. config_do_restart) call setup_land_ice_test_case(domain)
-
!
! Initialize core
!
@@ -133,9 +131,24 @@
type (dm_info) :: dminfo
integer :: err, err_tmp
+ integer :: iCell, iLevel, nCells, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
+ lowerSurface, bedTopography, layerThicknessFractions
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+
err = 0
- call land_ice_compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+ nCells = mesh % nCells
+ nVertLevels = mesh % nVertLevels
+ layerThicknessFractions => mesh % layerThicknessFractions % array
+ thickness => block % state % time_levs(1) % state % thickness % array
+ upperSurface => block % state % time_levs(1) % state % upperSurface % array
+ lowerSurface => block % state % time_levs(1) % state % lowerSurface % array
+ bedTopography => block % state % time_levs(1) % state % bedTopography % array
+ layerThickness => block % state % time_levs(1) % state % layerThickness % array
+
+
+
call compute_mesh_scaling(mesh)
call mpas_rbf_interp_initialize(mesh)
@@ -147,7 +160,25 @@
block % state % time_levs(1) % state % uReconstructZonal % array, &
block % state % time_levs(1) % state % uReconstructMeridional % array &
)
-
+
+
+ ! Initialize state ==== eventually this should be moved to its own subroutine/module
+ ! no ice shelves -- lower surface same as bed topography
+ ! Eventually this should consider floatation!
+ lowerSurface(:) = bedTopography(:)
+
+ ! as always, upper surface is the lower surface plus the thickness
+ upperSurface(:) = lowerSurface(:) + thickness(:)
+
+ ! set up layerThickness from thickness using the layerThicknessFractions
+ do iCell = 1, nCells
+ do iLevel = 1, nVertLevels
+ layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
+ end do
+ end do
+ ! end initalize state ====
+
+
call land_ice_vel_block_init(mesh, err_tmp)
err = ior(err, err_tmp)
Deleted: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_test_cases.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_test_cases.F 2012-03-07 22:50:57 UTC (rev 1603)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_test_cases.F 2012-03-07 23:27:50 UTC (rev 1604)
@@ -1,240 +0,0 @@
-module land_ice_test_cases
-
- use mpas_grid_types
- use mpas_configure
- use mpas_constants
-
- contains
-
-
- subroutine setup_land_ice_test_case(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Configure grid metadata and model state for the test case
- ! specified in the namelist
- !
- ! Output: block - a subset (not necessarily proper) of the model domain to be
- ! initialized
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i
- type (block_type), pointer :: block_ptr
-
- if (config_test_case == 0) then
- 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'
-
-!!! block_ptr => domain % blocklist
-!!! do while (associated(block_ptr))
-!!! call sw_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
-!!! do i=2,nTimeLevs
-!!! call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
-!!! end do
-
-!!! block_ptr => block_ptr % next
-!!! end do
-
-!whl - removed remaining shallow water test cases
- else if (config_test_case == 1) then
- ! the dome test case
- write(0,*) 'Setting up land ice test case 1'
- write(0,*) ' -- Evolution of an isothermal dome'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call land_ice_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else
- write(0,*) 'Only test case 0 is currently supported.'
- stop
- end if
-
- end subroutine setup_land_ice_test_case
-
-
- subroutine land_ice_test_case_1(mesh, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup land ice test case 1: Isothermal dome
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: mesh
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: r0 = 60000.0*sqrt(0.125) !meters
- real (kind=RKIND), parameter :: h0 = 2000.0*sqrt(0.125) !meters
- real (kind=RKIND), parameter :: x0 = 30000.0 !meters
- real (kind=RKIND), parameter :: y0 = 30000.0 !meters
-
- integer :: iCell, iLevel, nCells, nVertLevels, index_temperature
- real (kind=RKIND) :: r
-
- real (kind=RKIND), dimension(:), pointer :: xCell, yCell, thickness, upperSurface, &
- lowerSurface, bedTopography, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
-
-
- if(mesh % on_a_sphere) then
- print *, "Land ice test case 1 only supports planar meshes."
- stop
- end if
-
- nCells = mesh % nCells
- nVertLevels = mesh % nVertLevels
-
- xCell => mesh % xCell % array
- yCell => mesh % yCell % array
- layerThicknessFractions => mesh % layerThicknessFractions % array
- thickness => state % thickness % array
- upperSurface => state % upperSurface % array
- lowerSurface => state % lowerSurface % array
- bedTopography => state % bedTopography % array
- normalVelocity => state % normalVelocity % array
- layerThickness => state % layerThickness % array
- tracers => state % tracers % array
-
- index_temperature = state % index_temperature
-
- ! zero velocity everywhere
- normalVelocity(:,:) = 0.0
-
- ! flat bed at sea level
- bedTopography(:) = 0.0
-
- ! constant, arbitrary temperature
- tracers(index_temperature,:,:) = -20.0 ! degrees C
-
- do iCell = 1, nCells
- r = sqrt((xCell(iCell)-x0)**2 + (yCell(iCell)-y0)**2)
- if(r < r0) then
- thickness(iCell) = h0*sqrt(1 - (r/r0)**2)
- else
- thickness(iCell) = 0.0
- end if
- end do
-
- ! no ice shelves -- lower surface same as bed topography
- lowerSurface(:) = bedTopography(:)
-
- ! as always, upper surface is the lower surface plus the thickness
- upperSurface(:) = lowerSurface(:) + thickness(:)
-
- ! set up layerThickness from thickness using the layerThicknessFractions
- do iCell = 1, nCells
- do iLevel = 1, nVertLevels
- layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
- end do
- end do
-
-
- end subroutine land_ice_test_case_1
-
-!!! 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
-
-!!! type (mesh_type), intent(inout) :: grid
-!!! type (state_type), intent(inout) :: state
-
-!!! real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
-!!! real (kind=RKIND), parameter :: h0 = 1000.0
-!!! real (kind=RKIND), parameter :: theta_c = 0.0
-!!! real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
-!!! real (kind=RKIND), parameter :: alpha = pii/4.0
-
-!!! integer :: iCell, iEdge, iVtx
-!!! real (kind=RKIND) :: r, u, v
-!!! real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
-!!! grid % xCell % array = grid % xCell % array * a
-!!! grid % yCell % array = grid % yCell % array * a
-!!! grid % zCell % array = grid % zCell % array * a
-!!! grid % xVertex % array = grid % xVertex % array * a
-!!! grid % yVertex % array = grid % yVertex % array * a
-!!! grid % zVertex % array = grid % zVertex % array * a
-!!! grid % xEdge % array = grid % xEdge % array * a
-!!! grid % yEdge % array = grid % yEdge % array * a
-!!! grid % zEdge % array = grid % zEdge % array * a
-!!! grid % dvEdge % array = grid % dvEdge % array * a
-!!! grid % dcEdge % array = grid % dcEdge % array * a
-!!! grid % areaCell % array = grid % areaCell % array * a**2.0
-!!! grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
-!!! grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
-!!! allocate(psiVertex(grid % nVertices))
-!!! do iVtx=1,grid % nVertices
-!!! psiVertex(iVtx) = -a * u0 * ( &
-!!! sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
-!!! cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
-!!! )
-!!! end do
-!!! do iEdge=1,grid % nEdges
-!!! state % u % array(1,iEdge) = -1.0 * ( &
-!!! psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
-!!! psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
-!!! ) / grid%dvEdge%array(iEdge)
-!!! end do
-!!! deallocate(psiVertex)
-
- !
- ! Initialize cosine bell at (theta_c, lambda_c)
- !
-!!! do iCell=1,grid % nCells
-!!! r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
-!!! if (r < a/3.0) then
-!!! state % h % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
-!!! else
-!!! state % h % array(1,iCell) = 0.0
-!!! end if
-!!! end do
-
-!!! end subroutine sw_test_case_1
-
-
-!!! 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
-
-!!! real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
-!!! real (kind=RKIND) :: arg1
-
-!!! arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
-!!! cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
-!!! sphere_distance = 2.*radius*asin(arg1)
-
-!!! end function sphere_distance
-
-
-end module land_ice_test_cases
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-03-07 22:50:57 UTC (rev 1603)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-03-07 23:27:50 UTC (rev 1604)
@@ -69,1261 +69,4 @@
end subroutine land_ice_timestep
- subroutine land_ice_rk4(domain, dt)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Advance model state forward in time by the specified time step using
- ! 4th order Runge-Kutta
- !
- ! 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
-
- type (domain_type), intent(inout) :: domain
- real (kind=RKIND), intent(in) :: dt
-
- integer :: iCell, k
- type (block_type), pointer :: block
- type (state_type) :: provis
-
- integer :: rk_step
-
- real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
-
- block => domain % blocklist
- call mpas_allocate_state(provis, &
- block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
- block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels, &
- block % mesh % nTracers)
-
- !
- ! Initialize time_levs(2) with state at current time
- ! Initialize first RK state
- ! Couple tracers time_levs(2) with h in time-levels
- ! Initialize RK weights
- !
- block => domain % blocklist
- do while (associated(block))
-
- block % state % time_levs(2) % state % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:)
-
- block % state % time_levs(2) % state % layerThickness % array(:,:) = block % state % time_levs(1) % state % layerThickness % array(:,:)
-
- do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % layerThickness % array(k,iCell)
- end do
- end do
-
- call mpas_copy_state(provis, block % state % time_levs(1) % state)
-
- block => block % next
- end do
-
- rk_weights(1) = dt/6.
- rk_weights(2) = dt/3.
- rk_weights(3) = dt/3.
- rk_weights(4) = dt/6.
-
- rk_substep_weights(1) = dt/2.
- rk_substep_weights(2) = dt/2.
- rk_substep_weights(3) = dt
- rk_substep_weights(4) = 0.
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! BEGIN RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do rk_step = 1, 4
-
-! --- update halos for diagnostic variables
-
- block => domain % blocklist
- do while (associated(block))
-!!! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % pv_edge % array(:,:), &
-!!! block % mesh % nVertLevels, block % mesh % nEdges, &
-!!! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
-!!! if (config_h_mom_eddy_visc4 > 0.0) then
-!!! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % divergence % array(:,:), &
-!!! block % mesh % nVertLevels, block % mesh % nCells, &
-!!! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-!!! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % vorticity % array(:,:), &
-!!! block % mesh % nVertLevels, block % mesh % nVertices, &
-!!! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
-!!! end if
-
- block => block % next
- end do
-
-! --- compute tendencies
-
- block => domain % blocklist
- do while (associated(block))
- call land_ice_compute_tend(block % tend, provis, block % mesh)
- call land_ice_compute_scalar_tend(block % tend, provis, block % mesh)
- call land_ice_enforce_boundary_edge(block % tend, block % mesh)
- block => block % next
- end do
-
-! --- update halos for prognostic variables
-
- block => domain % blocklist
- do while (associated(block))
-! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % normalVelocity % array(:,:), &
-! block % mesh % nVertLevels, block % mesh % nEdges, &
-! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % layerThickness % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call mpas_dmpar_exch_halo_field3d_real(domain % dminfo, block % tend % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
-
-! --- compute next substep state
-
- if (rk_step < 4) then
- block => domain % blocklist
- do while (associated(block))
-! provis % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:) &
-! + rk_substep_weights(rk_step) * block % tend % normalVelocity % array(:,:)
- provis % layerThickness % array(:,:) = block % state % time_levs(1) % state % layerThickness % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % layerThickness % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- provis % tracers % array(:,k,iCell) = ( &
- block % state % time_levs(1) % state % layerThickness % array(k,iCell) * &
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % layerThickness % array(k,iCell)
- end do
- end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- provis % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:)
- end if
- call land_ice_compute_solve_diagnostics(dt, provis, block % mesh)
- block => block % next
- end do
- end if
-
-!--- accumulate update (for RK4)
-
- block => domain % blocklist
- do while (associated(block))
-! block % state % time_levs(2) % state % normalVelocity % array(:,:) = block % state % time_levs(2) % state % normalVelocity % array(:,:) &
-! + rk_weights(rk_step) * block % tend % normalVelocity % array(:,:)
- block % state % time_levs(2) % state % layerThickness % array(:,:) = block % state % time_levs(2) % state % layerThickness % array(:,:) &
- + rk_weights(rk_step) * block % tend % layerThickness % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
- end do
- end do
- block => block % next
- end do
-
- end do
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! END RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
- !
- ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
- !
- block => domain % blocklist
- do while (associated(block))
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- / block % state % time_levs(2) % state % layerThickness % array(k,iCell)
- end do
- end do
-
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % state % time_levs(2) % state % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:)
- end if
-
- call land_ice_compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
-
- call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % normalVelocity% array, &
- block % state % time_levs(2) % state % uReconstructX % array, &
- block % state % time_levs(2) % state % uReconstructY % array, &
- block % state % time_levs(2) % state % uReconstructZ % array, &
- block % state % time_levs(2) % state % uReconstructZonal % array, &
- block % state % time_levs(2) % state % uReconstructMeridional % array &
- )
-
- block => block % next
- end do
-
- call mpas_deallocate_state(provis)
-
- end subroutine land_ice_rk4
-
-
- subroutine land_ice_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
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
-
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, areaTriangle, &
- meshScalingDel2, meshScalingDel4
-!!! real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge
-!!! real (kind=RKIND), dimension(:,:), pointer :: vh
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, &
- layerThicknessEdge, layerThickness, normalVelocity, tend_layerThickness, &
-!!! layerThicknessEdge, h, u, v, tend_h, tend_u, &
-!!! circulation, vorticity, ke, pv_edge, divergence, &
- layerThicknessVertex
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r
-!!! real (kind=RKIND) :: u_diffusion
-
-!!! real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
-!!! real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
-!!! real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-!!! real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
- real (kind=RKIND) :: ke_edge
-
-
-!!! h => s % h % array
-!!! u => s % u % array
-!!! v => s % v % array
- layerThickness => s % layerThickness % array
- normalVelocity => s % normalVelocity % array
-
- layerThicknessEdge => s % layerThicknessEdge % array
-!!! circulation => s % circulation % array
-!!! vorticity => s % vorticity % array
-!!! divergence => s % divergence % array
-!!! ke => s % ke % array
-!!! pv_edge => s % pv_edge % array
-!!! vh => s % vh % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
-!!! h_s => grid % h_s % array
-!!! fVertex => grid % fVertex % array
-!!! fEdge => grid % fEdge % array
-
- tend_layerThickness => tend % layerThickness % array
-! tend_u => tend % u % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
-!!! u_src => grid % u_src % array
-
- meshScalingDel2 => grid % meshScalingDel2 % array
- meshScalingDel4 => grid % meshScalingDel4 % array
-
-
- !
- ! Compute height tendency for each cell
- !
- tend_layerThickness(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge)
- tend_layerThickness(k,cell1) = tend_layerThickness(k,cell1) - flux
- tend_layerThickness(k,cell2) = tend_layerThickness(k,cell2) + flux
- end do
- end do
- do iCell=1,grid % nCellsSolve
- do k=1,nVertLevels
- tend_layerThickness(k,iCell) = tend_layerThickness(k,iCell) / areaCell(iCell)
- end do
- end do
-
-#ifdef LANL_FORMULATION
- !
- ! Compute u (normal) velocity tendency for each edge (cell face)
- !
-!!! tend_u(:,:) = 0.0
-!!! do iEdge=1,grid % nEdgesSolve
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! vertex1 = verticesOnEdge(1,iEdge)
-!!! vertex2 = verticesOnEdge(2,iEdge)
-
-!!! do k=1,nVertLevels
-!!! q = 0.0
-!!! do j = 1,nEdgesOnEdge(iEdge)
-!!! eoe = edgesOnEdge(j,iEdge)
-!!! workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
-!!! q = q + weightsOnEdge(j,iEdge) * normalVelocity(k,eoe) * workpv * layerThicknessEdge(k,eoe)
-!!! end do
-
-!!! tend_normalVelocity(k,iEdge) = &
-!!! q &
-!!! - ( ke(k,cell2) - ke(k,cell1) + &
-!!! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
-!!! ) / dcEdge(iEdge)
-!!! end do
-!!! end do
-
-
-#endif
-
-#ifdef NCAR_FORMULATION
- !
- ! Compute u (normal) velocity tendency for each edge (cell face)
- !
-!!! tend_u(:,:) = 0.0
-!!! do iEdge=1,grid % nEdgesSolve
-!!! vertex1 = verticesOnEdge(1,iEdge)
-!!! vertex2 = verticesOnEdge(2,iEdge)
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-
-!!! do k=1,nVertLevels
-!!! vorticity_abs = fEdge(iEdge) + (circulation(k,vertex1) + circulation(k,vertex2)) / &
-!!! (areaTriangle(vertex1) + areaTriangle(vertex2))
-
-!!! workpv = 2.0 * vorticity_abs / (h(k,cell1) + h(k,cell2))
-
-!!! tend_normalVelocity(k,iEdge) = workpv * vh(k,iEdge) - &
-!!! (ke(k,cell2) - ke(k,cell1) + &
-!!! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
-!!! ) / &
-!!! dcEdge(iEdge)
-!!! end do
-!!! end do
-#endif
-
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! only valid for visc == constant
-!!! if (config_h_mom_eddy_visc2 > 0.0) then
-!!! do iEdge=1,grid % nEdgesSolve
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! vertex1 = verticesOnEdge(1,iEdge)
-!!! vertex2 = verticesOnEdge(2,iEdge)
-
-!!! do k=1,nVertLevels
-!!! u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-!!! -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
-!!! u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
-!!! tend_normalVelocity(k,iEdge) = tend_normalVelocity(k,iEdge) + u_diffusion
-!!! end do
-!!!! end do
-!!! end if
-
- !
- ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="red">abla^4 u
- ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- ! applied recursively.
- ! strictly only valid for h_mom_eddy_visc4 == constant
- !
-!!! if (config_h_mom_eddy_visc4 > 0.0) then
-!!! allocate(delsq_divergence(nVertLevels, nCells+1))
-!!! allocate(delsq_u(nVertLevels, nEdges+1))
-!!! allocate(delsq_circulation(nVertLevels, nVertices+1))
-!!! allocate(delsq_vorticity(nVertLevels, nVertices+1))
-
-!!! delsq_u(:,:) = 0.0
-
- ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
-!!! do iEdge=1,grid % nEdges
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! vertex1 = verticesOnEdge(1,iEdge)
-!!! vertex2 = verticesOnEdge(2,iEdge)
-
-!!! do k=1,nVertLevels
-
-!!! delsq_normalVelocity(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-!!! -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
-!!! end do
-!!! end do
-
- ! vorticity using </font>
<font color="red">abla^2 u
-!!! delsq_circulation(:,:) = 0.0
-!! do iEdge=1,nEdges
-!!! vertex1 = verticesOnEdge(1,iEdge)
-!!! vertex2 = verticesOnEdge(2,iEdge)
-!!! do k=1,nVertLevels
-!!! delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
-!!! - dcEdge(iEdge) * delsq_u(k,iEdge)
-!!! delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
-!!! + dcEdge(iEdge) * delsq_u(k,iEdge)
-!!! end do
-!!! end do
-!!! do iVertex=1,nVertices
-!!! r = 1.0 / areaTriangle(iVertex)
-!!! do k=1,nVertLevels
-!!! delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
-!!! end do
-!!! end do
-
- ! Divergence using </font>
<font color="red">abla^2 u
-!!! delsq_divergence(:,:) = 0.0
-!!! do iEdge=1,nEdges
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! do k=1,nVertLevels
-!!! delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
-!!! + delsq_u(k,iEdge)*dvEdge(iEdge)
-!!! delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
-!!! - delsq_u(k,iEdge)*dvEdge(iEdge)
-!!! end do
-!!! end do
-!!! do iCell = 1,nCells
-!!! r = 1.0 / areaCell(iCell)
-!!! do k = 1,nVertLevels
-!!! delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
-!!! end do
-!!! end do
-
- ! Compute - \kappa </font>
<font color="red">abla^4 u
- ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="red">abla^2 u) )
-!!! do iEdge=1,grid % nEdgesSolve
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! vertex1 = verticesOnEdge(1,iEdge)
-!!! vertex2 = verticesOnEdge(2,iEdge)
-
-!!! do k=1,nVertLevels
-
-!!! u_diffusion = ( delsq_divergence(k,cell2) &
-!!! - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
-!!! -( delsq_vorticity(k,vertex2) &
-!!! - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
-
-!!! u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
-!!! tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
-
-!!! end do
-!!! end do
-
-!!! deallocate(delsq_divergence)
-!!! deallocate(delsq_u)
-!!! deallocate(delsq_circulation)
-!!! deallocate(delsq_vorticity)
-
-!!! end if
-
- ! Compute u (velocity) tendency from wind stress (u_src)
-!!! if(config_wind_stress) then
-!!! do iEdge=1,grid % nEdges
-!!! tend_u(1,iEdge) = tend_u(1,iEdge) &
-!!! + u_src(1,iEdge)/rho_ref/layerThicknessEdge(1,iEdge)
-!!! end do
-!!! endif
-
-!!! if (config_bottom_drag) then
-!!! do iEdge=1,grid % nEdges
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
-!!! ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
-!!! + ke(1,cellsOnEdge(2,iEdge)))
-
-!!! tend_u(1,iEdge) = tend_u(1,iEdge) &
-!!! - 1.0e-3*normalVelocity(1,iEdge) &
-!!! *sqrt(2.0*ke_edge)/layerThicknessEdge(1,iEdge)
-!!! end do
-!!! endif
-
- end subroutine land_ice_compute_tend
-
- subroutine land_ice_compute_scalar_tend(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
- real (kind=RKIND) :: flux, tracer_edge, r
- real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
- integer, dimension(:,:), pointer :: boundaryEdge
- real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
-
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
- integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThicknessEdge
-
- normalVelocity => s % normalVelocity % array
- layerThicknessEdge => s % layerThicknessEdge % array
- dcEdge => grid % dcEdge % array
- deriv_two => grid % deriv_two % array
- dvEdge => grid % dvEdge % array
- tracers => s % tracers % array
- cellsOnEdge => grid % cellsOnEdge % array
- boundaryCell=> grid % boundaryCell % array
- boundaryEdge=> grid % boundaryEdge % array
- areaCell => grid % areaCell % array
- tracer_tend => tend % tracers % array
-
- coef_3rd_order = 0.
- if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
- if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
-
- tracer_tend(:,:,:) = 0.0
-
- if (config_tracer_adv_order == 2) then
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * tracer_edge
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- end do
- end do
- end if
- end do
-
- else if (config_tracer_adv_order == 3) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (normalVelocity(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * normalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- else u <= 0:
- else
- flux = dvEdge(iEdge) * normalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
- end do
-
- else if (config_tracer_adv_order == 4) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- flux = dvEdge(iEdge) * normalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
- end do
-
- endif ! if (config_tracer_adv_order == 2 )
-
- !
- ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="red">abla \phi)
- !
-!!! if ( config_h_tracer_eddy_diff2 > 0.0 ) then
-
- !
- ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
- !
-!!! allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
-!!! boundaryMask = 1.0
-!!! where(boundaryEdge.eq.1) boundaryMask=0.0
-
-!!! do iEdge=1,grid % nEdges
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! invAreaCell1 = 1.0/areaCell(cell1)
-!!! invAreaCell2 = 1.0/areaCell(cell2)
-
-!!! do k=1,grid % nVertLevels
-!!! do iTracer=1, grid % nTracers
-!!! ! \kappa_2 </font>
<font color="red">abla \phi on edge
-!!! tracer_turb_flux = config_h_tracer_eddy_diff2 &
-!!! *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
-
- ! div(h \kappa_2 </font>
<font color="red">abla \phi) at cell center
-!!! flux = dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
-!!! tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
-!!! tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
-!!! end do
-!!! end do
-
-!!! end do
-
-!!! deallocate(boundaryMask)
-
-!!! end if
-
- !
- ! tracer tendency: del4 horizontal tracer diffusion, &
- ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="red">abla \phi)])
- !
-!!! if ( config_h_tracer_eddy_diff4 > 0.0 ) then
-
- !
- ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
- !
-!!! allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
-!!! boundaryMask = 1.0
-!!! where(boundaryEdge.eq.1) boundaryMask=0.0
-
-!!! allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
-
-!!! delsq_tracer(:,:,:) = 0.
-
- ! first del2: div(h </font>
<font color="red">abla \phi) at cell center
-!!! do iEdge=1,grid % nEdges
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-
-!!! do k=1,grid % nVertLevels
-!!! do iTracer=1, grid % nTracers
-!!! delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
-!!! + dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
-!!! delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
-!!! - dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
-!!! end do
-!!! end do
-
-!!! end do
-
-!!! do iCell = 1, grid % nCells
-!!! r = 1.0 / grid % areaCell % array(iCell)
-!!! do k=1,grid % nVertLevels
-!!! do iTracer=1,grid % nTracers
-!!! delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
-!!! end do
-!!! end do
-!!! end do
-
- ! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
-!!! do iEdge=1,grid % nEdges
-!!! cell1 = grid % cellsOnEdge % array(1,iEdge)
-!!! cell2 = grid % cellsOnEdge % array(2,iEdge)
-!!! invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
-!!! invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
-
-!!! do k=1,grid % nVertLevels
-!!! do iTracer=1,grid % nTracers
-!!! tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
-!!! flux = dvEdge(iEdge) * tracer_turb_flux
-!!! tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
-!!! tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
-!!! end do
-!!! enddo
-
-!!! end do
-
-!!! deallocate(delsq_tracer)
-!!! deallocate(boundaryMask)
-
-!!! end if
-
- end subroutine land_ice_compute_scalar_tend
-
-
- subroutine land_ice_compute_solve_diagnostics(dt, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: dt
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux !!!, vorticity_abs, workpv
-
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, areaTriangle
-!!! real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge
-!!! real (kind=RKIND), dimension(:,:), pointer :: vh
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, layerThicknessEdge, &
- layerThickness, normalVelocity, tend_layerThickness, &
-!!! real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, layerThicknessEdge, h, u, v, tend_h, tend_u, &
-!!! circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, &
-!!! gradPVn, gradPVt, divergence, vorticity_cell, &
- layerThicknessVertex
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, &
- edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, h1, h2
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
-
- layerThickness => s % layerThickness % array
- normalVelocity => s % normalVelocity % array
-!!! v => s % v % array
-!!! vh => s % vh % array
- layerThicknessEdge => s % layerThicknessEdge % array
- layerThicknessVertex => s % layerThicknessVertex % array
- tend_layerThickness => s % layerThickness % array
-!!! tend_u => s % u % array
-!!! circulation => s % circulation % array
-!!! vorticity => s % vorticity % array
-!!! divergence => s % divergence % array
-!!! ke => s % ke % array
-!!! pv_edge => s % pv_edge % array
-!!! pv_vertex => s % pv_vertex % array
-!!! pv_cell => s % pv_cell % array
-!!! vorticity_cell => s % vorticity_cell % array
-!!! gradPVn => s % gradPVn % array
-!!! gradPVt => s % gradPVt % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
-!!! h_s => grid % h_s % array
-!!! fVertex => grid % fVertex % array
-!!! fEdge => grid % fEdge % array
- deriv_two => grid % deriv_two % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- boundaryEdge => grid % boundaryEdge % array
- boundaryCell => grid % boundaryCell % array
-
- !
- ! Find those cells that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- enddo
- enddo
-
- !
- ! Compute height on cell edges at velocity locations
- ! Namelist options control the order of accuracy of the reconstructed layerThicknessEdge value
- !
-
- coef_3rd_order = 0.
- if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
- if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
- if (config_thickness_adv_order == 2) then
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- layerThicknessEdge(k,iEdge) = 0.5 * (layerThickness(k,cell1) + layerThickness(k,cell2))
- end do
- end if
- end do
-
- else if (config_thickness_adv_order == 3) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * layerThickness(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * layerThickness(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (normalVelocity(k,iEdge) > 0) then
- layerThicknessEdge(k,iEdge) = &
- 0.5*(layerThickness(k,cell1) + layerThickness(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else normalVelocity <= 0:
- else
- layerThicknessEdge(k,iEdge) = &
- 0.5*(layerThickness(k,cell1) + layerThickness(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- else if (config_thickness_adv_order == 4) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * layerThickness(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * layerThickness(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- layerThicknessEdge(k,iEdge) = &
- 0.5*(layerThickness(k,cell1) + layerThickness(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- endif ! if(config_thickness_adv_order == 2)
-
- !
- ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
- ! used to when reading for edges that do not exist
- !
- normalVelocity(:,nEdges+1) = 0.0
-
- !
- ! Compute circulation and relative vorticity at each vertex
- !
-!!! circulation(:,:) = 0.0
-!!! do iEdge=1,nEdges
-!!! do k=1,nVertLevels
-!!! circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * normalVelocity(k,iEdge)
-!!! circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * normalVelocity(k,iEdge)
-!!! end do
-!!! end do
-!!! do iVertex=1,nVertices
-!!! do k=1,nVertLevels
-!!! vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
-!!! end do
-!!! end do
-
-
- !
- ! Compute the divergence at each cell center
- !
-!!! divergence(:,:) = 0.0
-!!! do iEdge=1,nEdges
-!!! cell1 = cellsOnEdge(1,iEdge)
-!!! cell2 = cellsOnEdge(2,iEdge)
-!!! if (cell1 <= nCells) then
-!!! do k=1,nVertLevels
-!!! divergence(k,cell1) = divergence(k,cell1) + normalVelocity(k,iEdge)*dvEdge(iEdge)
-!!! enddo
-!!! endif
-!!! if(cell2 <= nCells) then
-!!! do k=1,nVertLevels
-!!! divergence(k,cell2) = divergence(k,cell2) - normalVelocity(k,iEdge)*dvEdge(iEdge)
-!!! enddo
-!!! end if
-!!! end do
-!!! do iCell = 1,nCells
-!!! r = 1.0 / areaCell(iCell)
-!!! do k = 1,nVertLevels
-!!! divergence(k,iCell) = divergence(k,iCell) * r
-!!! enddo
-!!! enddo
-
- !
- ! Compute kinetic energy in each cell
- !
-!!! ke(:,:) = 0.0
-!!! do iCell=1,nCells
-!!! do i=1,nEdgesOnCell(iCell)
-!!! iEdge = edgesOnCell(i,iCell)
-!!! do k=1,nVertLevels
-!!! ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * normalVelocity(k,iEdge)**2.0
-!!! end do
-!!! end do
-!!! do k=1,nVertLevels
-!!! ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
-!!! end do
-!!! end do
-
- !
- ! Compute v (tangential) velocities
- !
-!!! v(:,:) = 0.0
-!!! do iEdge = 1,nEdges
-!!! do i=1,nEdgesOnEdge(iEdge)
-!!! eoe = edgesOnEdge(i,iEdge)
-!!! do k = 1,nVertLevels
-!!! v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * normalVelocity(k, eoe)
-!!! end do
-!!! end do
-!!! end do
-
-#ifdef NCAR_FORMULATION
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
-!!! vh(:,:) = 0.0
-!!! do iEdge=1,grid % nEdgesSolve
-!!! do j=1,nEdgesOnEdge(iEdge)
-!!! eoe = edgesOnEdge(j,iEdge)
-!!! do k=1,nVertLevels
-!!! vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * normalVelocity(k,eoe) * layerThicknessEdge(k,eoe)
-!!! end do
-!!! end do
-!!! end do
-#endif
-
-
- !
- ! Compute height at vertices, pv at vertices, and average pv to edge locations
- ! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
- !
- do iVertex = 1,nVertices
- do k=1,nVertLevels
- layerThicknessVertex(k,iVertex) = 0.0
- do i=1,grid % vertexDegree
- layerThicknessVertex(k,iVertex) = layerThicknessVertex(k,iVertex) + layerThickness(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- layerThicknessVertex(k,iVertex) = layerThicknessVertex(k,iVertex) / areaTriangle(iVertex)
-
-!!! pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / layerThicknessVertex(k,iVertex)
- end do
- end do
-
-
- !
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
- !
-!!! do iEdge = 1,nEdges
-!!! do k = 1,nVertLevels
-!!! gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
-!!! dvEdge(iEdge)
-!!! enddo
-!!! enddo
-
- !
- ! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells )
- !
-!!! pv_edge(:,:) = 0.0
-!!! do iVertex = 1,nVertices
-!!! do i=1,grid % vertexDegree
-!!! iEdge = edgesOnVertex(i,iVertex)
-!!! do k=1,nVertLevels
-!!! pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
-!!! end do
-!!! end do
-!!! end do
-
-
- !
- ! Modify PV edge with upstream bias.
- !
-!!! do iEdge = 1,nEdges
-!!! do k = 1,nVertLevels
-!!! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
-!!! enddo
-!!! enddo
-
-
- !
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
-!!! pv_cell(:,:) = 0.0
-!!! vorticity_cell(:,:) = 0.0
-!!! do iVertex = 1, nVertices
-!!! do i=1,grid % vertexDegree
-!!! iCell = cellsOnVertex(i,iVertex)
-!!! if (iCell <= nCells) then
-!!! do k = 1,nVertLevels
-!!! pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
-!!! vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
-!!! enddo
-!!! endif
-!!! enddo
-!!! enddo
-
-
- !
- ! Compute gradient of PV in normal direction
- ! ( this computes gradPVn for all edges bounding real cells )
- !
-!!! gradPVn(:,:) = 0.0
-!!! do iEdge = 1,nEdges
-!!! if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
-!!! do k = 1,nVertLevels
-!!! gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
-!!! dcEdge(iEdge)
-!!! enddo
-!!! endif
-!!! enddo
-
- ! Modify PV edge with upstream bias.
- !
-!!! do iEdge = 1,nEdges
-!!! do k = 1,nVertLevels
-!!! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * normalVelocity(k,iEdge) * dt * gradPVn(k,iEdge)
-!!! enddo
-!!! enddo
-
- !
- ! set pv_edge = fEdge / layerThicknessEdge at boundary points
- !
- ! if (maxval(boundaryEdge).ge.0) then
- ! do iEdge = 1,nEdges
- ! cell1 = cellsOnEdge(1,iEdge)
- ! cell2 = cellsOnEdge(2,iEdge)
- ! do k = 1,nVertLevels
- ! if(boundaryEdge(k,iEdge).eq.1) then
- ! v(k,iEdge) = 0.0
- ! if(cell1.gt.0) then
- ! h1 = layerThickness(k,cell1)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
- ! layerThicknessEdge(k,iEdge) = h1
- ! else
- ! h2 = layerThickness(k,cell2)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
- ! layerThicknessEdge(k,iEdge) = h2
- ! endif
- ! endif
- ! enddo
- ! enddo
- ! endif
-
-
- end subroutine land_ice_compute_solve_diagnostics
-
-
- subroutine land_ice_enforce_boundary_edge(tend, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Enforce any boundary conditions on the normal velocity at each edge
- !
- ! Input: grid - grid metadata
- !
- ! Output: tend_u set to zero at boundaryEdge == 1 locations
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (mesh_type), intent(in) :: grid
-
- integer, dimension(:,:), pointer :: boundaryEdge
- real (kind=RKIND), dimension(:,:), pointer :: tend_u
- integer :: nCells, nEdges, nVertices, nVertLevels
- integer :: iEdge, k
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- boundaryEdge => grid % boundaryEdge % array
-!!! tend_u => tend % u % array
-
- if(maxval(boundaryEdge).le.0) return
-
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
-
- if(boundaryEdge(k,iEdge).eq.1) then
-!!! tend_u(k,iEdge) = 0.0
- endif
-
- enddo
- enddo
-
- end subroutine land_ice_enforce_boundary_edge
-
-
end module land_ice_time_integration
Modified: branches/land_ice_projects/test_cases/dome/visualize_dome.py
===================================================================
--- branches/land_ice_projects/test_cases/dome/visualize_dome.py 2012-03-07 22:50:57 UTC (rev 1603)
+++ branches/land_ice_projects/test_cases/dome/visualize_dome.py 2012-03-07 23:27:50 UTC (rev 1604)
@@ -55,6 +55,8 @@
xCell = f.variables['xCell']
yCell = f.variables['yCell']
temperature = f.variables['temperature']
+lowerSurface = f.variables['lowerSurface']
+upperSurface = f.variables['upperSurface']
vert_levs = f.dimensions['nVertLevels']
@@ -92,6 +94,19 @@
plt.draw()
fig = plt.figure(2)
+ax = fig.add_subplot(121, aspect='equal')
+plt.scatter(xCell[:], yCell[:], 80, lowerSurface[time_slice,:], marker='h', edgecolors='none')
+plt.colorbar()
+plt.title('lower surface at time ' + str(time_slice) )
+plt.draw()
+ax = fig.add_subplot(122, aspect='equal')
+plt.scatter(xCell[:], yCell[:], 80, upperSurface[time_slice,:], marker='h', edgecolors='none')
+plt.colorbar()
+plt.title('upper surface at time ' + str(time_slice) )
+plt.draw()
+
+
+fig = plt.figure(3)
for templevel in xrange(0,vert_levs+2):
ax = fig.add_subplot(3,4,templevel+1, aspect='equal')
var_slice = temperature[time_slice,:,templevel]
</font>
</pre>