倒算法计算大气热源的数据与程序问题

云一林 · 发表于 2020-2-8 22:20:41

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

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

x

谁能帮我看看这个倒算法计算大气热源程序哪里出错了么？程序来源于气象家园论坛，Fortran运行没有报错但是出不来结果。不知道是不是用的数据不对，用的数据是：

pres.surf.mon.mean.dat（地面气压数据，单位pa,高度上只有1层）

uwnd.mon.mean.dat(经向风，单位m/s，高度上有27层）

vwnd.mon.mean.dat（纬向风，单位m/s，高度上有27层）

tmp.mon.mean.dat（温度，单位K，高度上有27层）

shum.mon.mean.dat（比湿度，单位kg/kg，高度上有27层）

omega.mon.mean.dat（p坐标上的垂直速度，单位pa/s，高度上有27层）

所用数据来源于jra55，以下是各个数据的ctl：

dset E:\lunwen\pres.surf.mon.mean.dat

undef 1e+20

xdef 144 linear 0 2.5

ydef 73 linear -90 2.5

zdef 1 linear 0 1

tdef 672 linear 00Z01JAN1958 1mo

vars 1

pres 0 99 *Pressure

Endvars

dset E:\lunwen\uwnd.mon.mean.dat

undef 1e+20

xdef 144 linear 0 2.5

ydef 73 linear -90 2.5

zdef 27 levels 1000 975 950 925 900 875 850 825 800 775 750 700 650 600 550 500 450 400 350 300 250 225 200 175 150 125 100

tdef 672 linear 00Z01JAN1958 1mo

vars 1

uwnd 27 99 * u-component of wind

Endvars

dset E:\lunwen\vwnd.mon.mean.dat

undef 1e+20

xdef 144 linear 0 2.5

ydef 73 linear -90 2.5

zdef 27 levels 1000 975 950 925 900 875 850 825 800 775 750 700 650 600 550 500 450 400 350 300 250 225 200 175 150 125 100

tdef 672 linear 00Z01JAN1958 1mo

vars 1

vwnd 27 99 * v-component of wind

Endvars

dset E:\lunwen\tmp.mon.mean.dat

undef 1e+20

xdef 144 linear 0 2.5

ydef 73 linear -90 2.5

zdef 27 levels 1000 975 950 925 900 875 850 825 800 775 750 700 650 600 550 500 450 400 350 300 250 225 200 175 150 125 100

tdef 672 linear 00Z01JAN1958 1mo

vars 1

tmp 27 99 *Temperature

Endvars

dset E:\lunwen\shum.mon.mean.dat

undef 1e+20

xdef 144 linear 0 2.5

ydef 73 linear -90 2.5

zdef 27 levels 1000 975 950 925 900 875 850 825 800 775 750 700 650 600 550 500 450 400 350 300 250 225 200 175 150 125 100

tdef 672 linear 00Z01JAN1958 1mo

vars 1

shum 27 99 *Specific humidity

Endvars

dset E:\lunwen\omega.mon.mean.dat

undef 1e+20

xdef 144 linear 0 2.5

ydef 73 linear -90 2.5

zdef 27 levels 1000 975 950 925 900 875 850 825 800 775 750 700 650 600 550 500 450 400 350 300 250 225 200 175 150 125 100

tdef 672 linear 00Z01JAN1958 1mo

vars 1

omega 27 99 *Pressure vertical velocity

Endvars

经过grads检验，以上ctl是正确的数据描述文件。以下为fortran程序：

!the yrev must be noted.

implicit none

integer,parameter :: NX=144,NY=73,NL=27,NT=672

real ,parameter :: PI=3.1415926,RR=6371229.0

real ,parameter :: RAD=PI/180.0,DD=RR*RAD*2.5,ck=0.286,gama=0.0065

real ,parameter :: rd=287.0,g=9.80616, dk=g/(rd*gama), WTt=0.0

real ,parameter :: DT=24*3600.0, CP=1004.0, lh=2501000.0

!INPUT

logical :: yrev=.false.,smooth=.FALSE. !yrev=.true. 数据中y是反向的，如果yrev=.false.,则不改变方向

real :: pres(NX,NY)

real :: u(nx,ny,nl),v(nx,ny,nl),t(nx,ny,nl),W(NX,NY,NL),q(NX,NY,NL) ! H(NX,NY,NL)

!TEMP

REAL :: temp(nx,73,nl),OROG1(NX,73)

REAL :: gtmp(nx,ny,nl)

real :: TH(NX,NY,NL),tf(nx,ny,nl)

real :: ft1(nx,ny,nl),ft3(nx,ny,nl)

REAL :: FF(NX,NY),GG(NX,NY)

real :: ILAY(NX,NY)

REAL :: PP(NL)=(/1000.,975.,950.,925.,900.,875.,850.,825.,800.,775.,750.,700.,650.,600.,550.,500.,450.,400.,350.,300.,250.,225.,200.,175.,150.,125.,100./)

!OUTPUT

real :: ft(nx,ny,nl),fv(nx,ny,nl),fw(nx,ny,nl),q2(nx,ny)

real :: Q1(NX,NY),HT(NX,NY,NL),ffG(nx,ny,nl)

real :: wg(nx,ny,nl),ug(nx,ny,nl),tg(nx,ny,nl)

!PROCESS

integer :: I,J,K,L,ip

INTEGER :: IREC

real :: dx,dy,dp,ss,sx,si,sj,l1

CHARACTER*42 :: FILEPATH='E:\lunwen\'

open(11,file=FILEPATH//'pres.surf.mon.mean.dat',form='unformatted',access='direct',action='read',recl=nx*73)

OPEN(13,FILE=FILEPATH//'uwnd.mon.mean.dat',FORM='UNFORMATTED',ACCESS='DIRECT',action='read',RECL=NX*73)

OPEN(14,FILE=FILEPATH//'vwnd.mon.mean.dat',FORM='UNFORMATTED',ACCESS='DIRECT',action='read',RECL=NX*73)

OPEN(15,FILE=FILEPATH//'tmp.mon.mean.dat',FORM='UNFORMATTED',ACCESS='DIRECT',action='read',RECL=NX*73)

OPEN(16,FILE=FILEPATH//'shum.mon.mean.dat',FORM='UNFORMATTED',ACCESS='DIRECT',action='read',RECL=NX*73)

OPEN(17,FILE=FILEPATH//'omega.mon.mean.dat',FORM='UNFORMATTED',ACCESS='DIRECT',action='read',RECL=NX*73)

!write data into:

OPEN(20,FILE=FILEPATH//'Q\Q1TNWNS.grd',FORM='binary')

OPEN(21,FILE=FILEPATH//'Q\q1vnwNS.grd',FORM='binary')

OPEN(22,FILE=FILEPATH//'Q\q1wnwNS.grd',FORM='binary')

OPEN(24,FILE=FILEPATH//'Q\q2tnwNS.grd',FORM='binary')

OPEN(25,FILE=FILEPATH//'Q\q2vnwNS.grd',FORM='binary')

OPEN(26,FILE=FILEPATH//'Q\q2wnwNS.grd',FORM='binary')

OPEN(27,FILE=FILEPATH//'Q\q1&2nwNS.grd',FORM='binary')

OPEN(28,FILE=FILEPATH//'Q\q1l17nwNS.grd',FORM='binary')

OPEN(29,FILE=FILEPATH//'Q\q2l17nwNS.grd',FORM='binary')

LOOPDAY: DO K=1,nt

print *,'Begin time = ',k

HT =-9.99E33 !Q1

ffG =-9.99e33 !Q2

tg =-9.99e33 !(1.1)

ug =-9.99e33 !(1.2)

wg =-9.99e33 !(1.3)

ft =-9.99e33 !(2.1)

fv =-9.99e33 !(2.2)

fw =-9.99e33 !(2.3)

!!!!!get data from files

read (11,rec=k)((orog1(i,j),i=1,nx),j=1,73)

pres(:,1:71)=orog1(:,2:72)

if (yrev) then

pres(:,1:71)=pres(:,71:1:-1)

endif

LOOPLEV1: DO L=1,NL

IREC=L+NL*(K-1)

READ(13,REC=IREC)((temp(I,J,L),I=1,NX),J=1,73)

u(:,1:71,L)=temp(:,2:72,l)

READ(14,REC=IREC)((temp(I,J,L),I=1,NX),J=1,73)

v(:,1:71,l)=temp(:,2:72,l)

READ(15,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

t(:,1:71,l)=temp(:,2:72,l)

READ(17,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

w(:,1:71,l)=temp(:,2:72,l)

if (yrev) then

! h(:,1:71,l) = h(:,71:1:-1,l)

u(:,1:71,l) = u(:,71:1:-1,l)

v(:,1:71,l) = v(:,71:1:-1,l)

t(:,1:71,l) = t(:,71:1:-1,l)

w(:,1:71,l) = w(:,71:1:-1,l)

endif

if (smooth) then

FF(:,:)=u(:,:,l)

call smth9(nx,ny,FF,GG)

u(:,:,l)=gg(:,:)

FF(:,:)=v(:,:,l)

call smth9(nx,ny,FF,GG)

v(:,:,l)=GG(:,:)

FF(:,:)=t(:,:,l)

call smth9(nx,ny,FF,GG)

t(:,:,l)=GG(:,:)

FF(:,:)=w(:,:,l)

call smth9(nx,ny,FF,GG)

w(:,:,l)=GG(:,:)

endif

ENDDO looplev1

loopq1:do l=1,nl

irec=l+(k-1)*nl

READ(16,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

q(:,1:71,l)=temp(:,2:72,l)

if (yrev) then

q(:,1:71,l) = q(:,71:1:-1,l)

endif

if (smooth) then

FF(:,:)=q(:,:,l)

call smth9(nx,ny,FF,GG)

q(:,:,l)=GG(:,:)

endif

enddo loopq1

IF(K > 1)THEN

looplev2: DO L=1,NL

irec=l+nl*(k-2)

READ(15,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

tf(:,1:71,l)=temp(:,2:72,l)

if (yrev) then

tf(:,1:71,l) = tf(:,71:1:-1,l)

endif

if (smooth) then

FF(:,:)=tf(:,:,l)

call smth9(nx,ny,FF,GG)

tf(:,:,l)=GG(:,:)

endif

ENDDO looplev2

loopq2:do l=1,nl

irec=l+nl*(k-2)

READ(16,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

ft1(:,1:71,l)=temp(:,2:72,l)

if (yrev) then

ft1(:,1:71,l) = ft1(:,71:1:-1,l)

endif

if (smooth) then

FF(:,:)=ft1(:,:,l)

call smth9(nx,ny,FF,GG)

ft1(:,:,l)=GG(:,:)

endif

enddo loopq2

ENDIF

IF(K < NT)THEN

looplev3: DO L=1,NL

irec=l+nl*k

READ(15,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

th(:,1:71,l)=temp(:,2:72,l)

if (yrev) then

th(:,1:71,l) = th(:,71:1:-1,l)

endif

if (smooth) then

FF(:,:)=TH(:,:,l)

call smth9(nx,ny,FF,GG)

th(:,:,l)=GG(:,:)

endif

ENDDO looplev3

loopq3:do l=1,nl

irec=l+nl*k

READ(16,REC=irec)((temp(I,J,L),I=1,NX),J=1,73)

ft3(:,1:71,l)=temp(:,2:72,l)

if(yrev) then

ft3(:,1:71,l) = ft3(:,71:1:-1,l)

endif

if (smooth) then

FF(:,:)=ft3(:,:,l)

call smth9(nx,ny,FF,GG)

ft3(:,:,l)=GG(:,:)

endif

enddo loopq3

ENDIF

! 计算高于地形的最低层

do j=1,ny

do i=1,nx

do l=1,nl

l1=pp(l)*100.0

if(pres(i,j) >= l1)then

ip=l

exit

endif

ENDDO

ILAY(I,J)=IP

enddo

! 计算温度的局地变化 (1.1,2.1)

IF(K == 1)THEN

do j=1,ny

do i=1,nx

ip=ilay(i,j)

do l=ip,nl

TG(i,j,l)=(TH(i,j,l)-T(i,j,l))/DT

enddo

do l=ip,nl

fT(i,j,l)=(ft3(i,j,l)-q(i,j,l))/DT

enddo

ELSE IF(K == NT)THEN

do j=1,ny

do i=1,nx

ip=ilay(i,j)

do l=ip,nl

TG(i,j,l)=(T(i,j,l)-Tf(i,j,l))/DT

enddo

do l=ip,nl

fT(i,j,l)=(q(i,j,l)-fT1(i,j,l))/DT

enddo

ELSE

do j=1,ny

do i=1,nx

ip=ilay(i,j)

do l=ip,nl

TG(i,j,l)=(TH(i,j,l)-Tf(i,j,l))/2.0/DT

enddo

do l=ip,nl

fT(i,j,l)=(fT3(i,j,l)-fT1(i,j,l))/2.0/DT

enddo

ENDIF

! 计算温度平流---中央差分 (1.2,2.2)

do i=2,nx-1

do j=2,ny-1

dx=cos((-87.5+2.5*(j-1))*RAD)

dy=2*dd

ip=ilay(i,j)

DO L=ip,NL

UG(I,J,L)=(U(I,J,L)*(T(I+1,J,L)-T(I-1,J,L))/dx &

+(V(I,J,L)*(T(I,J+1,L)-T(I,J-1,L))))/dy

ENDDO

do l=ip,nl

fv(I,J,L)=(U(I,J,L)*(q(I+1,J,L)-q(I-1,J,L))/dx &

+(V(I,J,L)*(q(I,J+1,L)-q(I,J-1,L))))/dy

enddo

ENDDO

do j=2,ny-1

dx=cos((-87.5+2.5*(j-1))*RAD)

dy=2*dd

ip=ilay(1,j)

DO L=ip,NL

UG(1,J,L)=(U(1,J,L)*(T(2,J,L)-T(nx,J,L))/dx &

+(V(1,J,L)*(T(1,J+1,L)-T(1,J-1,L))))/dy

ENDDO

do l=ip,nl

fv(1,J,L)=(U(1,J,L)*(q(2,J,L)-q(nx,J,L))/dx &

+(V(1,J,L)*(q(1,J+1,L)-q(1,J-1,L))))/dy

enddo

ip=ilay(nx,j)

DO L=ip,NL

UG(nx,J,L)=(U(nx,J,L)*(T(1,J,l)-T(nx-1,J,L))/dx &

+(V(NX,J,L)*(T(NX,J+1,L)-T(NX,J-1,L))))/dy

enddo

do l=ip,nl

fv(nx,J,L)=(U(nx,J,L)*(q(1,J,l)-q(nx-1,J,L))/dx+ &

+(V(NX,J,L)*(q(NX,J+1,L)-q(NX,J-1,L))))/dy

enddo

ENDDO

! 计算位温

DO J=1,NY

DO I=1,NX

DO L=1,NL

gtmp(i,j,l)=t(i,j,l)*(1000.0/pp(l))**ck

ENDDO

enddo

do i=1,nx

do j=2,ny-1

IP=ILAY(I,J)

! 计算位温垂直输送 (1.3,2.3)

do l=ip+1,nl-1

dp=pp(l+1)-pp(l-1)

wg(i,j,l)=W(I,J,L)*(GTMP(I,J,L+1)-GTMP(I,J,L-1))/dp

ENDDO

do l=ip+1,7

dp=pp(l+1)-pp(l-1)

fw(i,j,l)=W(I,J,L)*(q(I,J,L+1)-q(I,J,L-1))/DP

enddo

dp=pp(ip+1)-pp(ip)

WG(I,J,IP)=W(I,J,IP)*(GTMP(I,J,IP+1)-GTMP(I,J,IP))/DP

fW(I,J,IP)=W(I,J,IP)*(q(I,J,IP+1)-q(I,J,IP))/DP

dp=pp(nl)-pp(nl-1)

WG(I,J,NL)=W(I,J,NL)*(GTMP(I,J,NL)-GTMP(I,J,NL-1))/DP

dp=pp(27)-pp(26)

fW(I,J,27)=W(I,J,27)*(q(I,J,27)-q(I,J,26))/DP

ENDDO

! 计算加热率

DO I=1,NX

DO J=2,NY-1

IP=ILAY(I,J)

DO L=IP,NL

wg(i,j,l)=WG(I,J,L)*(PP(L)/1000.0)**CK

HT(I,J,L)=TG(I,J,L)+UG(I,J,L)+WG(I,J,L)

ENDDO

do l=ip,nl

ft(i,j,l)=-ft(i,j,l)*lh/cp

fv(i,j,l)=-fv(i,j,l)*lh/cp

fw(i,j,l)=-fw(i,j,l)*lh/cp

ffG(I,J,L)=fT(I,J,L)+fv(I,J,L)+fW(I,J,L)

enddo

SS=0.0

sx=0.0

do L=IP,NL-1

dp=pp(l+1)-pp(l)

SI=(HT(I,J,L)+HT(I,J,L+1))*DP/2.0

SS=SS-SI

ENDDO

do l=ip,7

dp=pp(l+1)-pp(l)

Sj=(ffG(I,J,L)+ffG(I,J,L+1))*DP/2.0

Sx=Sx-Sj

enddo

Q1(I,J)=SS/G*100.0*CP

q2(i,j)=sx/g*100.0*cp

ENDDO

! write(30,rec=irec)(((gtmp(i,j,l),i=1,nx),j=2,ny-1),l=1,nl)

write(20)(((tg(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(21)(((ug(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(22)(((wg(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(24)(((ft(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(25)(((fv(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(26)(((fw(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(28)(((ht(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

write(29)(((ffG(I,J,L),I=1,NX),J=2,NY-1),l=1,nl)

!

write(27)((q1(I,J),I=1,NX),J=2,NY-1)

write(27)((q2(I,J),I=1,NX),J=2,NY-1)

WRITE(6,*)'end TIME=',K

ENDDO LOOPDAY

CLOSE(12);CLOSE(13);CLOSE(14);CLOSE(15);close(16)

CLOSE(20);CLOSE(21);CLOSE(22);close(24);close(25)

close(26);close(27);close(28);close(29)

stop

end

subroutine smth9(nx,ny,ff,gg)

real :: ff(nx,ny),gg(nx,ny)

integer:: i,j

real :: s

do i=2,nx-1

do j=2,ny-1

s=0.5

gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(i+1,j)+ff(i-1,j)+ff(i,j+1) &

+ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(i+1,J+1)+ff(i-1,j-1)+ &

ff(i+1,j-1)+ff(i-1,j+1)-4*ff(i,j))

enddo

i=1

do j=2,ny-1

gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(i+1,j)+ff(nx,j)+ff(i,j+1) &

+ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(i+1,J+1)+ff(nx,j-1)+ &

ff(i+1,j-1)+ff(nx,j+1)-4*ff(i,j))

enddo

i=nx

do j=2,ny-1

gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(1,j)+ff(i-1,j)+ff(i,j+1) &

+ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(1,J+1)+ff(i-1,j-1)+ &

ff(i,j-1)+ff(i-1,j+1)-4*ff(i,j))

enddo

gg(:,1)=ff(:,1)

gg(:,ny)=ff(:,ny)

return

end

修改部分已经标记出来，fortran程序运行结果：

forrtl: severe (29): file not found, unit 11, file E:\lunwen\

pres.surf.mon.mean.dat

Image PC Routine Line Source

q1q2.exe 0041DCD9 Unknown Unknown Unknown

Incrementally linked image--PC correlation disabled.

Press any key to continue

翻译为找不到文件pres.surf.mon.mean.dat，为什么会找不到，是因为这个文件高度层与别的文件不同吗？只有一层，还是说我使用的原始数据不正确？？请问谁做出来了可以告诉我一声吗？是数据用的不对还是说fortran程序还有些地方我没有改正？？

Houckin · 发表于 2020-4-12 19:45:46

请问楼主搞定了没

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

[求助] 倒算法计算大气热源的数据与程序问题

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