<p><b>mpetersen@lanl.gov</b> 2010-11-09 10:55:35 -0700 (Tue, 09 Nov 2010)</p><p>BRANCH COMMIT: Add equation of state flag to ocean core.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/topography2_mrp/namelist.input.ocean 2010-11-08 19:06:52 UTC (rev 602)
+++ branches/ocean_projects/topography2_mrp/namelist.input.ocean 2010-11-09 17:55:35 UTC (rev 603)
@@ -41,6 +41,9 @@
config_vert_viscosity = 1.0e-4
config_vert_diffusion = 1.0e-4
/
+&eos
+ config_eos_type = 'linear'
+/
&advection
config_vert_tracer_adv = 'upwind'
config_tracer_adv_order = 2
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-08 19:06:52 UTC (rev 602)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-09 17:55:35 UTC (rev 603)
@@ -29,6 +29,7 @@
namelist real vmix config_vmixTanhDiffMin 1.0e-5
namelist real vmix config_vmixTanhZMid -100
namelist real vmix config_vmixTanhZWidth 100
+namelist character eos config_eos_type linear
namelist character advection config_vert_tracer_adv 'centered'
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F 2010-11-08 19:06:52 UTC (rev 602)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F 2010-11-09 17:55:35 UTC (rev 603)
@@ -25,18 +25,6 @@
type (block_type), pointer :: block_ptr
type (dm_info) :: dminfo
- integer :: iTracer
- real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
- real (kind=RKIND) :: centerx, centery
- integer :: nCells, nEdges, nVertices, nVertLevels
-
- integer, dimension(:), pointer :: &
- maxLevelCell, maxLevelEdge, maxLevelVertex
-
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -98,125 +86,6 @@
block_ptr => block_ptr % next
end do
- ! Initialize z-level grid variables from h, read in from input file.
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- h => block_ptr % state % time_levs(1) % state % h % array
- u => block_ptr % state % time_levs(1) % state % u % array
- rho => block_ptr % state % time_levs(1) % state % rho % array
- tracers => block_ptr % state % time_levs(1) % state % tracers % array
- u_src => block_ptr % mesh % u_src % array
- xCell => block_ptr % mesh % xCell % array
- yCell => block_ptr % mesh % yCell % array
- latCell => block_ptr % mesh % latCell % array
- lonCell => block_ptr % mesh % lonCell % array
-
- hZLevel => block_ptr % mesh % hZLevel % array
- zMidZLevel => block_ptr % mesh % zMidZLevel % array
- zTopZLevel => block_ptr % mesh % zTopZLevel % array
-
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
-
-!mrp temp, to agree with previous test
-!h(1,:) = 500.0
-!h(2,:) = 1250.0
-!h(3,:) = 3250.0
-
- pi=3.1415
- ! Tracer blob in Central Pacific, away from boundaries:
- !latCenter=pi/16; lonCenter=9./8.*pi
-
- ! Tracer blob in Central Pacific, near boundaries:
- latCenter=pi*2./16; lonCenter=13./16.*pi
-
- if (config_vert_grid_type.eq.'zlevel') then
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- hZLevel = h(:,1)
-
- ! hZLevel should be in the grid.nc and restart.nc file
- ! the following lines are in mpas_interface and should be
- ! removed from here. mrp 101019
- zTopZLevel(1) = 0.0
- do iLevel = 1,nVertLevels
- zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
- zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
- enddo
-
- endif
-
- if (config_vert_grid_type.eq.'zlevel') then
- ! Set tracers, if not done in grid.nc file
- !tracers = 0.0
- centerx = 1.0e6
- centery = 1.0e6
- dist = 2.0e5
- do iCell = 1,nCells
- dist = sqrt( (latCell(iCell)-latCenter)**2 + (lonCell(iCell)-lonCenter)**2)
- do iLevel = 1,nVertLevels
- ! for 20 layer test
- ! tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 5.0 ! temperature
- ! tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-
- ! for x3, 25 layer test
- !tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 10.0 ! temperature
- !tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = &
- ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
-
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = block_ptr % mesh % lonCell % array(iCell)
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = block_ptr % mesh % latCell % array(iCell)
-
- ! place tracer blob with radius dist at (centerx,centery)
- !if ( sqrt( (xCell(iCell)-centerx)**2 &
- ! + (yCell(iCell)-centery)**2) &
- ! .lt.dist) then
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- !else
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 0.0
- !endif
-
- rho(iLevel,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) &
- + 7.6e-4*tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell))
-
- enddo
- enddo
-
- endif
-
- ! print some diagnostics
- print '(10a)', 'ilevel',&
- ' rho ',&
- ' min u max u ',&
- ' min u_src max u_src ', &
- ' min h max h ',&
- ' hZLevel zMidZlevel', &
- ' zTopZlevel'
- do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
- minval(u(iLevel,1:nEdges)), maxval(u(iLevel,1:nEdges)), &
- minval(u_src(iLevel,1:nEdges)), maxval(u_src(iLevel,1:nEdges)), &
- minval(h(iLevel,1:nCells)), maxval(h(iLevel,1:nCells)), &
- hZLevel(iLevel),zMidZlevel(iLevel),zTopZlevel(iLevel)
- enddo
-
- print '(10a)', 'itracer ilevel min tracer max tracer'
- do iTracer=1,block_ptr % state % time_levs(1) % state % num_tracers
- do iLevel = 1,nVertLevels
- print '(2i5,20es12.4)', iTracer,ilevel, &
- minval(tracers(itracer,iLevel,1:nCells)), maxval(tracers(itracer,iLevel,1:nCells))
- enddo
- enddo
-
- block_ptr => block_ptr % next
- end do
-
end subroutine setup_sw_test_case
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-11-08 19:06:52 UTC (rev 602)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-11-09 17:55:35 UTC (rev 603)
@@ -1477,22 +1477,9 @@
!
! equation of state
!
- ! For a isopycnal model, density should remain constant.
+ ! For an isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
-
-! Jacket and McDougal: need to add flag and organize in subroutine.
-! dont forget that eos is called after initial conditions in test_cases right now.
- !call equation_of_state(s,grid)
-
-! linear is current default
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- ! Linear equation of state, for the time being
- rho(k,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(s % index_temperature,k,iCell) &
- + 7.6e-4*tracers(s % index_salinity,k,iCell))
- end do
- end do
+ call equation_of_state(s,grid)
endif
!
@@ -1634,15 +1621,73 @@
end subroutine enforce_boundaryEdge
+
subroutine equation_of_state(s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
!
! Input: grid - grid metadata
+ ! s - state: tracers
!
- ! Output: s - computed diagnostics
+ ! Output: s - state: computed density
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: nCells, iCell, k
+ type (dm_info) :: dminfo
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
+
+ if (config_eos_type.eq.'linear') then
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(s % index_temperature,k,iCell) &
+ + 7.6e-4*tracers(s % index_salinity,k,iCell))
+ end do
+ end do
+
+ elseif (config_eos_type.eq.'jm') then
+
+ call equation_of_state_jm(s, grid)
+
+ else
+ print *, ' Incorrect choice of config_eos_type:',&
+ config_eos_type
+ call dmpar_abort(dminfo)
+ endif
+
+ end subroutine equation_of_state
+
+
+ subroutine equation_of_state_jm(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! The UNESCO equation of state computed using the
+ ! potential-temperature-based bulk modulus from Jackett and
+ ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
implicit none
type (state_type), intent(in) :: s
@@ -1830,6 +1875,7 @@
end do
deallocate(pRefEOS)
- end subroutine equation_of_state
+ end subroutine equation_of_state_jm
+
end module time_integration
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-08 19:06:52 UTC (rev 602)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-09 17:55:35 UTC (rev 603)
@@ -53,11 +53,11 @@
call dmpar_abort(dminfo)
endif
-
! Initialize z-level grid variables from h, read in from input file.
block => domain % blocklist
do while (associated(block))
+ h => block % state % time_levs(1) % state % h % array
hZLevel => block % mesh % hZLevel % array
zMidZLevel => block % mesh % zMidZLevel % array
zTopZLevel => block % mesh % zTopZLevel % array
@@ -79,6 +79,11 @@
if (config_vert_grid_type.eq.'zlevel') then
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ ! Upon restart, the correct hZLevel should be in restart.nc
+ if (.not. config_do_restart) hZLevel = h(:,1)
+
! hZLevel should be in the grid.nc and restart.nc file,
! and h for k=1 must be specified there as well.
@@ -97,16 +102,6 @@
endif
maxLevelCell(nCells+1) = 0
- ! mrp insert topography mesa for testing
- ! do iCell = 1,nCells
- ! if (latCell(iCell)> 0.0/180.*3.14.and.&
- ! latCell(iCell)< 30.0/180.*3.14.and. &
- ! lonCell(iCell)>180.0/180.*3.14.and. &
- ! lonCell(iCell)<210.0/180.*3.14) then
- ! maxLevelCell(iCell) = 10
- ! endif
- ! end do
-
! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -253,20 +248,20 @@
type (domain_type), intent(inout) :: domain
integer, intent(in) :: itimestep
real (kind=RKIND), intent(in) :: dt
- type (block_type), pointer :: block_ptr
+ type (block_type), pointer :: block
call timestep(domain, dt)
if (mod(itimestep, config_stats_interval) == 0) then
- block_ptr => domain % blocklist
- if(associated(block_ptr % next)) then
+ block => domain % blocklist
+ if(associated(block % next)) then
write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
'that there is only one block per processor.'
end if
call timer_start("global diagnostics")
call computeGlobalDiagnostics(domain % dminfo, &
- block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
+ block % state % time_levs(2) % state, block % mesh, &
itimestep, dt)
call timer_stop("global diagnostics")
end if
</font>
</pre>