（求助）关于湿q矢量，在线等

shenqiankai · 发表于 2017-3-14 16:54:00

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

在论坛里面下载了http://bbs.06climate.com/forum.php?mod=viewthread&tid=44512的湿q矢量的脚本来改的，相应的公式也找到了，但是算出来的位数不对，现在把相应的gs贴出来，请高手帮忙看一下哪里出错了。
'reinit'
'open D:\sq\pm\170222\fnl_20170221_12_00.ctl'

'set mpdset cnworld'
'set map 1 1 4'
tt=1
while(tt<25)
'set t 'tt''
'set lon 40 170'
'set lat 0 90'
'set lev 700'
'set xlopts 1 5 0.16'
'set ylopts 1 5 0.16'
'set xlint 2'
'set ylint 2'
'set grads off'
'set grid off'
*
'define dx=2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0'
'define dy=2.0*6370949.0*3.14159/180.0'
'define f=2*7.292e-5*sin(lat*3.1415/180)'
'define r=287'
'define a=17.1543'
'define b=36'
'define t=tmpprs'
'define rh=rhprs'
'define prs=lev'
if(t>263.0);
'define a1=0.622*6.11*EXP(17.26*(t-273.16)/(t-35.86))'
'define b1=prs-0.378*6.11*EXP(17.26*(t-273.16)/(t-35.86))'
'define  qs=a1/b1'
endif
if(t<=263.0);
'define a1=0.622*6.11*EXP(21.87*(TMPprs-273.16)/(TMPprs-7.66))'
'define b1=lev-0.378*6.11*EXP(21.87*(TMPprs-273.16)/(TMPprs-7.66))'
'define  qs=a1/b1'
endif
'define  q=qs*rh/100'

'define dux=cdiff(ugrdprs,x)'
'define dvx=cdiff(vgrdprs,x)'
'define du=ugrdprs(z-1)-ugrdprs(z+1)'
'define dv=vgrdprs(z-1)-vgrdprs(z+1)'
'define dp=100*(lev(z-1)-lev(z+1))'
'define qx1=(f*((dv/dp)*(dux/dx)-(du/dp)*(dvx/dx)))/2'
*
'define p=lev'
'define cp=(1+0.86*q)*1004'
'define pt=t*pow((1000/prs),r/cp)'
'define ht=r*pow(p/1000,r/cp)/p'
'define dux=cdiff(ugrdprs,x)'
'define dvx=cdiff(vgrdprs,x)'
'define dptx=cdiff(pt,x)'
'define dpty=cdiff(pt,y)'
'define qx2=-ht*((dux/dx)*(dptx/dx)+(dvx/dx)*(dpty/dy))/2'

*
'define p=lev'
'define cp=(1+0.86*q)*1004'
'define l=2.5*1e6-2323*(t-273.16)'
'define tv=t*(1+0.61*qs)'
'define c=a*(273.16-b)/pow((t-b),2))'
'define qp=(c*r*tv-cp)*qs/((c*l*qs-cp)*p)'
'define lhe=l*vvelprs*qp'
'define dlhex=cdiff(lhe,x)'
'define qx3=-(dlhex/dx)*r/(cp*p*2)'
'define qx=qx1+qx2+qx3'
*
'define duy=cdiff(ugrdprs,y)'
'define dvy=cdiff(vgrdprs,y)'
'define du=ugrdprs(z-1)-ugrdprs(z+1)'
'define dv=vgrdprs(z-1)-vgrdprs(z+1)'
'define dp=100*(lev(z-1)-lev(z+1))'
'define qy1=(f*((dv/dp)*(duy/dy)-(du/dp)*(dvy/dy)))/2'
*
'define p=lev'
'define cp=(1+0.86*q)*1004'
'define ht=r*pow(p/1000,r/cp)/p'
'define duy=cdiff(ugrdprs,y)'
'define dvy=cdiff(vgrdprs,y)'
'define dptx=cdiff(pt,x)'
'define dpty=cdiff(pt,y)'
'define qy2=-ht*((duy/dy)*(dptx/dx)+(dvy/dy)*(dpty/dy))/2'
*
'define p=lev'
'define cp=(1+0.86*q)*1004'
'define l=2.5*1e6-2323*(t-273.16)'
'define tv=t*(1+0.61*qs)'
'define c=a*(273.16-b)/pow((t-b),2))'
'define qp=(c*r*tv-cp)*qs/((c*l*qs-cp)*p)'
'define lhe=l*vvelprs*qp'
'define dlhey=cdiff(lhe,y)'
'define qy3=-(dlhey/dy)*r/(cp*p*2)'
'define qy=qy1+qy2+qy3'
'd mag(qx,qy)*1e12'

哈库拉玛塔塔000 · 发表于 2017-4-27 14:03:16

楼主解决了吗？？？
if(t>263.0);
'define a1=0.622*6.11*EXP(17.26*(t-273.16)/(t-35.86))'
'define b1=prs-0.378*6.11*EXP(17.26*(t-273.16)/(t-35.86))'
'define qs=a1/b1'
endif
if(t<=263.0);
'define a1=0.622*6.11*EXP(21.87*(TMPprs-273.16)/(TMPprs-7.66))'
'define b1=lev-0.378*6.11*EXP(21.87*(TMPprs-273.16)/(TMPprs-7.66))'
'define qs=a1/b1'
endif

还有楼主程序中这部分，为什么要这样写呢，我看到源程序中没有这部分啊？？？

磨人的小妖精儿 · 发表于 2017-5-6 11:16:44

里面变成加cp试试

磨人的小妖精儿 · 发表于 2017-5-6 11:21:17

你下载这个计算qs有些问题

磨人的小妖精儿 · 发表于 2017-5-6 11:21:54

C$DEBUG
c$LARGE
C    CALCULATE  gaimoist-agestrophic Q VECT AND DIVERGENCE OF Q VECT

   parameter(m=220,n=218,kx=19,tx=49,d=15000)

   REAL IP(KX),sqx(M,N,KX,tx),sqy(M,N,KX,tx),sqd(M,N,KX,tx)
   real u(M,N,KX,tx),ux(M,N,KX,tx),uy(m,n,kx,tx),up(m,n,kx,tx)
   real v(M,N,KX,tx),vx(M,N,KX,tx),vy(m,n,kx,tx),vp(m,n,kx,tx)
   real ct(M,N,KX,tx),ctx(M,N,KX,tx),cty(m,n,kx,tx),ctp(m,n,kx,tx)
   real hh(KX),W(M,N,KX,tx),AA(M,N,KX,tx)
   real td(M,N,KX,tx),t(M,N,KX,tx),ee(m,n,kx,tx),es(m,n,kx,tx)
   real q(M,N,KX,tx),qs(M,N,KX,tx),bs(M,N,KX,tx),bb(m,n,kx,tx)
   real sqxx(M,N,KX,tx),sqyy(M,N,KX,tx),tv(m,n,kx,tx)
   real c(m,n,kx,tx),ttx(m,n,kx,tx),tty(m,n,kx,tx)
   real aax(M,N,KX,tx),aay(m,n,kx,tx)
real fqx(m,n,kx),fqy(m,n,kx),tpp(m,n,kx)

   character    c2*2
      character    dh(tx)*2

   DATA dh/'01','02','03','04','05','06','07','08',
   * '09','10','11','12','13','14', '15','16','17','18','19',
   * '20','21','22','23','24','25','26','27','28','29',
   * '30','31','32','33','34','35','36','37','38','39',
   * '40','41','42','43','44','45','46','47','48','49'/

   DATA IP/ 1000.00,950.00,900.00,850.00,
   &    800.00,750.00,700.00,650.00,600.00,550.00,500.00,450.00,
   &    400.00,350.00,300.00,250.00,200.00,150.00,100.00/

real a,b,CP,L,f
   parameter(a=17.1543, b=36.,R=287.05,CP=1004.)
   parameter(L=2500000., f=7.29e-5)
open(14,file='d:\1\t.grd',form='binary')
read(14)((((t(i,j,k,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
close(14)
open(1,file='d:\1\td.grd',form='binary')
read(1)((((td(i,j,k,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
close(1)

   OPEN(2,FILE='d:\1\u.grd',form='binary')

   READ(2)((((u(I,J,K,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
   close(2)

   OPEN(32,FILE='d:\1\v.grd',form='binary')
   READ(32)((((v(I,J,K,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
   close(32)

   OPEN(30,FILE='d:\1\w.grd',form='binary')
   READ(30)((((w(I,J,K,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
c do 28  it=1,tx
c do 28 k=1,kx
c    do 28 j=1,n
c    do 28 i=1,m
c if(w(i,j,k,it).ne.1.e35)then
c w(i,j,k,it)=-1.293*10*w(i,j,k,it)
c else
c w(i,j,k,it)=0.0
c endif
c28    continue
   close(30)
c       do 25 it=1,tx
c    do 25 j=1,n-1
C    do 25 i=1,m-1
C    do 25 k=1,kx

C if(t(i,j,k,tx).eq.1.e35.or.t(i,j+1,k,tx).eq.1.e35.or.
C    $t(i+1,j,k,tx).eq.1.e35.or.t(i+1,j+1,k,tx).eq.1.e35)then
C CT(i,j,k,tx)=1.e35
C    else
C    T(i,j,k,tx)=(t(i,j,k,tx)+t(i,j+1,k,tx)+t(i+1,j,k,tx)
C & +t(i+1,j+1,k,tx))/4
C    CT(I,J,K,tx)=T(I,J,K,tx)*(1000./IP(K))**0.286
C endif
C25    continue
   do 25 it=1,tx
   do 25 j=1,n
   do 25 i=1,m
   do 25 k=1,kx

if(t(i,j,k,it).eq.1.e35)then
CT(i,j,k,it)=1.e35
   else
   CT(I,J,K,it)=T(I,J,K,it)*(1000./IP(K))**0.286
endif
25    continue

   call jsp(ct,ctp)

   do 29 k=1,kx
29 hH(K)=R/IP(K)*(IP(K)/1000.)**0.286
do 30 it=1,tx
   do 30 k=1,kx
   DO 30 j=1,n
   DO 30 i=1,m
if(u(i,j,k,it).eq.1.e35.or.v(i,j,k,it).eq.1.e35.or.
   $ w(i,j,k,it).eq.1.e35.or.td(i,j,k,it).eq.1.e35.or.
   $ t(i,j,k,it).eq.1.e35)then
   aa(i,j,k,it)=1.e35
endif

if(u(i,j,k,it).ne.1.e35.and.v(i,j,k,it).ne.1.e35.and.
   $ w(i,j,k,it).ne.1.e35.and.td(i,j,k,it).ne.1.e35.and.
   $  t(i,j,k,it).ne.1.e35)then
   ee(i,j,k,it)=6.11*exp(a*(Td(I,J,K,it)-273.16)/(Td(I,J,K,it)-b))
   es(i,j,k,it)=6.11*exp(a*(T(I,J,K,it)-273.16)/(T(I,J,K,it)-b))

   Q(I,J,K,it)=0.622*ee(I,J,K,it)/(IP(K)-0.378*ee(i,j,k,it))
   QS(I,J,K,it)=0.622*es(I,J,K,it)/(IP(K)-0.378*es(i,j,k,it))
   bs(i,j,k,it)=ee(i,j,k,it)/es(i,j,k,it)

   C(I,J,K,it)=a*(273.16-b)/(T(I,J,K,it)-b)**2

   TV(I,J,K,it)=T(I,J,K,it)*(1+0.61*QS(I,J,K,it))

   BB(I,J,K,it)=(C(I,J,K,it)*R*TV(I,J,K,it)-CP)
   $ *QS(I,J,K,it)/((C(I,J,K,it)*L*QS(I,J,K,it)+CP)*ip(k))

   if(ctp(i,j,k,it).lt.0.and.bs(i,j,k,it).ge.0.8.and.
   & w(i,j,k,it).lt.0)then
   AA(I,J,K,it)=L*R*w(I,J,K,it)/(CP*Ip(K))*BB(I,J,K,it)/100
   end if

   if(ctp(i,j,k,it).ge.0.or.bs(i,j,k,it).le.0.8.or.
   &    w(i,j,k,it).ge.0)then
   aa(i,j,k,it)=1.e35
   end if
endif
30 CONTINUE
C    print *,cc(100,50,5)
C pause
   call jsp(u,up)
call jsp(v,vp)
   call jsx(aa,aax)
   call jsy(aa,aay)
   call jsx(ct,ctx)
   call jsy(ct,cty)
   call jsx(u,ux)
   call jsy(u,uy)
   call jsx(v,vx)
   call jsy(v,vy)
call jsx(t,ttx)
   call jsy(t,tty)

do 35 it=1,tx
   do 35 k=1,kx
   do 35 j=1,n
   do 35 i=1,m
sqx(i,j,k,it)=0.5*(f*(vp(i,j,k,it)*ux(i,j,k,it)-
   $up(i,j,k,it)*vx(i,j,k,it))-hh(k)*(ux(i,j,k,it)*ctx(i,j,k,it)
   $+vx(i,j,k,it)*cty(i,j,k,it)))-0.5*aax(i,j,k,it)
   sqy(i,j,k,it)=0.5*(f*(vp(i,j,k,it)*uy(i,j,k,it)-
   $up(i,j,k,it)*vy(i,j,k,it))-hh(k)*(uy(i,j,k,it)*ctx(i,j,k,it)
   $+vy(i,j,k,it)*cty(i,j,k,it)))-0.5*aay(i,j,k,it)
c    fs(i,j,k,it)=sqx(i,j,k,it)*ttx(i,j,k,it)+sqy(i,j,k,it)
c    $*tty(i,j,k,it)

35 continue
   call jsx(sqx,sqxx)
   call jsy(sqy,sqyy)

do 37 it=1,tx
   do 37 k=1,kx
   do 37 j=1,n
   do 37 i=1,m
c    fs(i,j,k,it)=fs(i,j,k,it)*1.0E16

c    sqx(i,j,k,it)=sqx(i,j,k,it)
c    sqy(i,j,k,it)=sqy(i,j,k,it)
c    sqxx(i,j,k,it)=sqxx(i,j,k,it)
c    sqyy(i,j,k,it)=sqyy(i,j,k,it)
   sqd(i,j,k,it)=2*(sqxx(i,j,k,it)+sqyy(i,j,k,it))

37 continue
c    print *, sqx
C    THE UNITE OF sQX AND sQY IS 'E-10m/(hpa.S**3)'
C    THE UNITE OF sQD IS 'E-14/(HPA.S**3)'
ccccccccccccccccccccccccccccccccccccccccccccccc
c    print *,ct
      do tt=1,tx
   open(11,file='d:\1\'//dh(tt)//'_fqx.grd',
   &form='binary')
   do 11 K=1,kx
   do 11 j=1,n
   do 11 i=1,m
fqx(I,J,K)=sqx(I,J,K,tt)
   write(11) fqx(I,J,K)
11    CONTINUE
   CLOSE(11)

   open(12,file='d:\1\'//dh(tt)//'_fqy.grd',
   &form='binary')
   do 12 K=1,kx
   do 12 j=1,n
   do 12 i=1,m
   fqy(I,J,K)=sqy(I,J,K,tt)
   write(12) fqy(I,J,K)
12    CONTINUE
   CLOSE(12)

   open(13,file='d:\1\'//dh(tt)//'_tc.grd',
   &form='binary')
   do 13 K=1,kx
   do 13 j=1,n
   do 13 i=1,m
      tpp(I,J,K)=ct(I,J,K,tt)
   write(13) tpp(I,J,K)
13    CONTINUE
   CLOSE(13)
   enddo

cccccccccccccccccccccccccccccccccccccccccccccccccccccc
   OPEN(20,FILE='d:\1\sqd.dat',form='binary')
c OPEN(21,FILE='d:\22\sq.dat',form='binary')
c OPEN(22,FILE='d:\22\fs.dat',form='binary')

c    write(14) ip(k)
   WRITE(20)((((sqd(I,J,K,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
c    WRITE(22)((((fs(I,J,K,it),i=1,m),j=1,n),k=1,kx),it=1,tx)
c close(22)
c    WRITE(21)((((sqx(I,J,K,tx),i=1,m),j=1,n),k=1,kx),it=1,tx)
c write(21)((((sqy(i,j,k,tx),i=1,m),j=1,n),k=1,kx),it=1,tx)
close(20)
c close(21)
c3000 format(f8.2)
c4000 FORMAT(41F8.2)


   END

      subroutine  jsx(ss,ssx)
   parameter(m=220,n=218,kx=19,tx=49)
   DIMENSION ss(M,N,KX,tx),ssx(m,n,kx,tx)

   D=8000
do 101 it=1,49
   do 100 k=1,kx
   do 100 j=1,n
   ssx(1,j,k,it)=(ss(2,j,k,it)-ss(1,j,k,it))/d
   ssx(m,j,k,it)=(ss(m,j,k,it)-ss(m-1,j,k,it))/d
100 continue
   do 101 k=1,kx
   do 101 j=1,n
   do 101 i=2,m-1
if(ss(i+1,j,k,it).eq.1.e35.or.ss(i-1,j,k,it).eq.1.e35)then
ssx(i,j,k,it)=0.0
else
   ssx(i,j,k,it)=(ss(i+1,j,k,it)-ss(i-1,j,k,it))/(2*d)
   endif
101 continue
   return
   end

      subroutine  jsy(ss,ssy)
   parameter(m=220,n=218,kx=19,tx=49)
      DIMENSION ss(M,N,KX,tx),ssy(m,n,kx,tx)

   D=8000
do 103 it=1,49
   do 102 k=1,kx
   do 102 i=1,m
   ssy(i,1,k,it)=(ss(i,2,k,it)-ss(i,1,k,it))/d
   ssy(i,n,k,it)=(ss(i,n,k,it)-ss(i,n-1,k,it))/d
102 continue
   do 103 k=1,kx
   do 103 i=1,m
   do 103 j=2,n-1
if(ss(i,j+1,k,it).eq.1.e35.or.ss(i,j-1,k,it).eq.1.e35)then
ssy(i,j,k,it)=0.0
else
   ssy(i,j,k,it)=(ss(i,j+1,k,it)-ss(i,j-1,k,it))/(2*d)
   endif
103 continue
   return
   end

      subroutine  jsp(ss,ssy)
   parameter(m=220,n=218,kx=19,tx=49)
   DIMENSION ss(M,N,KX,tx), ssy(m,n,kx,tx)
      real ip(kx)
   DATA IP/1000.00,950.00,900.00,850.00,
   & 800.00,750.00,700.00,650.00,600.00,550.00,500.00,450.00,
   &400.00,350.00,300.00,250.00,200.00,150.00,100.00/

   D=8000
do 105 it=1,49
   do 13 i=1,m
   do 13 j=1,n
if(ss(i,j,1,it).eq.1.e35.or.ss(i,j,2,it).eq.1.e35)then
ssy(i,j,1,it)=0.0
else
   ssy(i,j,1,it)=(ss(i,j,2,it)-ss(i,j,1,it))/(ip(2)-ip(1))
endif
13 continue
   do 14 i=1,m
   do 14 j=1,n
if(ss(i,j,19,it).eq.1.e35.or.ss(i,j,18,it).eq.1.e35)then
   ssy(i,j,19,it)=0.0
else
   ssy(i,j,19,it)=(ss(i,j,19,it)-ss(i,j,18,it))/(ip(19)-ip(18))
endif
14 continue

   do 15 j=1,n
   do 15 i=1,m
do 15 k=2,kx-1

if(ss(i,j,k+1,it).eq.1.e35.or.ss(i,j,k-1,it).eq.1.e35)then
ssy(i,j,k,it)=0.0
else
   ssy(i,j,k,it)=(ss(i,j,k+1,it)-ss(i,j,k-1,it))/(ip(k+1)-ip(k-1))
   endif
15 continue
105 continue
   return
   end

shenqiankai · 发表于 2017-5-17 16:32:23

磨人的小妖精儿发表于 2017-5-6 11:21
C$DEBUG
c$LARGE
C CALCULATE gaimoist-agestrophic Q VECT AND DIVERGENCE OF Q VECT

你这个是fortran程序吗

磨人的小妖精儿 · 发表于 2017-5-29 09:50:00

可以自己编matlab

泡沫学天气 · 发表于 2017-6-5 09:45:41

请问楼主最后问题解决了么

泡沫学天气 · 发表于 2017-6-5 09:46:01

请问楼主最后问题解决了么

泡沫学天气 · 发表于 2017-6-5 09:48:48

磨人的小妖精儿发表于 2017-5-6 11:21
C$DEBUG
c$LARGE
C CALCULATE gaimoist-agestrophic Q VECT AND DIVERGENCE OF Q VECT

用grads画 fnl资料请问您有gs文件么？或者这个程序的说明感激不尽

		自动登录	找回密码
密码			立即注册

[分享资料] （求助）关于湿q矢量，在线等

登录后查看更多精彩内容~