<p><b>mpetersen@lanl.gov</b> 2012-10-11 16:58:41 -0600 (Thu, 11 Oct 2012)</p><p>Revise basin routine in basin_pbc_mrp branch to create and write the bottomDepth and refBottomDepth variables needed for partial bottom cells.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin_pbc_mrp/src/basin.F
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/src/basin.F 2012-10-11 22:02:45 UTC (rev 2207)
+++ branches/ocean_projects/basin_pbc_mrp/src/basin.F 2012-10-11 22:58:41 UTC (rev 2208)
@@ -48,7 +48,7 @@
real, allocatable, dimension(:) :: areaCell, areaTriangle, dcEdge, dvEdge, angleEdge
real, allocatable, dimension(:,:) :: kiteAreasOnVertex, weightsOnEdge
-real, allocatable, dimension(:) :: fEdge, fVertex, h_s, work1
+real, allocatable, dimension(:) :: fEdge, fVertex, bottomDepth, work1
real, allocatable, dimension(:,:) :: u_src
real, allocatable, dimension(:,:,:) :: u, v, h
real, allocatable, dimension(:,:,:) :: rho
@@ -90,8 +90,14 @@
logical, parameter :: real_bathymetry=.true.
logical, parameter :: eliminate_inland_seas=.true.
logical, parameter :: l_woce = .true.
+logical, parameter :: write_OpenDX_flag = .true.
+! Set the number of top layers that are not allowed to have land, usually three.
+integer, parameter :: top_layers_without_land = 3
+
+
+
! Step 3: Specify some Parameters
real (kind=8), parameter :: &
h_total_max = 2000.0, &
@@ -107,6 +113,7 @@
! new grid variables
+real, dimension(nVertLevelsMOD) :: hZLevel, refBottomDepth
integer :: nCellsNew, nEdgesNew, nVerticesNew
integer :: maxEdgesNew, maxEdges2New, TWONew, vertexDegreeNew, nVertLevelsNew
integer, allocatable, dimension(:) :: indexToCellIDNew, indexToEdgeIDNew, indexToVertexIDNew
@@ -121,7 +128,7 @@
real, allocatable, dimension(:) :: areaCellNew, areaTriangleNew, dcEdgeNew, dvEdgeNew, angleEdgeNew
real, allocatable, dimension(:,:) :: kiteAreasOnVertexNew, weightsOnEdgeNew, normalsNew
-real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, h_sNew, hZLevel, referenceBottomDepth
+real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, bottomDepthNew
real, allocatable, dimension(:,:) :: u_srcNew
real, allocatable, dimension(:,:) :: windStressMonthlyNew
real, allocatable, dimension(:,:,:) :: uNew, vNew, hNew
@@ -132,7 +139,6 @@
! mapping variables
integer, allocatable, dimension(:) :: kmt, maxLevelCellNew
integer, allocatable, dimension(:) :: cellMap, edgeMap, vertexMap
-real, allocatable, dimension(:) :: depthCell
! work variables
integer :: i,j,jNew,k,jEdge,jEdgeNew,iVertex1New,iVertex2New,iCell1New,iCell2New
@@ -151,6 +157,7 @@
write(6,*)
! get depth profile for later
+write(6,*) ' calling get_dz'
call get_dz
! get grid
@@ -340,9 +347,10 @@
call write_graph
! write OpenDx file
-write(6,*) ' calling write_OpenDX'
-write(6,*)
-call write_OpenDX( on_a_sphere, &
+if (write_OpenDX_flag) then
+ write(6,*) ' calling write_OpenDX'
+ write(6,*)
+ call write_OpenDX( on_a_sphere, &
nCellsNew, &
nVerticesNew, &
nEdgesNew, &
@@ -365,11 +373,11 @@
areaCellNew, &
maxLevelCellNew, &
meshDensityNew, &
- depthCell, &
+ bottomDepthNew, &
temperatureNew(1,1,:), &
kiteAreasOnVertexNew )
+endif
-
!do iCell=1,nCellsNew
!ulon = 1.0; ulat = 0.0
!xin = xCellNew(iCell); yin = yCellNew(iCell); zin = zCellNew(iCell)
@@ -436,10 +444,12 @@
fEdgeNew(:) = 0.0
fVertexNew(:) = 0.0
- h_sNew(:) = 0.0
+ bottomDepthNew(:) = 0.0
uNew(:,:,:) = 0.0
vNew(:,:,:) = 0.0
+!!!!!!!!!!!! mrp need to define hZLevel in get_dz, below.
+
! basin-mod
! setting for three levels - Set h values for isopycnal system
write(6,*) ' setting three levels for isopycnal system'
@@ -447,16 +457,18 @@
hNew(1,1,:) = 500.0
hNew(1,2,:) = 1250.0
hNew(1,3,:) = 3250.0
- h_sNew(:) = -( hNew(1,1,:) + hNew(1,2,:) + hNew(1,3,:) )
+ bottomDepthNew(:) = -( hNew(1,1,:) + hNew(1,2,:) + hNew(1,3,:) )
else
hNew(1,1,:) = 3250.0
- h_sNew(:) = -( hNew(1,1,:) )
+ bottomDepthNew(:) = -( hNew(1,1,:) )
endif
hZLevel = hNew(1,:,1)
- referenceBottomDepth(1) = hZLevel(1)
+
+!!!!!!!!!!! remove these lines mrp
+ refBottomDepth(1) = hZLevel(1)
do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
+ refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
end do
! basin-mod
@@ -699,7 +711,6 @@
iNoData = 0
do iCell=1,nCellsNew
hNew(1,:,iCell) = dz(:)
- hZLevel = dz
if(mod(iCell,100).eq.0) write(6,*) 'l_woce t and s',iCell
rlon = lonCellNew(iCell)/dtr
rlat = latCellNew(iCell)/dtr
@@ -1086,7 +1097,7 @@
allocate(fEdge(nEdges))
allocate(fVertex(nVertices))
-allocate(h_s(nCells))
+allocate(bottomDepth(nCells))
allocate(work1(nCells))
allocate(u_src(nVertLevels,nEdges))
allocate(u(1,nVertLevels,nEdges))
@@ -1105,7 +1116,7 @@
cellsOnCell=0; verticesOnCell=0; verticesOnEdge=0
edgesOnVertex=0; cellsOnVertex=0; kiteAreasOnVertex=0
-fEdge=0; fVertex=0; h_s=0; u_src=0; work1=0
+fEdge=0; fVertex=0; bottomDepth=0; u_src=0; work1=0
u=0; v=0; h=0; rho=0
@@ -1149,7 +1160,7 @@
kiteAreasOnVertex, &
fEdge, &
fVertex, &
- h_s, &
+ bottomDepth, &
u, &
v, &
h &
@@ -1194,7 +1205,7 @@
write(6,*) ' kiteAreasOnVertex : ', minval(kiteAreasOnVertex), maxval(kiteAreasOnVertex)
write(6,*) ' fEdge : ', minval(fEdge), maxval(fEdge)
write(6,*) ' fVertex : ', minval(fVertex), maxval(fVertex)
-write(6,*) ' h_s : ', minval(h_s), maxval(h_s)
+write(6,*) ' bottomDepth : ', minval(bottomDepth), maxval(bottomDepth)
write(6,*) ' u : ', minval(u), maxval(u)
write(6,*) ' v : ', minval(v), maxval(v)
write(6,*) ' h : ', minval(h), maxval(h)
@@ -1275,7 +1286,7 @@
kiteAreasOnVertexNew, &
fEdgeNew, &
fVertexNew, &
- h_sNew, &
+ bottomDepthNew, &
boundaryEdgeNew, &
boundaryVertexNew, &
u_srcNew, &
@@ -1292,7 +1303,7 @@
temperatureRestoreMonthlyNew, &
salinityRestoreMonthlyNew, &
hZLevel, &
- referenceBottomDepth &
+ refBottomDepth &
)
call write_netcdf_finalize
@@ -1391,8 +1402,14 @@
write(6,*) minval(ztopo), maxval(ztopo)
do iCell=1,nCells
+
+ ! Convert from radians to degrees
rlon = lonCell(iCell) / dtr
rlat = latCell(iCell) / dtr
+
+ ! Find nearest coordinate in topo file.
+ ! This is 1/30th degree topo data, so multiply degrees by 30 and
+ ! round to get index.
ix = nint((rlon+180)*30) + nx
ix = mod(ix,nx)+1
iy = nint((rlat+90 )*30)
@@ -1401,8 +1418,10 @@
zdata = ztopo(ix,iy)
+ ! zdata is less than zero for ocean points.
if(zdata.lt.0.0) then
zdata = -zdata
+ bottomDepth(iCell) = zdata
r = 0
kmt_flag=.false.
do k=1,nVertLevelsMod
@@ -1414,8 +1433,15 @@
endif
endif
enddo
- if(kmt(iCell).eq.0) kmt(iCell)=nVertLevelsMod
+
+ ! zdata is deeper than deepest cell
+ if (kmt(iCell).eq.0) then
+ kmt(iCell)=nVertLevelsMod
+ bottomDepth(iCell) = refBottomDepth(nVertLevelsMod)
+ endif
+
! write(6,*) kmt(iCell)
+
endif
! if(zdata.lt.0.0) kmt(iCell) = nVertLevelsMod
@@ -1435,7 +1461,10 @@
iCell1 = cellsOnCell(j,iCell)
kmt_neighbor_max = max(kmt_neighbor_max,kmt(iCell1))
enddo
- kmt(iCell) = min(kmt(iCell),kmt_neighbor_max)
+ if (kmt(iCell).gt.kmt_neighbor_max) then
+ kmt(iCell) = kmt_neighbor_max
+ bottomDepth(iCell) = refBottomDepth(kmt(iCell))
+ endif
enddo
if(eliminate_inland_seas) then
@@ -1445,10 +1474,12 @@
KMT)
endif
-if(real_bathymetry) then
- where(kmt.eq.1) kmt=3
- where(kmt.eq.2) kmt=3
-endif
+! do not allow land or PBCs in top layers
+k = top_layers_without_land
+where(kmt.gt.0.and.kmt.le.k)
+ bottomDepth = refBottomDepth(k)
+ kmt=k
+endwhere
end subroutine define_kmt
@@ -1561,18 +1592,15 @@
allocate(areaCellNew(nCellsNew))
allocate(areaTriangleNew(nVerticesNew))
allocate(maxLevelCellNew(nCellsNew))
-allocate(depthCell(nCellsNew))
allocate(fEdgeNew(nEdgesNew))
allocate(fVertexNew(nVerticesNew))
-allocate(h_sNew(nCellsNew))
+allocate(bottomDepthNew(nCellsNew))
allocate(u_srcNew(nVertLevelsNew,nEdgesNew))
allocate(windStressMonthlyNew(12,nEdgesNew))
allocate(uNew(1,nVertLevelsNew,nEdgesNew))
allocate(vNew(1,nVertLevelsNew,nEdgesNew))
allocate(hNew(1,nVertLevelsNew,nCellsNew))
-allocate(hZLevel(nVertLevelsNew))
-allocate(referenceBottomDepth(nVertLevelsNew))
allocate(rhoNew(1,nVertLevelsNew,nCellsNew))
allocate(temperatureNew(1,nVertLevelsNew,nCellsNew))
allocate(salinityNew(1,nVertLevelsNew,nCellsNew))
@@ -1589,7 +1617,7 @@
xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
-fEdgeNew=0; fVertexNew=0; h_sNew=0; u_srcNew=0; windStressMonthlyNew = 0
+fEdgeNew=0; fVertexNew=0; bottomDepthNew=0; u_srcNew=0; windStressMonthlyNew = 0
uNew=0; vNew=0; hNew=0; rhoNew=0
temperatureNew=0; salinityNew=0; tracer1New=0;
@@ -1611,7 +1639,7 @@
meshDensityNew(jNew)=meshDensity(i)
areaCellNew(jNew)=areaCell(i)
maxLevelCellNew(jNew) = kmt(i)
- depthCell(jNew) = kmt(i)
+ bottomDepthNew(jNew) = bottomDepth(i)
endif
enddo
@@ -1878,6 +1906,17 @@
subroutine get_dz
integer k
+if(.not.real_bathymetry) then
+
+!!!!!!!!!! mrp add cases here, and a step to change this.
+
+ write(6,*) ' equally spaced layers'
+ do i = 1, nVertLevelsMOD
+ hZLevel(i) = h_total_max / nVertLevelsMOD
+ end do
+
+else
+
dz( 1) = 1001.244 ! 5.006218 10.01244
dz( 2) = 1011.258 ! 15.06873 20.12502
dz( 3) = 1031.682 ! 25.28342 30.44183
@@ -1921,11 +1960,21 @@
dz = dz / 100.0
- write(6,*)
- do k=1,40
- write(6,*) k,dz(k)
+ hZLevel = dz
+
+endif
+
+ refBottomDepth(1) = hZLevel(1)
+ do k = 2,nVertLevelsMod
+ refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
+ end do
+
+ write(6,*) ' k hZLevel refBottomDepth'
+ do k=1,nVertLevelsMod
+ write(6,'(i5,2f10.2)') k,hZLevel(k), refBottomDepth(k)
enddo
write(6,*)
end subroutine get_dz
+
end program map_to_basin
Modified: branches/ocean_projects/basin_pbc_mrp/src/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/basin_pbc_mrp/src/module_write_netcdf.F 2012-10-11 22:02:45 UTC (rev 2207)
+++ branches/ocean_projects/basin_pbc_mrp/src/module_write_netcdf.F 2012-10-11 22:58:41 UTC (rev 2208)
@@ -50,7 +50,7 @@
integer :: wrVarIDkiteAreasOnVertex
integer :: wrVarIDfEdge
integer :: wrVarIDfVertex
- integer :: wrVarIDh_s
+ integer :: wrVarIDbottomDepth
integer :: wrVarIDu
integer :: wrVarIDboundaryEdge
integer :: wrVarIDboundaryVertex
@@ -67,7 +67,7 @@
integer :: wrVarIDtemperatureRestoreMonthly
integer :: wrVarIDsalinityRestoreMonthly
integer :: wrVarIDhZLevel
- integer :: wrVarIDreferenceBottomDepth
+ integer :: wrVarIDrefBottomDepth
integer :: wrLocalnCells
integer :: wrLocalnEdges
@@ -226,7 +226,7 @@
dimlist( 1) = wrDimIDnVertices
nferr = nf_def_var(wr_ncid, 'fVertex', NF_DOUBLE, 1, dimlist, wrVarIDfVertex)
dimlist( 1) = wrDimIDnCells
- nferr = nf_def_var(wr_ncid, 'h_s', NF_DOUBLE, 1, dimlist, wrVarIDh_s)
+ nferr = nf_def_var(wr_ncid, 'bottomDepth', NF_DOUBLE, 1, dimlist, wrVarIDbottomDepth)
dimlist( 1) = wrDimIDnCells
nferr = nf_def_var(wr_ncid, 'temperatureRestore', NF_DOUBLE, 1, dimlist, wrVarIDtemperatureRestore)
dimlist( 1) = wrDimIDnCells
@@ -242,7 +242,7 @@
dimlist( 1) = wrDimIDnVertLevels
nferr = nf_def_var(wr_ncid, 'hZLevel', NF_DOUBLE, 1, dimlist, wrVarIDhZLevel)
dimlist( 1) = wrDimIDnVertLevels
- nferr = nf_def_var(wr_ncid, 'referenceBottomDepth', NF_DOUBLE, 1, dimlist, wrVarIDreferenceBottomDepth)
+ nferr = nf_def_var(wr_ncid, 'refBottomDepth', NF_DOUBLE, 1, dimlist, wrVarIDrefBottomDepth)
dimlist( 1) = wrDimIDnVertLevels
dimlist( 2) = wrDimIDnEdges
dimlist( 3) = wrDimIDTime
@@ -337,7 +337,7 @@
kiteAreasOnVertex, &
fEdge, &
fVertex, &
- h_s, &
+ bottomDepth, &
boundaryEdge, &
boundaryVertex, &
u_src, &
@@ -354,7 +354,7 @@
temperatureRestoreMonthly, &
salinityRestoreMonthly, &
hZLevel, &
- referenceBottomDepth &
+ refBottomDepth &
)
implicit none
@@ -401,7 +401,7 @@
real (kind=8), dimension(:,:), intent(in) :: kiteAreasOnVertex
real (kind=8), dimension(:), intent(in) :: fEdge
real (kind=8), dimension(:), intent(in) :: fVertex
- real (kind=8), dimension(:), intent(in) :: h_s
+ real (kind=8), dimension(:), intent(in) :: bottomDepth
integer, dimension(:,:), intent(in) :: boundaryEdge
integer, dimension(:,:), intent(in) :: boundaryVertex
real (kind=8), dimension(:,:), intent(in) :: u_src
@@ -418,7 +418,7 @@
real (kind=8), dimension(:,:), intent(in) :: temperatureRestoreMonthly
real (kind=8), dimension(:,:), intent(in) :: salinityRestoreMonthly
real (kind=8), dimension(:), intent(in) :: hZLevel
- real (kind=8), dimension(:), intent(in) :: referenceBottomDepth
+ real (kind=8), dimension(:), intent(in) :: refBottomDepth
integer :: nferr
@@ -609,7 +609,7 @@
start1(1) = 1
count1( 1) = wrLocalnCells
- nferr = nf_put_vara_double(wr_ncid, wrVarIDh_s, start1, count1, h_s)
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDbottomDepth, start1, count1, bottomDepth)
start2(2) = 1
count2( 1) = wrLocalnVertLevels
@@ -662,7 +662,7 @@
start1(1) = 1
count1( 1) = wrLocalnVertLevels
- nferr = nf_put_vara_double(wr_ncid, wrVarIDreferenceBottomDepth, start1, count1, referenceBottomDepth)
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDrefBottomDepth, start1, count1, refBottomDepth)
start3(3) = time
count3( 1) = wrLocalnVertLevels
</font>
</pre>