任意两点间斜剖面做法总结【利用GRADS自带函数】

river

深度狙杀 发表于 2014-8-22 10:18
我现在想按照一篇文献里的做法，画沿着一条直线的剖面。
为了直观（如图2），想用‘set y 1 12’但那个c ...

不行，那个函数现在还不支持时间的剖面，所以你只能通过计算这条直线上的值，输出成二进制文件，编写ctl再画图

summer1217

很牛的帖子~

深度狙杀

river 发表于 2014-8-22 19:09
不行，那个函数现在还不支持时间的剖面，所以你只能通过计算这条直线上的值，输出成二进制文件，编写ctl ...

我的资料里时间有24个时次，但层次只有一层，而这个函数支持的是层次，我就改动了一下ctl吧时次当层次处理，应该不会影响结果

我想把时间轴调到横坐标。但用了 xyrev on 之后就画不了了 ，秒退了

river

深度狙杀 发表于 2014-8-23 14:26
我的资料里时间有24个时次，但层次只有一层，而这个函数支持的是层次，我就改动了一下ctl吧时次当层次处 ...

首先我很佩服你，你的这个想法非常好，是可以的，不会影响最终结果。其次批评你一下，出现问题了就需要详细描述一下，顺便把你的gs帖上来，不然你这么说我也不知道问题在哪。最后再给你点个赞

深度狙杀

river 发表于 2014-8-23 16:28
首先我很佩服你，你的这个想法非常好，是可以的，不会影响最终结果。其次批评你一下，出现问题了就需要详 ...

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深的平均

请你帮我支支招，谢谢

river

深度狙杀 发表于 2014-8-24 15:05
ctl 将时间维想象为垂直坐标（24个时次）
dset D:\Mechanism\Basic\diagnose\result\Zadv\za.grd
undef ...

我晕，这么专业的问题，你最好还是发个新帖来提问，这样做过的人就可以给你出主意了。你的set x 1 14是什么意思呢？

深度狙杀

river 发表于 2014-8-24 17:18
我晕，这么专业的问题，你最好还是发个新帖来提问，这样做过的人就可以给你出主意了。你的set x 1 14是什 ...

set x 这个地方我觉得只要是有变化的就可以，没有实际意义。我set x 1 2或者set x 1 10等都一样的，只要不固定就行。主要是为了实现coll2gr这个函数x与z两个维同时变化的目的，实际的坐标由xlab标出即可。

关于值偏大的问题我再阅读几篇文献，看看到底有没有错。还不行就发贴
谢谢你的指导

river

深度狙杀 发表于 2014-8-24 19:02
set x 这个地方我觉得只要是有变化的就可以，没有实际意义。我set x 1 2或者set x 1 10等都一样的，只要 ...

所以这是个问题，我的贴子里set z 1 12是和前面对应的，并不是没有意义的······

小预报员

楼主真好！

深度狙杀

river 发表于 2014-8-24 22:12
所以这是个问题，我的贴子里set z 1 12是和前面对应的，并不是没有意义的······

嗯 我再好好看看  谢谢