登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 skykiss水仙 于 2015-12-7 17:18 编辑
想调用GRAPES模式模拟的数据,并做一些简单的运算,但是出现了这样一个问题PGFIO-F-202/unformatted read/unit=20/conflicting specifiers.
代码如下,还请各位大虾指教
PROGRAM Total
IMPLICIT NONE
REAL, PARAMETER :: t_kelvin = 273.15
REAl, PARAMETER :: es_alpha = 6.112
REAl, PARAMETER :: es_beta = 17.67
REAl, PARAMETER :: es_gamma = 243.5
REAL,dimension(1:26),PARAMETER :: p_ncep = (/1000.0, 975.0, 950.0, 925.0, 900.0, 850.0, 800.0, &
750.0, 700.0, 650.0, 600.0, 550.0, 500.0, 450.0, &
400.0, 350.0, 300.0, 250.0, 200.0, 150.0, 100.0, &
70.0, 50.0, 30.0, 20.0, 10.0/)
integer,parameter :: nx=360,ny=180,nz=26
REAL,dimension(nx,ny,nz) :: h0,t0,u0,v0,q0,rh0
REAL,dimension(nx,ny) :: ps0,ts0,temp
REAL,dimension(nx,ny+1,nz) :: h,t,u,v,rh
REAL,dimension(nx,ny+1) :: ps,ts
REAL,dimension(nx,ny+1,nz) :: h_a,t_a,u_a,v_a,rh_a
REAL,dimension(nx,ny+1) :: ps_a,ts_a
REAL,dimension(:,:,:,:),ALLOCATABLE :: h_m,t_m,u_m,v_m,rh_m,h_e,t_e,u_e,v_e,rh_e, &
h_m1,t_m1,u_m1,v_m1,rh_m1,h_e1,t_e1,u_e1,v_e1,rh_e1, &
h_m2,t_m2,u_m2,v_m2,rh_m2,h_e2,t_e2,u_e2,v_e2,rh_e2, &
h_m3,t_m3,u_m3,v_m3,rh_m3,h_e3,t_e3,u_e3,v_e3,rh_e3, &
h_m4,t_m4,u_m4,v_m4,rh_m4
REAL,dimension(:,:,:),ALLOCATABLE :: ps_m,ts_m,ps_e,ts_e, &
ps_m1,ts_m1,ps_e1,ts_e1, &
ps_m2,ts_m2,ps_e2,ts_e2, &
ps_m3,ts_m3,ps_e3,ts_e3, &
ps_m4,ts_m4
integer :: i,j,k,nrec,kk,fnlfid,ii,postfid,wfid
REAL :: tc, es,qs
logical :: ncep_up,write_error
data write_error/.true./
data ncep_up/.false./
CHARACTER(len=10),DIMENSION(:,:),ALLOCATABLE :: itime
integer :: nt,nsamples
integer(4) :: IARGC
character(len=10) :: char_ntime,char_nsamples,icase
CHARACTER(LEN=100) :: path_fnl,path_forcast
character(len=100) :: fname_itime,fname_forcast,fname_fnl,fnamew
path_fnl='/public2/users/lzuhuangjp/data/FNL'
! ------------------------------------
IF(IARGC() /= 3) THEN
print*, 'Usage: <execute file> <fname_itime> <nsamples> <ntime>'
stop
ENDIF
CALL getarg(1,fname_itime)
CALL getarg(2,char_nsamples)
CALL getarg(3,char_ntime)
read(char_ntime,*)nt
read(char_nsamples,*)nsamples
ALLOCATE(h_m(nx,ny+1,nz,nt))
ALLOCATE(t_m(nx,ny+1,nz,nt))
ALLOCATE(u_m(nx,ny+1,nz,nt))
ALLOCATE(v_m(nx,ny+1,nz,nt))
ALLOCATE(rh_m(nx,ny+1,nz,nt))
ALLOCATE(ps_m(nx,ny+1,nt))
ALLOCATE(ts_m(nx,ny+1,nt))
ALLOCATE(h_e(nx,ny+1,nz,nt))
ALLOCATE(t_e(nx,ny+1,nz,nt))
ALLOCATE(u_e(nx,ny+1,nz,nt))
ALLOCATE(v_e(nx,ny+1,nz,nt))
ALLOCATE(rh_e(nx,ny+1,nz,nt))
ALLOCATE(ps_e(nx,ny+1,nt))
ALLOCATE(ts_e(nx,ny+1,nt))
ALLOCATE(h_m1(nx,ny+1,nz,nt))
ALLOCATE(t_m1(nx,ny+1,nz,nt))
ALLOCATE(u_m1(nx,ny+1,nz,nt))
ALLOCATE(v_m1(nx,ny+1,nz,nt))
ALLOCATE(rh_m1(nx,ny+1,nz,nt))
ALLOCATE(ps_m1(nx,ny+1,nt))
ALLOCATE(ts_m1(nx,ny+1,nt))
ALLOCATE(h_e1(nx,ny+1,nz,nt))
ALLOCATE(t_e1(nx,ny+1,nz,nt))
ALLOCATE(u_e1(nx,ny+1,nz,nt))
ALLOCATE(v_e1(nx,ny+1,nz,nt))
ALLOCATE(rh_e1(nx,ny+1,nz,nt))
ALLOCATE(ps_e1(nx,ny+1,nt))
ALLOCATE(ts_e1(nx,ny+1,nt))
ALLOCATE(h_m2(nx,ny+1,nz,nt))
ALLOCATE(t_m2(nx,ny+1,nz,nt))
ALLOCATE(u_m2(nx,ny+1,nz,nt))
ALLOCATE(v_m2(nx,ny+1,nz,nt))
ALLOCATE(rh_m2(nx,ny+1,nz,nt))
ALLOCATE(ps_m2(nx,ny+1,nt))
ALLOCATE(ts_m2(nx,ny+1,nt))
ALLOCATE(h_e2(nx,ny+1,nz,nt))
ALLOCATE(t_e2(nx,ny+1,nz,nt))
ALLOCATE(u_e2(nx,ny+1,nz,nt))
ALLOCATE(v_e2(nx,ny+1,nz,nt))
ALLOCATE(rh_e2(nx,ny+1,nz,nt))
ALLOCATE(ps_e2(nx,ny+1,nt))
ALLOCATE(ts_e2(nx,ny+1,nt))
ALLOCATE(h_m3(nx,ny+1,nz,nt))
ALLOCATE(t_m3(nx,ny+1,nz,nt))
ALLOCATE(u_m3(nx,ny+1,nz,nt))
ALLOCATE(v_m3(nx,ny+1,nz,nt))
ALLOCATE(rh_m3(nx,ny+1,nz,nt))
ALLOCATE(ps_m3(nx,ny+1,nt))
ALLOCATE(ts_m3(nx,ny+1,nt))
ALLOCATE(h_e3(nx,ny+1,nz,nt))
ALLOCATE(t_e3(nx,ny+1,nz,nt))
ALLOCATE(u_e3(nx,ny+1,nz,nt))
ALLOCATE(v_e3(nx,ny+1,nz,nt))
ALLOCATE(rh_e3(nx,ny+1,nz,nt))
ALLOCATE(ps_e3(nx,ny+1,nt))
ALLOCATE(ts_e3(nx,ny+1,nt))
ALLOCATE(h_m4(nx,ny+1,nz,nt))
ALLOCATE(t_m4(nx,ny+1,nz,nt))
ALLOCATE(u_m4(nx,ny+1,nz,nt))
ALLOCATE(v_m4(nx,ny+1,nz,nt))
ALLOCATE(rh_m4(nx,ny+1,nz,nt))
ALLOCATE(ps_m4(nx,ny+1,nt))
ALLOCATE(ts_m4(nx,ny+1,nt))
OPEN(10,file='/gpfshome/huangjp/zhangmeng/program/time.txt',form='formatted',STATUS='unknown')
do i=1,nsamples
do j=1,nt
read(10,*) itime(i,j)
end do
end do
CLOSE(10)
path_forcast='/public2/users/lzuhuangjp/data/forcast_lib/post'
t_m=0.0
h_m=0.0
u_m=0.0
v_m=0.0
rh_m=0.0
ps_m=0.0
ts_m=0.0
t_m1=0.0
h_m1=0.0
u_m1=0.0
v_m1=0.0
rh_m1=0.0
ps_m1=0.0
ts_m1=0.0
t_m2=0.0
h_m2=0.0
u_m2=0.0
v_m2=0.0
rh_m2=0.0
ps_m2=0.0
ts_m2=0.0
t_m3=0.0
h_m3=0.0
u_m3=0.0
v_m3=0.0
rh_m3=0.0
ps_m3=0.0
ts_m3=0.0
t_m4=0.0
h_m4=0.0
u_m4=0.0
v_m4=0.0
rh_m4=0.0
ps_m4=0.0
ts_m4=0.0
do ii=1,nsamples
!!-------------------------------------read the postvar data----------------------------------------
fname_forcast=trim(path_forcast)//'/postvar'//trim(itime(ii,1))
postfid=20
OPEN(postfid,file=trim(fname_forcast),form='unformatted',access='sequential',STATUS='unknown')
write(*,*)'the transfer of:',itime(ii,1)
do kk=1,nt
do k=1,nz,1
read(postfid) ((u0(i,j,k),i=1,nx),j=1,ny)
enddo
do k=1,nz,1
read(postfid) ((v0(i,j,k),i=1,nx),j=1,ny)
enddo
do k=1,nz,1
read(postfid) ((t0(i,j,k),i=1,nx),j=1,ny)
enddo
do k=1,nz,1
read(postfid) ((h0(i,j,k),i=1,nx),j=1,ny)
enddo
do k=1,nz,1
read(postfid) ((q0(i,j,k),i=1,nx),j=1,ny)
enddo
do k=1,nz*6,1
read(postfid) ((temp(i,j),i=1,nx),j=1,ny)
enddo
read(postfid) ((ps0(i,j),i=1,nx),j=1,ny)
do k=1,3,1
read(postfid) ((temp(i,j),i=1,nx),j=1,ny)
enddo
read(postfid) ((ts0(i,j),i=1,nx),j=1,ny)
do k=1,5,1
read(postfid) ((temp(i,j),i=1,nx),j=1,ny)
enddo
do k=1,nz
do j=1,ny
do i=1,nx
tc=t0(i,j,k)-t_kelvin
es=es_alpha
qs=0.622*es/(p_ncep(k)-0.378*es)
rh0(i,j,k)=q0(i,j,k)*100.0/qs
if(rh0(i,j,k).lt.0.0) rh0(i,j,k)=0.0
if(rh0(i,j,k).gt.100.0) rh0(i,j,k)=100.0
enddo
enddo
enddo
DO k = 1,nz
DO i = 1,nx
DO j = 1,ny-1
u(i,j+1,k) =(u0(i,j,k)+u0(i,j+1,k))/2.0
v(i,j+1,k) =(v0(i,j,k)+v0(i,j+1,k))/2.0
t(i,j+1,k) =(t0(i,j,k)+t0(i,j+1,k))/2.0
h(i,j+1,k) =(h0(i,j,k)+h0(i,j+1,k))/2.0
rh(i,j+1,k) =(rh0(i,j,k)+rh0(i,j+1,k))/2.0
END DO
END DO
END DO
DO k = 1,nz
DO i = 1,nx
u(i,1,k) =sum(u0(:,1,k))/nx
v(i,1,k) =sum(v0(:,1,k))/nx
t(i,1,k) =sum(t0(:,1,k))/nx
h(i,1,k) =sum(h0(:,1,k))/nx
rh(i,1,k) =sum(rh0(:,1,k))/nx
u(i,ny+1,k) =sum(u0(:,ny,k))/nx
v(i,ny+1,k) =sum(v0(:,ny,k))/nx
t(i,ny+1,k) =sum(t0(:,ny,k))/nx
h(i,ny+1,k) =sum(h0(:,ny,k))/nx
rh(i,ny+1,k) =sum(rh0(:,ny,k))/nx
END DO
END DO
DO i = 1,nx
DO j = 1,ny-1
ps(i,j+1) = (ps0(i,j)+ps0(i,j+1))/2.0
ts(i,j+1) = (ts0(i,j)+ts0(i,j+1))/2.0
END DO
END DO
DO i = 1,nx
ps(i,1) =sum(ps0(:,1))/nx
ps(i,ny+1) =sum(ps0(:,ny))/nx
ts(i,1) =sum(ts0(:,1))/nx
ts(i,ny+1) =sum(ts0(:,ny))/nx
END DO
fnlfid=30
fname_fnl =trim(path_fnl)//'/FNL'//trim(itime(ii,kk))//'000'
call read_FNL_4d(fnlfid,fname_fnl,t_a,h_a,u_a,v_a,rh_a,ps_a,ts_a,nx,ny+1,nz,nt,ncep_up)
t_e1(:,:,:,kk)=(t_a-t)**2.0
h_e1(:,:,:,kk)=(h_a-h)**2.0
u_e1(:,:,:,kk)=(u_a-u)**2.0
v_e1(:,:,:,kk)=(v_a-v)**2.0
rh_e1(:,:,:,kk)=(rh_a-rh)**2.0
ps_e1(:,:,kk)=(ps_a/100.-ps)**2.0
ts_e1(:,:,kk)=(ts_a-ts)**2.0
t_e2(:,:,:,kk)=t_a
h_e2(:,:,:,kk)=h_a
u_e2(:,:,:,kk)=u_a
v_e2(:,:,:,kk)=v_a
rh_e2(:,:,:,kk)=rh_a
ps_e2(:,:,kk)=ps_a/100.0
ts_e2(:,:,kk)=ts_a
t_e3(:,:,:,kk)=t
h_e3(:,:,:,kk)=h
u_e3(:,:,:,kk)=u
v_e3(:,:,:,kk)=v
rh_e3(:,:,:,kk)=rh
ps_e3(:,:,kk)=ps
ts_e3(:,:,kk)=ts
end do
close(postfid)
t_m1=t_e1+t_m1
h_m1=h_e1+h_m1
u_m1=u_e1+u_m1
v_m1=v_e1+v_m1
rh_m1=rh_e1+rh_m1
ps_m1=ps_e1+ps_m1
ts_m1=ts_e1+ts_m1
t_m3=t_e2+t_m3
h_m3=h_e2+h_m3
u_m3=u_e2+u_m3
v_m3=v_e2+v_m3
rh_m3=rh_e2+rh_m3
ps_m3=ps_e2+ps_m3
ts_m3=ts_e2+ts_m3
t_m4=t_e3+t_m4
h_m4=h_e3+h_m4
u_m4=u_e3+u_m4
v_m4=v_e3+v_m4
rh_m4=rh_e3+rh_m4
ps_m4=ps_e3+ps_m4
ts_m4=ts_e3+ts_m4
t_m2=(t_m3-t_m4)**2.0
h_m2=(h_m3-h_m4)**2.0
u_m2=(u_m3-u_m4)**2.0
v_m2=(v_m3-v_m4)**2.0
rh_m2=(rh_m3-rh_m4)**2.0
ps_m2=(ps_m3-ps_m4)**2.0
ts_m2=(ts_m3-ts_m4)**2.0
end do
t_m1=t_m1/nsamples
h_m1=h_m1/nsamples
u_m1=u_m1/nsamples
v_m1=v_m1/nsamples
rh_m1=rh_m1/nsamples
ps_m1=ps_m1/nsamples
ts_m1=ts_m1/nsamples
t_m2=t_m2/(nsamples**2.0)
h_m2=h_m2/(nsamples**2.0)
u_m2=u_m2/(nsamples**2.0)
v_m2=v_m2/(nsamples**2.0)
rh_m2=rh_m2/(nsamples**2.0)
ps_m2=ps_m2/(nsamples**2.0)
ts_m2=ts_m2/(nsamples**2.0)
t_m=t_m1-t_m2
h_m=h_m1-h_m2
u_m=u_m1-u_m2
v_m=v_m1-v_m2
rh_m=rh_m1-rh_m2
ps_m=ps_m1-ps_m2
ts_m=ps_m1-ps_m2
if(write_error)then
fnamew ='Total_4d'
call write_FNL_4d(wfid,fnamew,t_m,h_m,u_m,v_m,rh_m,ps_m,ts_m,nx,ny+1,nz,nt,ncep_up)
end if
deallocate(itime)
deallocate(t_m)
deallocate(h_m)
deallocate(u_m)
deallocate(v_m)
deallocate(rh_m)
deallocate(ps_m)
deallocate(ts_m)
deallocate(t_e)
deallocate(h_e)
deallocate(u_e)
deallocate(v_e)
deallocate(rh_e)
deallocate(ps_e)
deallocate(ts_e)
deallocate(t_m1)
deallocate(h_m1)
deallocate(u_m1)
deallocate(v_m1)
deallocate(rh_m1)
deallocate(ps_m1)
deallocate(ts_m1)
deallocate(t_e1)
deallocate(h_e1)
deallocate(u_e1)
deallocate(v_e1)
deallocate(rh_e1)
deallocate(ps_e1)
deallocate(ts_e1)
deallocate(t_m2)
deallocate(h_m2)
deallocate(u_m2)
deallocate(v_m2)
deallocate(rh_m2)
deallocate(ps_m2)
deallocate(ts_m2)
deallocate(t_e2)
deallocate(h_e2)
deallocate(u_e2)
deallocate(v_e2)
deallocate(rh_e2)
deallocate(ps_e2)
deallocate(ts_e2)
deallocate(t_m3)
deallocate(h_m3)
deallocate(u_m3)
deallocate(v_m3)
deallocate(rh_m3)
deallocate(ps_m3)
deallocate(ts_m3)
deallocate(t_e3)
deallocate(h_e3)
deallocate(u_e3)
deallocate(v_e3)
deallocate(rh_e3)
deallocate(ps_e3)
deallocate(ts_e3)
deallocate(t_m4)
deallocate(h_m4)
deallocate(u_m4)
deallocate(v_m4)
deallocate(rh_m4)
deallocate(ps_m4)
deallocate(ts_m4)
END PROGRAM Total
subroutine write_FNL_4d(fid,fname,t,h,u,v,rh,ps,ts,nxx,nyy,nzz,nt,ncep_up)
implicit none
character(len=100),intent(in) :: fname
integer,intent(in) :: fid,nxx,nyy,nzz,nt
logical,intent(in) :: ncep_up
REAL,dimension(nxx,nyy,nzz,nt),intent(in) :: h,t,u,v,rh
REAL,dimension(nxx,nyy,nt),intent(in) :: ps,ts
integer :: i,j,k,nrec,kk,ks,ke,kinc
open(fid,file=trim(fname),access='direct',form='unformatted',status='unknown',recl=nxx*nyy*4)
nrec=1
do kk=1,nt
IF (ncep_up) THEN
ks = 1
ke = nzz
kinc = 1
ELSE
ks = nzz
ke = 1
kinc = -1
END IF
do k=ks,ke,kinc
write(fid,rec=nrec) ((t(i,j,k,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
do k=ks,ke,kinc
write(fid,rec=nrec) ((h(i,j,k,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
do k=ks,ke,kinc
write(fid,rec=nrec) ((u(i,j,k,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
do k=ks,ke,kinc
write(fid,rec=nrec) ((v(i,j,k,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
IF (ncep_up) THEN
ks = 1
ke = nzz-5
kinc = 1
ELSE
ks = nzz-5
ke = 1
kinc = -1
END IF
do k=ks,ke,kinc
write(fid,rec=nrec) ((rh(i,j,k,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
write(fid,rec=nrec) ((ps(i,j,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
write(fid,rec=nrec) ((ts(i,j,kk),i=1,nxx),j=1,nyy)
nrec=nrec+1
end do
close(fid)
end subroutine
subroutine read_FNL_4d(fid,fname,t,h,u,v,rh,ps,ts,nxx,nyy,nzz,nt,ncep_up)
implicit none
character(len=100),intent(in) :: fname
integer,intent(in) :: fid,nxx,nyy,nzz,nt
logical,intent(in) :: ncep_up
REAL,dimension(nxx,nyy,nzz),intent(out) :: h,t,u,v,rh
REAL,dimension(nxx,nyy),intent(out) :: ps,ts
integer :: i,j,k,nrec,kk,ks,ke,kinc
open(fid,file=trim(fname),access='direct',form='unformatted',status='unknown',recl=nxx*nyy*4)
nrec=1
do kk=1,nt
IF (ncep_up) THEN
ks = 1
ke = nzz
kinc = 1
ELSE
ks = nzz
ke = 1
kinc = -1
END IF
do k=ks,ke,kinc
read(fid,rec=nrec) ((t(i,j,k),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
do k=ks,ke,kinc
read(fid,rec=nrec) ((h(i,j,k),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
do k=ks,ke,kinc
read(fid,rec=nrec) ((u(i,j,k),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
do k=ks,ke,kinc
read(fid,rec=nrec) ((v(i,j,k),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
IF (ncep_up) THEN
ks = 1
ke = nzz-5
kinc = 1
ELSE
ks = nzz-5
ke = 1
kinc = -1
END IF
do k=ks,ke,kinc
read(fid,rec=nrec) ((rh(i,j,k),i=1,nxx),j=1,nyy)
nrec=nrec+1
enddo
read(fid,rec=nrec) ((ps(i,j),i=1,nxx),j=1,nyy)
nrec=nrec+1
read(fid,rec=nrec) ((ts(i,j),i=1,nxx),j=1,nyy)
nrec=nrec+1
end do
close(fid)
end subroutine
|