<p><b>sprice@lanl.gov</b> 2011-07-14 17:04:48 -0600 (Thu, 14 Jul 2011)</p><p>branch commit (land ice): More work on ice sheet mesh creation using altered version of basin.F. Still not working.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice/icesheet/src/icesheet.F
===================================================================
--- branches/land_ice/icesheet/src/icesheet.F 2011-07-14 15:23:34 UTC (rev 921)
+++ branches/land_ice/icesheet/src/icesheet.F 2011-07-14 23:04:48 UTC (rev 922)
@@ -1072,7 +1072,7 @@
implicit none
real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
!real (kind=4), allocatable, dimension(:,:) :: ztopo
-real (kind=4), allocatable, dimension(:,:) :: thk
+real (kind=4), allocatable, dimension(:,:) :: thck
integer :: nx, ny, inx, iny, ix, iy
!real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax, xmin, xmax
real :: pi, dtr, thkdata, rlon, rlat, r, ymin, ymax, xmin, xmax
@@ -1124,8 +1124,10 @@
endif
if(real_bathymetry) then
- nx = 301
- ny = 561
+! nx = 301 ! 5 km
+! ny = 561
+ nx = 76 ! 20 km
+ ny = 141
allocate(x(nx))
allocate(y(ny))
! allocate(ztopo(nx,ny))
@@ -1135,10 +1137,10 @@
! allocate(lon(nx,ny))
! allocate(lat(nx,ny))
- allocate(thk(nx,ny))
+ allocate(thck(nx,ny))
! lon = 0.0
! lat = 0.0
- thk = 0.0
+ thck = 0.0
call read_topo_init( inx, iny)
@@ -1146,13 +1148,16 @@
if(inx.ne.nx) stop ' nx topo'
if(iny.ne.ny) stop ' ny topo'
- call read_topo_fields(x,y,thk)
+ call read_topo_fields(x,y,thck)
call read_topo_finalize()
+ print *, 'min(x,y,thck), max(x,y,thck), val(x,y,thck) at 1,1'
write(6,*) minval(x), maxval(x), x(1)
write(6,*) minval(y), maxval(y), y(1)
- write(6,*) minval(thk), maxval(thk), thk(1,1)
+ write(6,*) minval(thck), maxval(thck), thck(1,1)
+ pause
+
do iCell=1,nCells
! ! calc. relative lat/lon
@@ -1169,39 +1174,41 @@
! ix = max(1,ix); ix = min(nx,ix)
! iy = max(1,iy); iy = min(ny,iy)
- ix = xCell(iCell) / 5.0d3
- ix = int( ix )
- iy = yCell(iCell) / 5.0d3
- iy = int( iy )
+ ix = nint( xCell(iCell) / 20.0d3 )
+ iy = nint( yCell(iCell) / 20.0d3 )
print *, 'xCell(iCell)', xCell(iCell)
- print *, 'ix', ix
print *, 'yCell(iCell)', yCell(iCell)
- print *, 'iy', iy
+ print *, 'ix,iy = ', ix,iy
- thkdata = thk(ix,iy)
+ !do i=1:nx; do j=1:ny
+ thkdata = thck(ix,iy)
+ print *, 'x1(ix,iy) = ', x(ix)
+ print *, 'y1(ix,iy) = ', y(iy)
+ print *, 'thck(ix,iy) = ', thkdata
+
if(thkdata .gt. 0)then
kmt(iCell) = nVertLevelsMod
-! elseif( thkdata .eq. 0)then
-! kmt(iCell) = 0
+ elseif( thkdata .eq. 0)then
+ kmt(iCell) = 0
endif
enddo
-! ! output no. of cells removed
-! allocate(work_kmt(nCells))
-! work_kmt = 0.0
-! where(kmt.eq.0) work_kmt=1.0
-! write(6,*) 'number of cells culled ',sum(work_kmt)
-! deallocate(work_kmt)
+ ! output no. of cells removed
+ allocate(work_kmt(nCells))
+ work_kmt = 0.0
+ where(kmt.eq.0) work_kmt=1.0
+ write(6,*) 'number of cells culled ',sum(work_kmt)
+ deallocate(work_kmt)
deallocate(x)
deallocate(y)
- deallocate(thk)
+ deallocate(thck)
endif
-! eliminate isolated ocean cells
+! eliminate isolated cells
do iCell=1,nCells
flag = .true.
if(kmt(iCell).ne.0) then
@@ -1630,8 +1637,8 @@
integer k
dz( 1) = 0.0 ! 5.006218 10.01244
- dz( 2) = 50.0 ! 15.06873 20.12502
- dz( 3) = 100.0 ! 25.28342 30.44183
+ dz( 2) = 2000.0 ! 15.06873 20.12502
+ dz( 3) = 4000.0 ! 25.28342 30.44183
! dz( 1) = 1001.244 ! 5.006218 10.01244
! dz( 2) = 1011.258 ! 15.06873 20.12502
! dz( 3) = 1031.682 ! 25.28342 30.44183
@@ -1673,7 +1680,8 @@
! dz(39) = 25000.00 ! 5124.999 5249.999
! dz(40) = 25000.00 ! 5374.999 5499.999
- dz = dz / 100.0
+! dz = dz / 100.0
+ dz = dz / 1.0
write(6,*)
! do k=1,40
Modified: branches/land_ice/icesheet/src/module_read_topo.F
===================================================================
--- branches/land_ice/icesheet/src/module_read_topo.F 2011-07-14 15:23:34 UTC (rev 921)
+++ branches/land_ice/icesheet/src/module_read_topo.F 2011-07-14 23:04:48 UTC (rev 922)
@@ -3,12 +3,10 @@
integer :: rd_ncid
integer :: rdDimIDnx
integer :: rdDimIDny
- integer :: rdVarIDz
+! integer :: rdVarIDz
integer :: rdVarIDx
integer :: rdVarIDy
-! integer :: rdVarIDthk
-! integer :: rdVarIDlat
-! integer :: rdVarIDlon
+ integer :: rdVarIDthk
integer :: rdLocalnx
integer :: rdLocalny
@@ -24,54 +22,47 @@
integer, intent(out) :: nx, ny
integer :: nferr
+
+ integer, dimension(3) :: dimids
- nferr = nf_open('topo/Greenland_5km_v1.1.nc', NF_SHARE, rd_ncid)
- write(6,*) ' nferr ', nferr, rd_ncid
+! nferr = nf_open('topo/Greenland_5km_v1.1.nc', NF_SHARE, rd_ncid)
+! nferr = nf_open('topo/gis_5km.180511.nc', NF_SHARE, rd_ncid)
+ nferr = nf_open('topo/gis_20km.180511.nc', NF_SHARE, rd_ncid)
+ write(6,*) ' nferr ncid ', nferr, rd_ncid
+
!
! Get IDs for variable dimensions
!
-! nferr = nf_inq_dimid(rd_ncid, 'x', rdDimIDnx)
-! write(6,*) ' nferr ', nferr, rdDimIDnx
-! nferr = nf_inq_dimlen(rd_ncid, rdDimIDnx, rdLocalnx)
-! write(6,*) ' nferr ', nferr, rdLocalnx
-! nferr = nf_inq_dimid(rd_ncid, 'y', rdDimIDny)
-! write(6,*) ' nferr ', nferr, rdDimIDny
-! nferr = nf_inq_dimlen(rd_ncid, rdDimIDny, rdLocalny)
-! write(6,*) ' nferr ', nferr, rdLocalny
-
nferr = nf_inq_dimid(rd_ncid, 'x1', rdDimIDnx)
- write(6,*) ' nferr ', nferr, rdDimIDnx
+ write(6,*) ' nferr xdim-id ', nferr, rdDimIDnx
nferr = nf_inq_dimlen(rd_ncid, rdDimIDnx, rdLocalnx)
- write(6,*) ' nferr ', nferr, rdLocalnx
+ write(6,*) ' nferr xdim-length ', nferr, rdLocalnx
nferr = nf_inq_dimid(rd_ncid, 'y1', rdDimIDny)
- write(6,*) ' nferr ', nferr, rdDimIDny
+ write(6,*) ' nferr ydim-id ', nferr, rdDimIDny
nferr = nf_inq_dimlen(rd_ncid, rdDimIDny, rdLocalny)
- write(6,*) ' nferr ', nferr, rdLocalny
+ write(6,*) ' nferr ydim-length ', nferr, rdLocalny
nx = rdLocalnx
ny = rdLocalny
-
- write(6,*) nx, ny
+ print *, 'x,y dims of .nc var = ', nx, ny
!
! Get IDs for variables
!
-! nferr = nf_inq_varid(rd_ncid, 'x', rdVarIDx)
-! write(6,*) ' nferr ', nferr, rdVarIDx
-! nferr = nf_inq_varid(rd_ncid, 'y', rdVarIDy)
-! write(6,*) ' nferr ', nferr, rdVarIDy
-! nferr = nf_inq_varid(rd_ncid, 'z', rdVarIDz)
-! write(6,*) ' nferr ', nferr, rdVarIDz
nferr = nf_inq_varid(rd_ncid, 'x1', rdVarIDx)
- write(6,*) ' nferr ', nferr, rdVarIDx
+ write(6,*) ' nferr xvar-id ', nferr, rdVarIDx
nferr = nf_inq_varid(rd_ncid, 'y1', rdVarIDy)
- write(6,*) ' nferr ', nferr, rdVarIDy
- nferr = nf_inq_varid(rd_ncid, 'thk', rdVarIDz)
- write(6,*) ' nferr ', nferr, rdVarIDz
+ write(6,*) ' nferr yvar-id ', nferr, rdVarIDy
+ nferr = nf_inq_varid(rd_ncid, 'thk', rdVarIDthk)
+ write(6,*) ' nferr thkvar-id ', nferr, rdVarIDthk
+
+ nferr = nf_inq_vardimid( rd_ncid, rdVarIDthk, dimids )
+ write(6,*) ' nferr thk-dim-ids ', nferr, dimids
+
end subroutine read_topo_init
@@ -83,9 +74,6 @@
include 'netcdf.inc'
real (kind=4), dimension(:), intent(out) :: x,y
-! real (kind=4), dimension(:,:), intent(out) :: z
-! real (kind=4), dimension(:,:), intent(out) :: x
-! real (kind=4), dimension(:,:), intent(out) :: y
real (kind=4), dimension(:,:), intent(out) :: thk
integer, dimension(1) :: start1, count1
@@ -98,41 +86,28 @@
start1(1) = 1
count1(1) = rdLocalnx
nferr = nf_get_vara_real(rd_ncid, rdVarIDx, start1, count1, x)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDx
+ write(6,*) ' nferr ncid xvar-id (get x-val) ', nferr, rd_ncid, rdVarIDx
start1(1) = 1
count1(1) = rdLocalny
nferr = nf_get_vara_real(rd_ncid, rdVarIDy, start1, count1, y)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDy
+ write(6,*) ' nferr ncid yvar-id (get y-val) ', nferr, rd_ncid, rdVarIDy
start2(1) = 1
start2(2) = 1
count2(1) = rdLocalnx
count2(2) = rdLocalny
- nferr = nf_get_vara_real(rd_ncid, rdVarIDz, start2, count2, thk)
- write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDz, rdLocalnx
-
-! start2(1) = 1
-! start2(2) = 1
-! count2(1) = rdLocalnx
-! count2(2) = rdLocalny
-! nferr = nf_get_vara_real(rd_ncid, rdVarIDx, start2, count2, x)
-! write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDx, rdLocalnx
-!
-! start2(1) = 1
-! start2(2) = 1
-! count2(1) = rdLocalnx
-! count2(2) = rdLocalny
-! nferr = nf_get_vara_real(rd_ncid, rdVarIDy, start2, count2, y)
-! write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDy, rdLocalnx
-!
-! start2(1) = 1
-! start2(2) = 1
-! count2(1) = rdLocalnx
-! count2(2) = rdLocalny
-! nferr = nf_get_vara_real(rd_ncid, rdVarIDthk, start2, count2, thk)
-! write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDthk, rdLocalnx
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+ count3(1) = rdlocalnx
+ count3(2) = rdlocalny
+ count3(3) = 1
+
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDthk, start3, count3, thk)
+ write(6,*) ' nferr ncid thkvar-id (get thk val) ', nferr, rd_ncid, rdVarIDthk
+
end subroutine read_topo_fields
Modified: branches/land_ice/icesheet/src/module_write_netcdf.F
===================================================================
--- branches/land_ice/icesheet/src/module_write_netcdf.F 2011-07-14 15:23:34 UTC (rev 921)
+++ branches/land_ice/icesheet/src/module_write_netcdf.F 2011-07-14 23:04:48 UTC (rev 922)
@@ -107,7 +107,7 @@
wrLocalnVertLevels = nVertLevels
wrLocalvertexDegree = vertexDegree
- nferr = nf_create('ocean.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), wr_ncid)
+ nferr = nf_create('icesheet.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), wr_ncid)
!
! Define dimensions
</font>
</pre>