<p><b>ringler@lanl.gov</b> 2011-04-26 16:52:40 -0600 (Tue, 26 Apr 2011)</p><p><br>
add code for creation of real ocean basin with WOCE T,S,u used as initial conditions<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin/src/Makefile
===================================================================
--- branches/ocean_projects/basin/src/Makefile 2011-04-22 03:33:28 UTC (rev 801)
+++ branches/ocean_projects/basin/src/Makefile 2011-04-26 22:52:40 UTC (rev 802)
@@ -45,11 +45,12 @@
utilities.o \
module_read_netcdf.o \
module_read_topo.o \
+ module_read_TS.o \
module_write_netcdf.o
all: map
-basin.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_read_topo.o
+basin.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_read_topo.o module_read_TS.o
map: $(OBJS)
$(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
Modified: branches/ocean_projects/basin/src/basin.F
===================================================================
--- branches/ocean_projects/basin/src/basin.F 2011-04-22 03:33:28 UTC (rev 801)
+++ branches/ocean_projects/basin/src/basin.F 2011-04-26 22:52:40 UTC (rev 802)
@@ -2,6 +2,7 @@
use read_netcdf
use read_topo
+use read_TS
use write_netcdf
use utilities
@@ -48,11 +49,17 @@
real, allocatable, dimension(:,:,:) :: circulation, vorticity, ke, rho
real, allocatable, dimension(:,:,:,:) :: tracers
-real, dimension(25) :: dz
+integer nlon, nlat, ndepth
+real(kind=4), allocatable, dimension(:) :: t_lon, t_lat, depth_t
+real(kind=4), allocatable, dimension(:,:) :: mTEMP, mSALT
+real(kind=4), allocatable, dimension(:,:,:) :: TEMP, SALT
+real(kind=4), allocatable, dimension(:,:) :: TAUX, TAUY
+real, dimension(40) :: dz
+
! Step 1: Set the number of Vertical levels, and Tracers
-integer, parameter :: nVertLevelsMOD = 20
-integer, parameter :: nTracersMod = 4
+integer, parameter :: nVertLevelsMOD = 40
+integer, parameter :: nTracersMod = 2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! basin-mod
@@ -66,22 +73,24 @@
! Step 2: Set if the grid is on a sphere or not, and it's radius
-!character (len=16) :: on_a_sphere = 'YES '
-!real*8, parameter :: sphere_radius = 6.37122e6
+character (len=16) :: on_a_sphere = 'YES '
+real*8, parameter :: sphere_radius = 6.37122e6
!real*8, parameter :: sphere_radius = 1.0
-character (len=16) :: on_a_sphere = 'NO '
-real*8, parameter :: sphere_radius = 0.0
+!character (len=16) :: on_a_sphere = 'NO '
+!real*8, parameter :: sphere_radius = 0.0
-logical, parameter :: real_bathymetry=.false.
+logical, parameter :: real_bathymetry=.true.
+logical, parameter :: l_woce = .true.
+
! Step 3: Specify some Parameters
real (kind=8), parameter :: &
h_total_max = 2000.0, &
u_max = 0.0, &
- u_src_max = 10.0, & ! max wind stress, N/m2
- beta = 2.0e-11, &
- f0 = 1.4e-4, &
+ u_src_max = 0.1, & ! max wind stress, N/m2
+ beta = 1.4e-11, &
+ f0 = -1.1e-4, &
omega = 7.29212e-5
real (kind=8) :: ymid, ytmp, ymax, xmid, xloc, yloc, pert, ymin, distance, r, c1(3), c2(3)
@@ -108,24 +117,49 @@
real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, h_sNew
real, allocatable, dimension(:,:) :: u_srcNew
real, allocatable, dimension(:,:,:) :: uNew, vNew, hNew, vhNew
-real, allocatable, dimension(:,:,:) :: circulationNew, vorticityNew, keNew, rhoNew, temperatureNew, salinityNew, tracer1, tracer2
+real, allocatable, dimension(:,:,:) :: circulationNew, vorticityNew, keNew, rhoNew, temperatureNew, salinityNew
+real, allocatable, dimension(:) :: temperatureRestoreNew, salinityRestoreNew
real, allocatable, dimension(:,:,:,:) :: tracersNew
! mapping variables
-integer, allocatable, dimension(:) :: kmt, maxLevelCell
+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
integer :: iCell, iCell1, iCell2, iCell3, iEdge, iVertex, iVertex1, iVertex2
-integer :: iCellNew, iEdgeNew, iVertexNew, jCell1, jCell2
-real :: xin, yin, zin, ulon, ulat, ux, uy, uz, rlon, rlat
+integer :: iCellNew, iEdgeNew, iVertexNew, ndata, jCell1, jCell2, jCell, iter
+real :: xin, yin, zin, ulon, ulat, ux, uy, uz, rlon, rlat, temp_t, temp_s
! get to work
write(6,*) ' starting'
write(6,*)
+write(6,*) ' getting woce t and s '
+call read_TS_init(nlon, nlat, ndepth)
+write(6,*) ' TS INIT ', nlon, nlat, ndepth
+allocate(t_lon(nlon), t_lat(nlat), depth_t(ndepth), TEMP(nlon,nlat,ndepth), SALT(nlon,nlat,ndepth))
+allocate(TAUX(nlon,nlat), TAUY(nlon,nlat))
+allocate(mTEMP(nlat,ndepth), mSALT(nlat,ndepth))
+call read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
+call read_TS_finalize()
+do k=1,ndepth
+ ndata = 0; temp_t=0; temp_s=0
+ do j=1,nlat
+ do i=1,nlon
+ if(TEMP(i,j,k).gt.-10.0) then
+ ndata = ndata + 1
+ temp_t = temp_t + TEMP(i,j,k)
+ temp_s = temp_s + SALT(i,j,k)
+ endif
+ enddo
+ enddo
+ mTEMP(:,k) = temp_t / float(ndata)
+ mSALT(:,k) = temp_s / float(ndata)
+ write(6,*) ndata,mTemp(1,k),mSalt(1,k)
+enddo
+
! get depth profile for later
call get_dz
@@ -201,7 +235,7 @@
cellsOnVertexNew, &
edgesOnCellNew, &
areaCellNew, &
- maxLevelCell, &
+ maxLevelCellNew, &
depthCell, &
kiteAreasOnVertexNew )
@@ -251,9 +285,10 @@
!Step 4: Check the Initial conditions routine get_init_conditions
subroutine get_init_conditions
implicit none
-real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r
+real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r, temp_t, temp_s
real :: dotProd
-integer :: iTracer
+integer :: iTracer, ix, iy, ndata, i, j, k, ixt, iyt, ncull, jcount, iNoData, kdata(nVertLevelsMod)
+logical :: flag_lat
pi = 4.0*atan(1.0)
dtr = pi/180.0
@@ -293,8 +328,8 @@
hNew(1,3,:) = 3250.0
h_sNew(:) = -( hNew(1,1,:) + hNew(1,2,:) + hNew(1,3,:) )
else
-! hNew(1,1,:) = 3250.0
-! h_sNew(:) = -( hNew(1,1,:) )
+ hNew(1,1,:) = 3250.0
+ h_sNew(:) = -( hNew(1,1,:) )
endif
! basin-mod
@@ -358,7 +393,6 @@
endif
enddo
else
- if(nVertLevelsMOD .eq. 3) then
ymin = minval(yCellNew)
ymax = maxval(yCellNew)
r = 3.0e5
@@ -367,29 +401,16 @@
ytmp = yCellNew(i)
pert = exp(-(ytmp-ymid)**2/((0.1*r)**2))
k=1
- salinityNew(1,k,:) = pert
- pert = exp(-(ytmp-ymid)**2/((1.0*r)**2))
- k=2
- salinityNew(1,k,:) = pert
- pert = exp(-(ytmp-ymid)**2/((10.0*r)**2))
- k=3
- salinityNew(1,k,:) = pert
+ if(nVertLevelsMOD .eq. 3) then
+ salinityNew(1,k,:) = pert
+ pert = exp(-(ytmp-ymid)**2/((1.0*r)**2))
+ k=2
+ salinityNew(1,k,:) = pert
+ pert = exp(-(ytmp-ymid)**2/((10.0*r)**2))
+ k=3
+ salinityNew(1,k,:) = pert
+ endif
enddo
- else
-
- ! Set tracer info for double-gyre basin:
- do k=1,nVertLevelsMod
- temperatureNew(1,k,:) = 31-k
- salinityNew(1,k,:) = 32 + k*0.32
- enddo
- tracer1=1.0
- tracer2=1.0
- tracer1(1,1:10,:) = 0.0
- do k=1,nVertLevelsMod,2
- tracer2(1, k ,:) = 0.
- enddo
- endif
-
endif
! basin-mod
@@ -477,16 +498,6 @@
if(real_bathymetry) then
-halfwidth = 10.0*dtr
-h_sNew = 0.0
-do iCell=1,nCellsNew
- rlon = lonCellNew(iCell) - pi
- rlat = latCellNew(iCell)
- hNew(:,:,iCell) = 100.0 + 10.0*exp(-(rlon**2+rlat**2)/halfwidth**2)
- tracersNew(:,1,:,iCell) = hNew(:,:,iCell)
- if(nTracersMod.gt.1) tracersNew(:,2,:,iCell) = 1.0
- write(50,*) hNew(1,1,iCell)
-enddo
u_srcNew = 0.0
do iEdge=1,nEdgesNew
@@ -495,46 +506,133 @@
zin = zEdgeNew(iEdge)
rlon = lonEdgeNew(iEdge)/dtr
rlat = latEdgeNew(iEdge)/dtr
- stress = -0.1*cos(3.0*latEdgeNew(iEdge))
- write(6,*) rlat, stress
- ulon = stress
- ulat = 0.0
+ ix = nint(rlon/0.1 - 0.05) + nlon/2 + 1
+ ix = mod(ix,nlon)+1
+ iy = nlat
+ do jcount=1,nlat
+ if(t_lat(jcount).gt.rlat) then
+ iy = jcount
+ exit
+ endif
+ enddo
+ !stress = -0.1*cos(3.0*latEdgeNew(iEdge))
+ !ulon = stress
+ !ulat = 0.0
+ ulon = TAUX(ix,iy)
+ ulat = TAUY(ix,iy)
+ write(6,*) rlon, t_lon(ix), rlat, t_lat(iy)
call transform_from_lonlat_to_xyz(xin,yin,zin,ulon,ulat,ux,uy,uz)
if(boundaryEdgeNew(1,iEdge).eq.1) then
u_srcNew(1,iEdge) = 0.0
else
iCell1 = cellsOnEdgeNew(1,iEdge)
iCell2 = cellsOnEdgeNew(2,iEdge)
- p(1) = xCellNew(iCell1); p(2) = yCellNew(iCell1); p(3) = zCellNew(iCell1)
+ p(1) = xCellNew(iCell1); p(2) = yCellNew(iCell1); p(3) = zCellNew(iCell1)
q(1) = xCellNew(iCell2); q(2) = yCellNew(iCell2); q(3) = zCellNew(iCell2)
q = q - p
call unit_vector_in_3space(q)
u_srcNew(1,iEdge) = ux*q(1) + uy*q(2) + uz*q(3)
endif
enddo
+
-!set tracers
-! set all layer temps to 10.
-! set simple salinity
-!tracersNew = 0.0
-!do k = 1,nVertLevelsMod
- !tracersNew(:,1,k,:) = 10.0
- !tracersNew(:,2,k,:) = 1.4 + k*0.6 ! salinity
-!enddo
+!set tracers at a first guess
+tracersNew = 0.0
+temperatureNew = -99.0
+salinityNew = -99.0
+do iCell=1,nCellsNew
+do k = 1,maxLevelCellNew(iCell)
+ temperatureNew(1,k,iCell) = 20.0 - 10.0*k/nVertLevelsMod
+ salinityNew(1,k,iCell) = 34.0 ! salinity
+enddo
+enddo
-!set rho
-!rhoNew = 1000.0
+! update T and S field with WOCE data
+if(l_woce) then
+iNoData = 0
+do iCell=1,nCellsNew
+ if(mod(iCell,100).eq.0) write(6,*) 'l_woce t and s',iCell
+ rlon = lonCellNew(iCell)/dtr
+ rlat = latCellNew(iCell)/dtr
+ ix = nint(rlon/0.1 - 0.05) + nlon/2 + 1
+ ix = mod(ix,nlon)+1
+ iy = nlat
+ do j=1,nlat
+ if(t_lat(j).gt.rlat) then
+ iy = j
+ exit
+ endif
+ enddo
+ do k=1,maxLevelCellNew(iCell)
+ ndata = 0; temp_t = 0; temp_s = 0; kdata(:) = 0
+ do i=-50,50
+ ixt = ix + 8*i
+ if(ixt.lt.1) then
+ ixt = ixt + nlon
+ elseif(ixt.gt.nlon) then
+ ixt = ixt - nlon
+ endif
+ do j=-50,50
+ iyt = iy + 8*j
+ flag_lat = .true.
+ if(iyt.lt.1.or.iyt.gt.nlat) then
+ iyt = 1
+ flag_lat = .false.
+ endif
+ if(TEMP(ixt,iyt,k).gt.-10.0.and.flag_lat) then
+ ndata = ndata + 1
+ temp_t = temp_t + TEMP(ixt,iyt,k)
+ temp_s = temp_s + SALT(ixt,iyt,k)
+ endif
+ enddo
+ enddo
-!rlon = 0
-!do k=1,nVertLevelsMod
- !hNew(:,k,:) = dz(k)
- !rlon = rlon + dz(k)
-!enddo
-!h_sNew(:) = -rlon
+ if(ndata.gt.0) then
+ temperatureNew(1,k,iCell) = temp_t / float(ndata)
+ salinityNEW(1,k,iCell) = temp_s / float(ndata)
+ kdata(k) = 1
+ else
+ if(k.eq.1) iNoData = iNoData + 1
+ if(k.ge.3) then
+ if(kdata(k-1).eq.1) maxLevelCellNew(iCell) = k-1
+ endif
+ endif
+ enddo
+enddo
+
+! do a couple of smoothing passes
+do iter=1,3
+do iCell=1,nCellsNew
+do k=1,maxLevelCellNew(iCell)
+ ndata=1
+ temp_t = temperatureNew(1,k,iCell)
+ temp_s = salinityNew(1,k,iCell)
+ do j=1,nEdgesOnCellNew(iCell)
+ jCell = cellsOnCellNew(j,iCell)
+ if(maxLevelCellNew(jCell).ge.k) then
+ temp_t = temp_t + temperatureNew(1,k,jCell)
+ temp_s = temp_s + salinityNew(1,k,jCell)
+ ndata = ndata + 1
+ endif
+ enddo
+ temperatureNew(1,k,iCell) = temp_t / ndata
+ salinityNew(1,k,iCell) = temp_s / ndata
+enddo
+enddo
+write(6,*) maxval(temperatureNew(1,1,:)),maxval(salinityNew(1,1,:))
+enddo
+
+write(6,*) iNoData, nCellsNew
+
+temperatureRestoreNew(:) = temperatureNew(1,1,:)
+salinityRestoreNew(:) = salinityNew(1,1,:)
+
endif
+endif
+
write(6,*) ' done get_init_conditions'
end subroutine get_init_conditions
@@ -944,7 +1042,7 @@
yVertexNew, &
zVertexNew, &
indexToVertexIDNew, &
- maxLevelCell, &
+ maxLevelCellNew, &
cellsOnEdgeNew, &
nEdgesOnCellNew, &
nEdgesOnEdgeNew, &
@@ -978,7 +1076,9 @@
rhoNew, &
tracersNew, &
temperatureNew, &
- salinityNew,tracer1,tracer2 &
+ salinityNew, &
+ temperatureRestoreNew, &
+ salinityRestoreNew &
)
call write_netcdf_finalize
@@ -1008,7 +1108,7 @@
real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
real (kind=4), allocatable, dimension(:,:) :: ztopo
integer :: nx, ny, inx, iny, ix, iy
-real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax, xmin, xmax
+real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax
real :: latmin, latmax, lonmin, lonmax
logical :: flag, kmt_flag
pi = 4.0*atan(1.0)
@@ -1032,22 +1132,12 @@
where(lonCell.lt.lonmin) kmt = 0
where(lonCell.gt.lonmax) kmt = 0
else
- ! solid boundary in y
ymin = minval(yCell)
write(6,*) ' minimum yCell ', ymin
ymax = maxval(yCell)
write(6,*) ' maximum yCell ', ymax
where(yCell.lt.1.001*ymin) kmt = 0
where(yCell.gt.0.999*ymax) kmt = 0
-
- ! solid boundary in x
- xmin = minval(xCell)
- write(6,*) ' minimum xCell ', xmin
- xmax = maxval(xCell)
- write(6,*) ' maximum xCell ', xmax
- where(xCell.lt.1.001*xmin) kmt = 0
- where(xCell.gt.0.999*xmax) kmt = 0
-
endif
@@ -1238,7 +1328,7 @@
allocate(angleEdgeNew(nEdgesNew))
allocate(areaCellNew(nCellsNew))
allocate(areaTriangleNew(nVerticesNew))
-allocate(maxLevelCell(nCellsNew))
+allocate(maxLevelCellNew(nCellsNew))
allocate(depthCell(nCellsNew))
allocate(fEdgeNew(nEdgesNew))
@@ -1256,9 +1346,10 @@
allocate(tracersNew(1,nTracersNew,nVertLevelsNew,nCellsNew))
allocate(temperatureNew(1,nVertLevelsNew,nCellsNew))
allocate(salinityNew(1,nVertLevelsNew,nCellsNew))
-allocate(tracer1(1,nVertLevelsNew,nCellsNew))
-allocate(tracer2(1,nVertLevelsNew,nCellsNew))
+allocate(temperatureRestoreNew(nCellsNew))
+allocate(salinityRestoreNew(nCellsNew))
+
xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0
xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
@@ -1268,6 +1359,10 @@
vorticityNew=0
tracersNew=0; temperatureNew=0; salinityNew=0
+temperatureRestoreNew = 0.0
+salinityRestoreNew = 0.0
+
+
do i=1,nCells
jNew = cellMap(i)
if(jNew.ne.0) then
@@ -1277,7 +1372,7 @@
latCellNew(jNew)=latCell(i)
lonCellNew(jNew)=lonCell(i)
areaCellNew(jNew)=areaCell(i)
- maxLevelCell(jNew) = kmt(i)
+ maxLevelCellNew(jNew) = kmt(i)
depthCell(jNew) = kmt(i)
endif
enddo
@@ -1510,34 +1605,56 @@
subroutine get_dz
+integer k
-dz( 1) = 1200.000000000000
-dz( 2) = 1595.000000000000
-dz( 3) = 2095.000000000000
-dz( 4) = 2721.000000000000
-dz( 5) = 3493.000000000000
-dz( 6) = 4431.000000000000
-dz( 7) = 5558.000000000000
-dz( 8) = 6889.000000000000
-dz( 9) = 8442.000000000000
-dz(10) = 10225.00000000000
-dz(11) = 12241.00000000000
-dz(12) = 14486.00000000000
-dz(13) = 16944.00000000000
-dz(14) = 19591.00000000000
-dz(15) = 22390.00000000000
-dz(16) = 25294.00000000000
-dz(17) = 28244.00000000000
-dz(18) = 31174.00000000000
-dz(19) = 34011.00000000000
-dz(20) = 36677.00000000000
-dz(21) = 39096.00000000000
-dz(22) = 41193.00000000000
-dz(23) = 42902.00000000000
-dz(24) = 44166.00000000000
-dz(25) = 44942.00000000000
+ dz( 1) = 1001.244 ! 5.006218 10.01244
+ dz( 2) = 1011.258 ! 15.06873 20.12502
+ dz( 3) = 1031.682 ! 25.28342 30.44183
+ dz( 4) = 1063.330 ! 35.75848 41.07513
+ dz( 5) = 1107.512 ! 46.61269 52.15025
+ dz( 6) = 1166.145 ! 57.98098 63.81171
+ dz( 7) = 1241.928 ! 70.02135 76.23099
+ dz( 8) = 1338.612 ! 82.92405 89.61711
+ dz( 9) = 1461.401 ! 96.92412 104.2311
+ dz(10) = 1617.561 ! 112.3189 120.4067
+ dz(11) = 1817.368 ! 129.4936 138.5804
+ dz(12) = 2075.558 ! 148.9582 159.3360
+ dz(13) = 2413.680 ! 171.4044 183.4728
+ dz(14) = 2863.821 ! 197.7919 212.1110
+ dz(15) = 3474.644 ! 229.4842 246.8575
+ dz(16) = 4320.857 ! 268.4617 290.0660
+ dz(17) = 5516.812 ! 317.6501 345.2342
+ dz(18) = 7230.458 ! 381.3865 417.5388
+ dz(19) = 9674.901 ! 465.9133 514.2878
+ dz(20) = 13003.92 ! 579.3074 644.3270
+ dz(21) = 17004.89 ! 729.3514 814.3759
+ dz(22) = 20799.33 ! 918.3725 1022.369
+ dz(23) = 23356.94 ! 1139.154 1255.939
+ dz(24) = 24527.19 ! 1378.574 1501.210
+ dz(25) = 24898.04 ! 1625.701 1750.191
+ dz(26) = 24983.22 ! 1875.107 2000.023
+ dz(27) = 24997.87 ! 2125.012 2250.002
+ dz(28) = 24999.79 ! 2375.000 2500.000
+ dz(29) = 24999.98 ! 2625.000 2749.999
+ dz(30) = 25000.00 ! 2874.999 2999.999
+ dz(31) = 25000.00 ! 3124.999 3249.999
+ dz(32) = 25000.00 ! 3374.999 3499.999
+ dz(33) = 25000.00 ! 3624.999 3749.999
+ dz(34) = 25000.00 ! 3874.999 3999.999
+ dz(35) = 25000.00 ! 4124.999 4249.999
+ dz(36) = 25000.00 ! 4374.999 4499.999
+ dz(37) = 25000.00 ! 4624.999 4749.999
+ dz(38) = 25000.00 ! 4874.999 4999.999
+ dz(39) = 25000.00 ! 5124.999 5249.999
+ dz(40) = 25000.00 ! 5374.999 5499.999
-dz = dz / 100.0
+ dz = dz / 100.0
+ write(6,*)
+ do k=1,40
+ write(6,*) k,dz(k)
+ enddo
+ write(6,*)
+
end subroutine get_dz
end program map_to_basin
Added: branches/ocean_projects/basin/src/module_read_TS.F
===================================================================
--- branches/ocean_projects/basin/src/module_read_TS.F (rev 0)
+++ branches/ocean_projects/basin/src/module_read_TS.F 2011-04-26 22:52:40 UTC (rev 802)
@@ -0,0 +1,165 @@
+module read_TS
+
+ integer :: rd_ncid, rd_ncids, rd_ncidu
+ integer :: rdDimIDt_lon
+ integer :: rdDimIDt_lat
+ integer :: rdDimIDdepth_t
+ integer :: rdVarIDt_lon
+ integer :: rdVarIDt_lat
+ integer :: rdVarIDdepth_t
+ integer :: rdVarIDTEMP
+ integer :: rdVarIDSALT
+ integer :: rdVarIDTAUX
+ integer :: rdVarIDTAUY
+
+ integer :: rdLocalt_lon
+ integer :: rdLocalt_lat
+ integer :: rdLocaldepth_t
+
+ contains
+
+ subroutine read_TS_init(nx, ny, nz)
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(out) :: nx, ny, nz
+
+ integer :: nferr, nferrs, nferru
+
+ nferr = nf_open('TS/woce_t_ann.3600x2431x42interp.r4.nc', NF_SHARE, rd_ncid)
+ write(6,*) ' nferr ', nferr, rd_ncid
+
+ !
+ ! Get IDs for variable dimensions
+ !
+ nferr = nf_inq_dimid(rd_ncid, 't_lon', rdDimIDt_lon)
+ write(6,*) ' nferr ', nferr, rdDimIDt_lon
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDt_lon, rdLocalt_lon)
+ write(6,*) ' nferr ', nferr, rdLocalt_lon
+ nferr = nf_inq_dimid(rd_ncid, 't_lat', rdDimIDt_lat)
+ write(6,*) ' nferr ', nferr, rdDimIDt_lat
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDt_lat, rdLocalt_lat)
+ write(6,*) ' nferr ', nferr, rdLocalt_lat
+ nferr = nf_inq_dimid(rd_ncid, 'depth_t', rdDimIDdepth_t)
+ write(6,*) ' nferr ', nferr, rdDimIDdepth_t
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDdepth_t, rdLocaldepth_t)
+ write(6,*) ' nferr ', nferr, rdLocaldepth_t
+
+ nx = rdLocalt_lon
+ ny = rdLocalt_lat
+ nz = rdLocaldepth_t
+
+ write(6,*) nx, ny, nz
+
+ !
+ ! Get IDs for variables
+ !
+ nferr = nf_inq_varid(rd_ncid, 't_lon', rdVarIDt_lon)
+ write(6,*) ' nferr ', nferr, rdVarIDt_lon
+ nferr = nf_inq_varid(rd_ncid, 't_lat', rdVarIDt_lat)
+ write(6,*) ' nferr ', nferr, rdVarIDt_lat
+ nferr = nf_inq_varid(rd_ncid, 'depth_t', rdVarIDdepth_t)
+ write(6,*) ' nferr ', nferr, rdVarIDdepth_t
+ nferr = nf_inq_varid(rd_ncid, 'TEMP', rdVarIDTEMP)
+ write(6,*) ' nferr ', nferr, rdVarIDTEMP
+
+ nferrs = nf_open('TS/woce_s_ann.3600x2431x42interp.r4.nc', NF_SHARE, rd_ncids)
+ nferrs = nf_inq_varid(rd_ncids, 'SALT', rdVarIDSALT)
+ write(6,*) ' nferrs ', nferrs, rdVarIDSALT
+
+ nferru = nf_open('TS/ws.old_ncep_1958-2000avg.interp3600x2431.nc', NF_SHARE, rd_ncidu)
+ nferru = nf_inq_varid(rd_ncidu, 'TAUX', rdVarIDTAUX)
+ nferru = nf_inq_varid(rd_ncidu, 'TAUY', rdVarIDTAUY)
+ write(6,*) ' nferru ', nferru, rdVarIDTAUX, rdVarIDTAUY
+
+ end subroutine read_TS_init
+
+ subroutine read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ real (kind=4), dimension(:), intent(out) :: t_lon, t_lat, depth_t
+ real (kind=4), dimension(:,:,:), intent(out) :: TEMP, SALT
+ real (kind=4), dimension(:,:), intent(out) :: TAUX, TAUY
+
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ integer :: nferr, nferrs, nferru
+
+ start1(1) = 1
+ count1(1) = rdLocalt_lon
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDt_lon, start1, count1, t_lon)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDt_lon
+
+ start1(1) = 1
+ count1(1) = rdLocalt_lat
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDt_lat, start1, count1, t_lat)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDt_lat
+
+ start1(1) = 1
+ count1(1) = rdLocaldepth_t
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDdepth_t, start1, count1, depth_t)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDdepth_t
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+ count3(1) = rdLocalt_lon
+ count3(2) = rdLocalt_lat
+ count3(3) = rdLocaldepth_t
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDTEMP, start3, count3, TEMP)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDTEMP
+ write(6,*) ' temperature' , minval(TEMP), maxval(TEMP)
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+ count3(1) = rdLocalt_lon
+ count3(2) = rdLocalt_lat
+ count3(3) = rdLocaldepth_t
+ nferrs = nf_get_vara_real(rd_ncids, rdVarIDSALT, start3, count3, SALT)
+ write(6,*) ' nferrs ', nferrs, rd_ncids, rdVarIDSALT
+ write(6,*) ' salinity' , minval(SALT), maxval(SALT)
+
+ start2(1) = 1
+ start2(2) = 1
+ count2(1) = rdLocalt_lon
+ count2(2) = rdLocalt_lat
+ nferru = nf_get_vara_real(rd_ncidu, rdVarIDTAUX, start2, count2, TAUX)
+ nferru = nf_get_vara_real(rd_ncidu, rdVarIDTAUY, start2, count2, TAUY)
+ write(6,*) ' nferru ', nferru, rd_ncidu, rdVarIDTAUX, rdVarIDTAUY
+ write(6,*) ' TAUX' , minval(TAUX), maxval(TAUX)
+ write(6,*) ' TAUY' , minval(TAUY), maxval(TAUY)
+
+
+ end subroutine read_TS_fields
+
+
+ subroutine read_TS_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr, nferrs, nferru
+
+ nferr = nf_close(rd_ncid)
+ write(6,*) ' nferr ', nferr
+
+ nferrs = nf_close(rd_ncids)
+ write(6,*) ' nferrs ', nferrs
+
+ nferru = nf_close(rd_ncidu)
+ write(6,*) ' nferru ', nferru
+
+
+ end subroutine read_TS_finalize
+
+end module read_TS
Modified: branches/ocean_projects/basin/src/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/basin/src/module_write_netcdf.F 2011-04-22 03:33:28 UTC (rev 801)
+++ branches/ocean_projects/basin/src/module_write_netcdf.F 2011-04-26 22:52:40 UTC (rev 802)
@@ -64,8 +64,8 @@
integer :: wrVarIDtracers
integer :: wrVarIDtemperature
integer :: wrVarIDsalinity
- integer :: wrVarIDtracer1
- integer :: wrVarIDtracer2
+ integer :: wrVarIDtemperatureRestore
+ integer :: wrVarIDsalinityRestore
integer :: wrLocalnCells
integer :: wrLocalnEdges
@@ -223,6 +223,10 @@
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)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'temperatureRestore', NF_DOUBLE, 1, dimlist, wrVarIDtemperatureRestore)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'salinityRestore', NF_DOUBLE, 1, dimlist, wrVarIDsalinityRestore)
dimlist( 1) = wrDimIDnVertLevels
dimlist( 2) = wrDimIDnEdges
dimlist( 3) = wrDimIDTime
@@ -277,14 +281,6 @@
dimlist( 2) = wrDimIDnCells
dimlist( 3) = wrDimIDTime
nferr = nf_def_var(wr_ncid, 'salinity', NF_DOUBLE, 3, dimlist, wrVarIDsalinity)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'tracer1', NF_DOUBLE, 3, dimlist, wrVarIDtracer1)
- dimlist( 1) = wrDimIDnVertLevels
- dimlist( 2) = wrDimIDnCells
- dimlist( 3) = wrDimIDTime
- nferr = nf_def_var(wr_ncid, 'tracer2', NF_DOUBLE, 3, dimlist, wrVarIDtracer2)
nferr = nf_put_att_text(wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, on_a_sphere)
@@ -349,7 +345,9 @@
rho, &
tracers, &
temperature, &
- salinity,tracer1,tracer2 &
+ salinity, &
+ temperatureRestore, &
+ salinityRestore &
)
implicit none
@@ -410,8 +408,9 @@
real (kind=8), dimension(:,:,:,:), intent(in) :: tracers
real (kind=8), dimension(:,:,:), intent(in) :: temperature
real (kind=8), dimension(:,:,:), intent(in) :: salinity
- real (kind=8), dimension(:,:,:), intent(in) :: tracer1
- real (kind=8), dimension(:,:,:), intent(in) :: tracer2
+ real (kind=8), dimension(:), intent(in) :: temperatureRestore
+ real (kind=8), dimension(:), intent(in) :: salinityRestore
+
integer :: nferr
integer, dimension(1) :: start1, count1
@@ -613,6 +612,14 @@
count2( 1) = wrLocalnVertLevels
count2( 2) = wrLocalnEdges
nferr = nf_put_vara_double(wr_ncid, wrVarIDu_src, start2, count2, u_src)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDtemperatureRestore, start1, count1, temperatureRestore)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDsalinityRestore, start1, count1, salinityRestore)
start3(3) = time
count3( 1) = wrLocalnVertLevels
@@ -680,18 +687,6 @@
count3( 2) = wrLocalnCells
count3( 3) = 1
nferr = nf_put_vara_double(wr_ncid, wrVarIDsalinity, start3, count3, salinity)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDtracer1, start3, count3, tracer1)
-
- start3(3) = time
- count3( 1) = wrLocalnVertLevels
- count3( 2) = wrLocalnCells
- count3( 3) = 1
- nferr = nf_put_vara_double(wr_ncid, wrVarIDtracer2, start3, count3, tracer2)
end subroutine write_netcdf_fields
</font>
</pre>