- 积分
- 315
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-4-27
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
program month_ALT
implicit none
include 'netcdf.inc'
integer*4 ncid, status, t, year ! file control
integer i, j, d, m, n
real delta, g
real, dimension(:), allocatable :: lon, lat, time, depth, new_depth, retime
real, dimension(:,:,:,:) :: SoilTemp(720,123,30,492), air(720,123,30,12), clm45_SoilTemp(720,123,290)
real, dimension(:,:,:) :: month_ALT(720,123,492)
integer :: tempid, latDimID, lonDimID, frTimeDimID, &
latVarID, lonVarID, frTimeVarID, ALTVarID
integer*4 :: start(10)
integer*4 :: count(10)
integer :: dimids(10)! allow up to 10 dimensions
integer :: dimid, xtype
integer :: ndim, natts, len
character(len = 31) :: dummy
!21
character(44) :: filename = '/p15/PCNMIP/rcn_data/jules/R01_JULES_1960.nc'
status = nf_open(filename,nf_nowrite,ncid)
!===========================================================
write(*,*) "read latitude"
status = nf_inq_var(ncid,2,dummy,xtype,ndim,dimids,natts)
allocate(lat(123))
do j = 1,ndim
status = nf_inq_dim(ncid,dimids(j),dummy,len)
if ( status /= nf_noerr ) write (*,*) nf_strerror(status)
start(j) = 1 ; count(j) = len
enddo
status = nf_get_vara_real(ncid, 2,start,count,lat)
! write(*,*) lat
!===============================================================
write(*,*) "read longitude"
status = nf_inq_var(ncid,1,dummy,xtype,ndim,dimids,natts)
allocate(lon(720))
do j = 1,ndim
status = nf_inq_dim(ncid,dimids(j),dummy,len)
if ( status /= nf_noerr ) write (*,*) nf_strerror(status)
start(j) = 1 ; count(j) = len
enddo
status = nf_get_vara_real(ncid, 1,start,count,lon)
! write(*,*) lon
!===================================================================
write(*,*) "read depth"
status = nf_inq_varid (ncid, "z_node", tempid)
status = nf_inq_var(ncid,tempid,dummy,xtype,ndim,dimids,natts)
allocate(depth(30))
do j = 1,ndim
status = nf_inq_dim(ncid,dimids(j),dummy,len)
if ( status /= nf_noerr ) write (*,*) nf_strerror(status)
start(j) = 1 ; count(j) = len
enddo
status = nf_get_vara_real(ncid, tempid,start,count,depth)
! write(*,*) depth
!=======================================================================
! 循环读取文件及SoilTemp
write(*,*) "read time and SoilTemp"
allocate(time(492))
allocate(retime(12))
do t = 0,40
year = t + 1960
write(filename(38:41),'(i4)') year
write(*,*) year
status = nf_open(filename,nf_nowrite,ncid)
!===============================================================
! write(*,*) "read time"
status = nf_inq_var(ncid,5,dummy,xtype,ndim,dimids,natts)
do j = 1,ndim
status = nf_inq_dim(ncid,dimids(j),dummy,len)
if ( status /= nf_noerr ) write (*,*) nf_strerror(status)
start(j) = 1 ; count(j) = len
enddo
status = nf_get_vara_real(ncid, 5,start,count,retime)
time(t*12+1:t*12+12) = retime(:)
!===================================================================
! write(*,*) "read SoilTemp"
status = nf_inq_varid (ncid, "SoilTemp", tempid)
status = nf_inq_var(ncid,tempid,dummy,xtype,ndim,dimids,natts)
do j = 1,ndim
status = nf_inq_dim(ncid,dimids(j),dummy,len)
if ( status /= nf_noerr ) write (*,*) nf_strerror(status)
start(j) = 1 ; count(j) = len
enddo
status = nf_get_vara_real(ncid,tempid,start,count,air)
! write(*,*) minval(air)
SoilTemp(:,:,:,t*12+1:t*12+12) = air(:,:,:,:)
enddo
write(*,*) maxval(SoilTemp)
!=============================================================================
! 创建新的深度数组
allocate(new_depth(290))
g = 0.05
do n = 1, 290
new_depth(n) = g + 0.01 * (n - 1)
end do
! write(*,*) 'new_depth = ', new_depth
!==============================================================================
! 土壤温度对深度线性插值
do t = 1,492
write(*,*) t
clm45_SoilTemp(:,:,1) = SoilTemp(:,:,1,t)
do i = 1,123
do j = 1,720
do d = 2,290
do n = 2, 30
! if ( SoilTemp(j,i,n,t) < 1.00000002E+20 .and. SoilTemp(j,i,n-1,t) < 1.00000002E+20 ) then
if ( (new_depth(d) - depth(n)) < 0 .and. (new_depth(d) - depth(n-1)) >0) then
clm45_SoilTemp(j,i,d) = SoilTemp(j,i,n-1,t) + (new_depth(d) - depth(n-1)) * &
((SoilTemp(j,i,n,t) - SoilTemp(j,i,n-1,t)) / (depth(n) - depth(n-1)))
end if
if ( clm45_SoilTemp(j,i,d) > 273.15 ) then
month_ALT(j,i,t) = new_depth(d)
end if
! end if
end do
end do
end do
end do
end do
|
-
编译出错
|