求fortran 计算整层水汽通量的程序

默戏

同求啊，好纠结

刘晶

grads就能算了

1234abc

同求同求同求

茴糖惹

grads就能算了

Alexander

哪位大佬写一下啊求助啊啊啊{:eb302:}{:eb302:}{:eb302:}{:eb302:}

15951676682

all.f90文档

15951676682

grads 用vint

小懒C

2020，拜蝙蝠所赐，以及不停的用Fortran和grads，写了这个代码，只是量值在某些地区略微有些不对，在定量研究上可能是需要改进的，不是很麻烦的话建议用grads与ncl等好了。

写在前面：
1.各位科研党，经过多张图计算，可以不必计算整层水汽通量，只画850hPa就好了，首先水汽输送都集中在低层，850hPa的图与整层的图在亚洲东部除量级不同外长的和双胞胎一样，而且画图也只是一个辅助你说明的东西，本人个人认为还是少点事只画850好

本人参照南信大PPT以及http://bbs.06climate.com/forum.p ... B%B2%E3%CB%AE%C6%FB
水汽通量单位相关http://bbs.06climate.com/forum.p ... 0&fromuid=44552

代码：该程序只有水汽通量计算，散度涡度读取部分程序可去

1. program main
   
2.   integer,parameter::nx=121,ny=37,nm=12,nz=6,nye=30
   
3.   !deltaphi 经度间距  deltalam 纬度间距  a 地球半径
   
4.   real,parameter::pi=3.14,deltaphi=2.5,deltalam=2.5,a=6371*10**3,g=9.8
   
5.   !i_ '程序中数据维数  cc 高空第cc层，与cc0配合使用
   
6.   integer ix,iy,iz,im,iye
   
7.   real a0
   
8.   character cc0(2)*4
   
9.   data cc0/'850','700'/
   
10.   real::cc(6)=(/1000,925,850,700,600,500/)
    
11.   
    
12.   ! u u风量   v v风量
    
13.   ! T 温度   Td 露点  T_Td 温度露点差
    
14.   ! e 实际水汽压  q 比湿
    
15.   ! qu 水汽通量x分量   qv 水汽通量y分量
    
16.   real u(nx,ny,nz,nm,nye),v(nx,ny,nz,nm,nye)
    
17.   real q(nx,ny,nz,nm,nye)
    
18.   
    
19.   real qu(nx,ny,nz,nm,nye),qv(nx,ny,nz,nm,nye)
    
20.   real quz(nx,ny,nz,nm),qvz(nx,ny,nz,nm)
    
21.   real qua(nx,ny,nm),qva(nx,ny,nm),pu
    
22.   real slp(nx,ny,nm,nye),slpz(nx,ny,nye)
    
23.   !div_q 水汽通量散度  div_qu qu带来的散度  div_qv qv带来的散度  div_qr 曲率项
    
24.   !vor_q 水汽通量涡度  vor_qu qu带来的涡度  vor_qv qv带来的涡度  vor_qr 曲率项
    
25.   !lon 经度  lat 纬度
    
26.   !deltax 经向距离间距  deltay 纬向距离间距)
    
27.   real::div_q(nx,ny,nz,nm,nye)=-9.99e+33,vor_q(nx,ny,nz,nm,nye)=-9.99e+33
    
28.   real::div_qa(nx,ny,nm,nye)=-9.99e+33,vor_qa(nx,ny,nm,nye)=-9.99e+33
    
29.   real div_qu,div_qv,div_qr
    
30.   real vor_qu,vor_qv,vor_qr
    
31.   
    
32.   real::lon(121)=(/((ix-1)*2.5+60,ix=1,nx)/),lat(37)=(/((iy-1)*2.5-10,iy=1,ny)/)
    
33.   open(10,file='E:\datiao2\part4circulation\vaportrans\qV3.grd',form='binary')
    
34.    do iye=1,30
    
35.        do im=1,12
    
36.        do iz=1,nz
    
37.       read(10) ((qu(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)
    
38.        end do
    
39.        do iz=1,nz
    
40.       read(10) ((qv(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)
    
41.        end do
    
42.        end do
    
43.    enddo
    
44.    close(10)
    
45.    open(11,file='E:\datiao2\part4circulation\vaportrans\divvor_q8.grd',form='binary')
    
46.    do iye=1,30
    
47.      do im=1,12
    
48.       do iz=1,nz
    
49.       read(11) ((div_q(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)
    
50.       end do
    
51.       do iz=1,nz
    
52.       read(11) ((vor_q(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)
    
53.       end do
    
54.      end do
    
55.    enddo
    
56.    close(11)
    
57.    div_q=div_q/10**8
    
58.    vor_q=vor_q/10**8
    
59.   open(12,file='E:\datiao2\part4circulation\vaportrans\slp.mon.mean.grd',form='binary')
    
60.     read(12) a0
    
61.   do iye=1,30
    
62.      do im=1,12
    
63.        read(12) ((slp(ix,iy,im,iye),ix=1,nx),iy=1,ny)
    
64.      enddo
    
65.   enddo
    
66.   write(*,*) slp(1,1,1,1)
    
67.   do ix=1,nx
    
68.       do iy=1,ny
    
69.           do iz=1,nz
    
70.               do im=1,nm
    
71.                   quz(ix,iy,iz,im)=sum(qu(ix,iy,iz,im,:))/nye
    
72.                   qvz(ix,iy,iz,im)=sum(qv(ix,iy,iz,im,:))/nye
    
73.                   slpz(ix,iy,im)=sum(slp(ix,iy,im,:))/nye
    
74.               enddo
    
75.           enddo
    
76.       enddo
    
77.   enddo
    
78.   
    
79.   
    
80.       do im=1,nm
    
81.           do iz=1,nz-1
    
82.               do iy=1,ny
    
83.                   do ix=1,nx
    
84.                       if(iz==1.and.slpz(ix,iy,im)>1000) then
    
85.                          pu=quz(ix,iy,iz,im)*(slpz(ix,iy,im)-cc(iz))*100
    
86.                       else if(iz>1) then
    
87.                           pu=quz(ix,iy,iz,im)*(cc(iz-1)-cc(iz))*100
    
88.                       endif
    
89.                       qua(ix,iy,im)=qua(ix,iy,im)+pu
    
90.                   enddo
    
91.               enddo
    
92.           enddo
    
93.       enddo
    
94.   write(*,*) qua(:,:,1)
    
95. 
    
96.       do im=1,nm
    
97.           do iz=1,nz
    
98.               do iy=1,ny
    
99.                   do ix=1,nx
    
100.                       if(iz==1.and.slpz(ix,iy,im)>1000) then
     
101.                          pu=qvz(ix,iy,iz,im)*(slpz(ix,iy,im)-cc(iz))*100
     
102.                       else if(iz>1) then
     
103.                           pu=qvz(ix,iy,iz,im)*(cc(iz-1)-cc(iz))*100
     
104.                       endif
     
105.                       qva(ix,iy,im)=qva(ix,iy,im)+pu
     
106.                   enddo
     
107.               enddo
     
108.           enddo
     
109.       enddo
     
110. 
     
111.   write(*,*) qva(:,:,1)
     
112.   open(10,file='E:\datiao2\part4circulation\vaportrans\quva.grd',form='binary')
     
113.        do im=1,12
     
114.       write(10) ((qua(ix,iy,im),ix=1,nx),iy=1,ny)
     
115.       write(10) ((qva(ix,iy,im),ix=1,nx),iy=1,ny)
     
116.        end do
     
117.    close(10)
     
118.     end
     
119.    

program main<br />   integer,parameter::nx=121,ny=37,nm=12,nz=6,nye=30<br />   !deltaphi 经度间距  deltalam 纬度间距  a 地球半径<br />   real,parameter::pi=3.14,deltaphi=2.5,deltalam=2.5,a=6371*10**3,g=9.8<br />   !i_ '程序中数据维数  cc 高空第cc层，与cc0配合使用<br />   integer ix,iy,iz,im,iye<br />   real a0<br />   character cc0(2)*4<br />   data cc0/'850','700'/<br />   real::cc(6)=(/1000,925,850,700,600,500/)<br />   <br />   ! u u风量   v v风量<br />   ! T 温度   Td 露点  T_Td 温度露点差<br />   ! e 实际水汽压  q 比湿<br />   ! qu 水汽通量x分量   qv 水汽通量y分量<br />   real u(nx,ny,nz,nm,nye),v(nx,ny,nz,nm,nye)<br />   real q(nx,ny,nz,nm,nye)<br />   <br />   real qu(nx,ny,nz,nm,nye),qv(nx,ny,nz,nm,nye)<br />   real quz(nx,ny,nz,nm),qvz(nx,ny,nz,nm)<br />   real qua(nx,ny,nm),qva(nx,ny,nm),pu<br />   real slp(nx,ny,nm,nye),slpz(nx,ny,nye)<br />   !div_q 水汽通量散度  div_qu qu带来的散度  div_qv qv带来的散度  div_qr 曲率项<br />   !vor_q 水汽通量涡度  vor_qu qu带来的涡度  vor_qv qv带来的涡度  vor_qr 曲率项<br />   !lon 经度  lat 纬度 <br />   !deltax 经向距离间距  deltay 纬向距离间距)<br />   real::div_q(nx,ny,nz,nm,nye)=-9.99e+33,vor_q(nx,ny,nz,nm,nye)=-9.99e+33<br />   real::div_qa(nx,ny,nm,nye)=-9.99e+33,vor_qa(nx,ny,nm,nye)=-9.99e+33<br />   real div_qu,div_qv,div_qr<br />   real vor_qu,vor_qv,vor_qr<br />   <br />   real::lon(121)=(/((ix-1)*2.5+60,ix=1,nx)/),lat(37)=(/((iy-1)*2.5-10,iy=1,ny)/)<br />   open(10,file='E:\datiao2\part4circulation\vaportrans\qV3.grd',form='binary')<br />    do iye=1,30<br />        do im=1,12<br />        do iz=1,nz<br />       read(10) ((qu(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)<br />        end do<br />        do iz=1,nz<br />       read(10) ((qv(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)<br />        end do<br />        end do<br />    enddo<br />    close(10)<br />    open(11,file='E:\datiao2\part4circulation\vaportrans\divvor_q8.grd',form='binary')<br />    do iye=1,30<br />      do im=1,12<br />       do iz=1,nz<br />       read(11) ((div_q(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)<br />       end do<br />       do iz=1,nz<br />       read(11) ((vor_q(ix,iy,iz,im,iye),ix=1,nx),iy=1,ny)<br />       end do<br />      end do<br />    enddo <br />    close(11)<br />    div_q=div_q/10**8<br />    vor_q=vor_q/10**8<br />   open(12,file='E:\datiao2\part4circulation\vaportrans\slp.mon.mean.grd',form='binary')<br />     read(12) a0<br />   do iye=1,30<br />      do im=1,12<br />        read(12) ((slp(ix,iy,im,iye),ix=1,nx),iy=1,ny)<br />      enddo<br />   enddo<br />   write(*,*) slp(1,1,1,1)<br />   do ix=1,nx<br />       do iy=1,ny<br />           do iz=1,nz<br />               do im=1,nm<br />                   quz(ix,iy,iz,im)=sum(qu(ix,iy,iz,im,:))/nye<br />                   qvz(ix,iy,iz,im)=sum(qv(ix,iy,iz,im,:))/nye<br />                   slpz(ix,iy,im)=sum(slp(ix,iy,im,:))/nye<br />               enddo<br />           enddo<br />       enddo<br />   enddo<br />   <br />   <br />       do im=1,nm<br />           do iz=1,nz-1<br />               do iy=1,ny<br />                   do ix=1,nx<br />                       if(iz==1.and.slpz(ix,iy,im)>1000) then <br />                          pu=quz(ix,iy,iz,im)*(slpz(ix,iy,im)-cc(iz))*100<br />                       else if(iz>1) then<br />                           pu=quz(ix,iy,iz,im)*(cc(iz-1)-cc(iz))*100<br />                       endif<br />                       qua(ix,iy,im)=qua(ix,iy,im)+pu<br />                   enddo<br />               enddo<br />           enddo<br />       enddo<br />   write(*,*) qua(:,:,1)<br /> <br />       do im=1,nm<br />           do iz=1,nz<br />               do iy=1,ny<br />                   do ix=1,nx<br />                       if(iz==1.and.slpz(ix,iy,im)>1000) then <br />                          pu=qvz(ix,iy,iz,im)*(slpz(ix,iy,im)-cc(iz))*100<br />                       else if(iz>1) then<br />                           pu=qvz(ix,iy,iz,im)*(cc(iz-1)-cc(iz))*100<br />                       endif<br />                       qva(ix,iy,im)=qva(ix,iy,im)+pu<br />                   enddo<br />               enddo<br />           enddo<br />       enddo<br /> <br />   write(*,*) qva(:,:,1)<br />   open(10,file='E:\datiao2\part4circulation\vaportrans\quva.grd',form='binary')<br />        do im=1,12<br />       write(10) ((qua(ix,iy,im),ix=1,nx),iy=1,ny)<br />       write(10) ((qva(ix,iy,im),ix=1,nx),iy=1,ny)<br />        end do<br />    close(10)<br />     end<br />    

只积分到500hPa这个程序，无计算散度的，散度依然用到的是grads画出的

图的填充都是水汽通量散度，矢量箭头为水汽通量