- 积分
- 746
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-10-7
- 最后登录
- 1970-1-1
|
发表于 2014-8-24 15:05:34
|
显示全部楼层
本帖最后由 深度狙杀 于 2014-8-24 15:13 编辑
ctl 将时间维想象为垂直坐标(24个时次)
dset D:\Mechanism\Basic\diagnose\result\Zadv\za.grd
undef 32767.
title anormal sst
xdef 360 linear -179. 1.0
ydef 416 linear -73.834 0.333
zdef 24 linear 2 10
tdef 1 linear 01jan1990 1mo
vars 1
zad 24 -999 sst
endvars
编写gs
'enable print D:\Mechanism\Basic\diagnose\result\Zadv\zad.gmf '
'set grads off'
'set csmooth on'
'set x 1'
'set y 1'
'set t 1'
'set z 4 16'
lon1 = 165.
lon2 = 245.
lat1 = -5.
lat2 = 25.0
lon = lon1
'collect 1 free'
while (lon <= lon2)
lat = lat1 + (lat2-lat1)*(lon-lon1)/(lon2-lon1)
'collect 1 gr2stn(zad,'lon','lat')'
lon = lon + 1
endwhile
'set x 1 14'
'set ylabs 5S,165E|0N,178E|5N,169W|10N,155W|15N,142W|20N,128W|25N,115W'
'set xlabs 90/04|90/06|90/08|90/10|90/12|91/02|91/04'
'set xyrev on'
'set gxout contour'
'set missconn on'
'set cthick 8'
'set cint 0.3'
'd coll2gr(1,13)'
'print'
'disable print'
;
PS:这个还未解决,又出现了新问题,编出的程序,计算值和论文相比明显偏大。
方法是根据海温诊断方程,计算每一项的贡献。
我的参数定义
module dat
implicit none
integer,parameter::nx=360,ny=416,nyear=2,nt=12
real,parameter::undef=32767.,month=30*24*60*60
real,parameter::pi=3.1416,r=6371000
real,dimension(nx,ny,nt)::u,v,temp,w,tgave
real,dimension(nx,ny,nt,nyear)::ua,va,wa,ta,q,tte,tga
real,dimension(nx,ny,nt,nyear)::zad,mad,vad,hq,total
integer::i,j,t,k
real::dy,sig; real,dimension(ny)::dx
end module
步长:
dy=2*pi*r*(0.333/360)
print*,dy
do j=1,ny
sig=-73.834+0.333*(j-1)
dx(j)=2*pi*r*cosd(sig)/360.
!print*,dx(j)
end do
纬向平流计算过程:边界采用前后差,中间采用中央差
zad=32767.
do k=1,nyear
do t=1,nt
!!!!!!!!!!!left boundary
do j=1,ny
if(u(1,j,t)/=undef.and.ta(2,j,t,k)/=undef.and.ta(1,j,t,k)/=undef.and.ua(1,j,t,k)/=undef.and.temp(2,j,t)/=undef.and.temp(1,j,t)/=undef)then
zad(1,j,t,k)=-u(1,j,t)*(ta(2,j,t,k)-ta(1,j,t,k))/(dx(j)) &
&-ua(1,j,t,k)*(temp(2,j,t)-temp(1,j,t))/(dx(j)) &
&-ua(1,j,t,k)*(ta(2,j,t,k)-ta(1,j,t,k))/(dx(j))
zad(1,j,t,k)=zad(1,j,t,k)*month
! write(*,*) zad(i,j,t,k)
else
zad(1,j,t,k)=undef
endif
!!!!!!!!!!!right boundary
if(u(nx,j,t)/=undef.and.ta(nx-1,j,t,k)/=undef.and.ta(nx,j,t,k)/=undef.and.ua(nx,j,t,k)/=undef.and.temp(nx-1,j,t)/=undef.and.temp(nx,j,t)/=undef)then
zad(nx,j,t,k)=-u(nx,j,t)*(ta(nx,j,t,k)-ta(nx-1,j,t,k))/(dx(j)) &
&-ua(nx,j,t,k)*(temp(nx,j,t)-temp(nx-1,j,t))/(dx(j)) &
&-ua(nx,j,t,k)*(ta(nx,j,t,k)-ta(nx-1,j,t,k))/(dx(j))
zad(nx,j,t,k)=zad(nx,j,t,k)*month
! write(*,*) zad(i,j,t,k)
else
zad(nx,j,t,k)=undef
endif
do i=2,nx-1
if(u(i,j,t)/=undef.and.ta(i+1,j,t,k)/=undef.and.ta(i-1,j,t,k)/=undef.and.ua(i,j,t,k)/=undef.and.temp(i+1,j,t)/=undef.and.temp(i-1,j,t)/=undef)then
zad(i,j,t,k)=-u(i,j,t)*(ta(i+1,j,t,k)-ta(i-1,j,t,k))/(2*dx(j)) &
&-ua(i,j,t,k)*(temp(i+1,j,t)-temp(i-1,j,t))/(2*dx(j)) &
&-ua(i,j,t,k)*(ta(i+1,j,t,k)-ta(i-1,j,t,k))/(2*dx(j))
zad(i,j,t,k)=zad(i,j,t,k)*month
! write(*,*) zad(i,j,t,k)
else
zad(i,j,t,k)=undef
endif
end do
end do
end do
end do
资料里海洋纬向流单位是m/s,所以计算出来后单位为K/S,我乘以month,然后得到K/month.感觉计算过程没错,就是最后整体偏大。郁闷!!!资料取自GODAS,按前人的做法,选取0-50m深的平均
请你帮我支支招,谢谢
|
-
采用xyrev后提示错误
-
有错误但可以画出图形,不能保存为gmf,且纵坐标只显示几个
-
不采用用xyrev可绘出图形
-
海温诊断方程
|