登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 lqouc 于 2014-5-29 08:48 编辑
老师您好，下面这段程序是南方气旋物理量诊断的，缺测值已经赋好了，数据对比过源文件了，应该没问题。就是计算散度，涡度时，量级不对
想请问一下：
我想到可能由于下面问题导致的：
（1）例如求地转风涡度时，格距应该用delta=4.0还是化成弧度，用sign=4.0*pi/180。vorg(i,j,iz,it)=9.8*((height(i+1,j,iz,it)+height(i-1,j,iz,it)-2*height(i,j,iz,it))/((cos(lat(j))*sign)**2)+(height(i+1,j,iz,it)+height(i-1,j,iz,it)-2*height(i,j,iz,it))/sign**2-(height(i,j+1,iz,it)-height(i,j-1,iz,it))/(2*sign*tan(lat(j))))/(2*Omega*sin(lat(j))*R**2)
(2)边界的缺测值已经赋上，grads里面undef也定义了相同的值，可是出图还是很奇怪
可是我处理了还是不行，能麻烦老师帮忙看一下吗？谢谢老师
下附fortran和grads的ctl文件
! read 20090416/high/height;temper;t-td;uv
! lon 32~160 interval: 4 degree
! lat 12~80 interval: 4 degree
! time 1~12: 2009.04.16.20:00~04.22.08:00 interval: 12 hours
! U,v,height,t: 1000,925,850,700,500,400,300,250,200,150,100 hpa
REAL,PARAMETER :: Omega=7.292e-5,R=6371,PI=3.1415926,Delta=4,DeltaT=12
INTEGER,PARAMETER :: nx=33,ny=18,nz=11,nt=12
REAL F,EE,AA,BB,sign
INTEGER WtTopLevel(NX,NY,NT)
REAl P(NZ),deltP(10),lat(ny),lon(nx),sigmadeltD
CHARACTER timefile(12)*12 ,levelfile(11)*4
real,allocatable :: temper(:,:,:,:),q(:,:,:,:),u(:,:,:,:), &
v(:,:,:,:),height(:,:,:,:),w(:,:,:,:), &
qu(:,:,:,:),qv(:,:,:,:),adq(:,:,:,:), &
adqv(:,:,:,:),theta(:,:,:,:),Wt(:,:,:,:), &
vorg(:,:,:,:),voro(:,:,:,:),div(:,:,:,:), &
deltD(:,:,:),ttd(:,:,:,:),thetaP(:,:,:,:)
allocate (temper(NX,NY,NZ,NT),q(NX,NY,nz,NT),u(NX,NY,NZ,NT), &
v(NX,NY,NZ,NT),height(NX,NY,NZ,NT),w(NX,NY,NZ,NT), &
qu(NX,NY,NZ,NT),qv(NX,NY,NZ,NT),adq(NX,NY,NZ,NT), &
adqv(NX,NY,NZ,NT),theta(NX,NY,NZ,NT),Wt(NX,NY,NZ,NT), &
vorg(NX,NY,NZ,NT),voro(NX,NY,NZ,NT),div(NX,NY,NZ,NT), &
deltD(NX,NY,NT),ttd(NX,NY,NZ,NT),thetaP(NX,NY,NZ,NT))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DATA P /1000,925,850,700,500,400,300,250,200,150,100/
DATA deltP/75, 75, 50, 200,100,100,50, 50, 50, 50/
DATA timefile/'09041620.000','09041708.000','09041720.000','09041808.000','09041820.000','09041908.000', &
'09041920.000','09042008.000','09042020.000','09042108.000','09042120.000','09042208.000'/
DATA levelfile/'1000','925','850','700','500','400','300','250','200','150','100'/
do i=1,nx
lon(i)=32+(i-1)*Delta
lon(i)=lon(i)*pi/180.
enddo
do i=1,ny
lat(i)=12+(i-1)*Delta
lat(i)=lat(i)*pi/180.
enddo
sign=4.0*pi/180
!--------------read temperature----------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='f:/micaps/temper/'//trim(levelfile(iz))//'/'//timefile(it))
do i=1,4
read(11,*)
enddo
do j=NY,1,-1
read(11,*)(temper(i,j,iz,it),i=1,NX) ! C
enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read u,v -------------------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='f:/micaps/uv/'//trim(levelfile(iz))//'/'//timefile(it))
do i=1,3
read(11,*)
enddo
do j=36,1,-1
if(j>NY) then
read(11,*)(u(i,j-ny,iz,it),i=1,NX)
else
read(11,*)(v(i,j,iz,it),i=1,NX)
endif
enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read t-td -------------------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='f:/micaps/t-td/'//trim(levelfile(iz))//'/'//timefile(it))
do i=1,4
read(11,*)
enddo
do j=NY,1,-1
read(11,*)(ttd(i,j,iz,it),i=1,NX) ! C
enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read height -------------------------
DO IZ=1,NZ ! level
DO IT=1,NT ! time
OPEN(11,file='f:/micaps/height/'//trim(levelfile(iz))//'/'//timefile(it))
do i=1,4
read(11,*)
enddo
do j=NY,1,-1
read(11,*)(height(i,j,iz,it),i=1,NX) ! C
enddo
CLOSE(11)
ENDDO
ENDDO
!@@@@@@@@@@@@@@@@@@@@ read end @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!======== t,td => specific humidity 比湿:
DO IZ=1,NZ ! level
DO IT=1,NT ! time
do i=1,nx
do j=1,ny
f=temper(i,j,iz,it)-ttd(i,j,iz,it)
ee=6.1078*exp(17.3*f/(273.16+f-35.86))
q(i,j,iz,it)=622*ee/(p(iz)-0.378*ee)
enddo
enddo
enddo
enddo
!######################## geostrophic wind vorticity 地转风涡度###########################
DO IZ=1,NZ ! level
DO IT=1,NT ! time
do i=1,nx
do j=1,ny
vorg(i,j,iz,it)=-9.99e8
enddo
enddo
do i=2,nx-1
do j=2,ny-1
vorg(i,j,iz,it)=9.8*((height(i+1,j,iz,it)+height(i-1,j,iz,it)-2*height(i,j,iz,it))/((cos(lat(j))*sign)**2)+&
(height(i+1,j,iz,it)+height(i-1,j,iz,it)-2*height(i,j,iz,it))/sign**2- &
(height(i,j+1,iz,it)-height(i,j-1,iz,it))/(2*sign*tan(lat(j))))/(2*Omega*sin(lat(j))*R**2)
enddo
enddo
enddo
enddo
!######################## vorticity of observed wind 实测风涡度 ########################
DO IZ=1,NZ ! level
DO IT=1,NT ! time
do i=1,nx
do j=1,ny
voro(i,j,iz,it)=-9.99e8
enddo
enddo
do i=2,nx-1
do j=2,ny-1
voro(i,j,iz,it)=((v(i+1,j,iz,it)-v(i-1,j,iz,it))/(cos(lat(j))*sign)- &
(u(i,j+1,iz,it)-u(i,j-1,iz,it))/sign+2*u(i,j,iz,it)*tan(lat(j)))/(2*R)
enddo
enddo
enddo
enddo
!######################## divergence 散度 ###########################
DO IZ=1,NZ ! level
DO IT=1,NT ! time
do i=1,nx
do j=1,ny
div(i,j,iz,it)=-9.99e8
enddo
enddo
do i=2,nx-1
do j=2,ny-1
div(i,j,iz,it)=((u(i+1,j,iz,it)-u(i-1,j,iz,it))/(cos(lat(j))*sign)+ &
(v(i,j+1,iz,it)-v(i,j-1,iz,it))/sign-2*v(i,j,iz,it)*tan(lat(j)))/(2*R)
enddo
enddo
enddo
enddo
print*,div(2,2,1,1)
!######################## vertical velocity 垂直速度 ###########################
!-------- correct w w修正 --------------
!- - - - - thermaldynamic method for w at the top level - - - - - - -
!- - - - - 100,150hpa: reference to 100,150,200 hpa
!############### vapor flux & vapor flux divergence 水汽通量和水汽通量散度 ##################
DO IZ=1,NZ ! level
DO IT=1,NT ! time
do i=1,nx
do j=1,ny
adqv(i,j,iz,it)=-9.99e8
enddo
enddo
do i=2,nx-1
do j=2,ny-1
qu(i,j,iz,it)=u(i,j,iz,it)*q(i,j,iz,it)/9.8
qv(i,j,iz,it)=v(i,j,iz,it)*q(i,j,iz,it)/9.8
AA=u(i,j,iz,it)*(q(i+1,j,iz,it)-q(i-1,j,iz,it))/(cos(lat(j))*sign)+v(i,j,iz,it)*(q(i,j+1,iz,it)-q(i,j-1,iz,it))/sign
adqv(i,j,iz,it)=(AA/(2*R)+q(i,j,iz,it)*div(i,j,iz,it))/9.8
enddo
enddo
enddo
enddo
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
OPEN(8,FILE='f:/out.grd',form='binary')
do it=1,nt
WRITE(8) (((u(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((v(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((q(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((temper(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((height(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((voro(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((vorg(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((div(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((qu(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((qv(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
WRITE(8) (((adqv(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
enddo
CLOSE(8)
do it=1,nt
OPEN(9,FILE='f:/div2.dat')
WRITE(9,*) (((div(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
close(9)
enddo
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
deallocate (temper,q,u,v,height,w,qu,qv,adq,adqv,theta,Wt, &
vorg,voro,div,deltD,ttd,thetaP)
END
ctl文件
dset f:\out.grd
undef -9.9900000E+08
xdef 33 linear 32 4
ydef 18 linear 12 4
zdef 11 levels 1000 925 850 700 500 400 300 250 200 150 100
tdef 12 linear 20:00Z16jul2009 12hr
vars 9
u 11 99 u wind
v 11 99 v wind
q 11 99 specific humidity
qu 11 99 qu wind
qv 11 99 qv wind
div 11 99 divergence
adqv 11 99 q divergence
voro 11 99 vorticity of observed wind
vorg 11 99 geostrophic wind vorticity
endvars
|
窥视卡
雷达卡
发表于 2014-5-23 20:58:36
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
