Fortran实习：计算温度露点差水平分布、位温，相当位温水平分布、涡度，散度等

WesleyMoob

program name
    implicit none
    integer,parameter :: nx=53,ny=29,nz=11,nt=2
    real,parameter :: delta=2.5,pi=3.1415926,a=6371000,a1=17.2693882,b1=35.86,a2=21.8745584
    real,parameter :: b2=7.66
    integer i,j,k,t,l
    real u(nx,ny,nz,nt),v(nx,ny,nz,nt),fi(ny),lat(ny),vor(nx,ny,nz,nt),div(nx,ny,nz,nt),omg(nx,ny,nz,nt)
    real omgg(nx,ny,nz,nt),dp(2:nz),zeta(nx,ny,nz,nt),ttd(nx,ny,nz,nt),temp(nx,ny,nz,nt),p(nz),td(nx,ny,nz,nt)
    real q(nx,ny,nz,nt),dtx(nx,ny,nz,nt),dty(nx,ny,nz,nt),dwx(nx,ny,nz,nt),dwy(nx,ny,nz,nt),ht(nx,ny,nz,nt)
    real e(nx,ny,nz,nt),a3(nx,ny,nz,nt),b3(nx,ny,nz,nt),tpl(nx,ny,nz,nt),wdpl(nx,ny,nz,nt),cp(nx,ny,nz,nt)
    real tl(nx,ny,nz,nt),el(nx,ny,nz,nt),ql(nx,ny,nz,nt),zl(nx,ny,nz,nt),se(nx,ny,nz,nt)
    CHARACTER timefile(2)*12 ,levelfile(11)*4
    real lon,lamda,dhx,dhy,duy,dvx,dux,dvy,M,pl
    DATA levelfile/'1000','925','850','700','500','400','300','250','200','150','100'/
    DATA timefile/'13052520.000','13052620.000'/
    p=(/1000,925,850,700,500,400,300,250,200,150,100/)
    !读取u,v数据(m/s)
    do t=1,nt
        do k=1,nz
            open(10,file='/Users/wesleymoob/Desktop/uv/'//trim(levelfile(k))//'/'//trim(timefile(t)))
                do l=1,3
                    read(10,*)
                enddo
                do j=ny,1,-1
                    read(10,*)(u(i,j,k,t),i=1,nx)        
                end do
                do j=ny,1,-1
                    read(10,*)(v(i,j,k,t),i=1,nx)        
                end do
                close(10)
        end do
    end do
    !计算经纬间隔
    lamda=pi/72.
    lon=a*lamda !dy
    do j=1,ny
        fi(j)=((10+(j-1)*delta)*pi)/180
        lat(j)=a*cos(fi(j))*lamda !dx
    end do
    !涡度
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    if(i==1 .and. j==1)then
                        duy=(u(i,j+1,k,t)-u(i,j,k,t))/lon
                        dvx=(v(i+1,j,k,t)-v(i,j,k,t))/lat(j)
                    else if(i==nx.and.j==ny)then
                        duy=(u(i,j,k,t)-u(i,j-1,k,t))/lon
                        dvx=(v(i,j,k,t)-v(i-1,j,k,t))/lat(j)
                    else if(i==1 .and. j==ny)then
                        duy=(u(i,j,k,t)-u(i,j-1,k,t))/lon
                        dvx=(v(i+1,j,k,t)-v(i,j,k,t))/lat(j)
                    else if(i==nx.and.j==1)then
                        duy=(u(i,j+1,k,t)-u(i,j,k,t))/lon
                        dvx=(v(i,j,k,t)-v(i-1,j,k,t))/lat(j)
                    else if(i==1 .and. j/=1 .and. j/=ny)then
                        duy=(u(i,j+1,k,t)-u(i,j-1,k,t))/(2*lon)
                        dvx=(v(i+1,j,k,t)-v(i,j,k,t))/lat(j)
                    else if(i==nx.and.j/=1.and.j/=ny)then
                        duy=(u(i,j+1,k,t)-u(i,j-1,k,t))/(2*lon)
                        dvx=(v(i,j,k,t)-v(i-1,j,k,t))/lat(j)
                    else if(j==1.and.i/=1.and.i/=nx)then
                        duy=(u(i,j+1,k,t)-u(i,j,k,t))/lon
                        dvx=(v(i+1,j,k,t)-v(i-1,j,k,t))/(2.*lat(j))
                    else if(j==ny.and.i/=1.and.i/=nx)then
                        duy=(u(i,j,k,t)-u(i,j-1,k,t))/lon
                        dvx=(v(i+1,j,k,t)-v(i-1,j,k,t))/(2.*lat(j))
                    else
                        duy=(u(i,j+1,k,t)-u(i,j-1,k,t))/(2.*lon)
                        dvx=(v(i+1,j,k,t)-v(i-1,j,k,t))/(2.*lat(j))
                    end if
                    vor(i,j,k,t)=dvx-duy+((u(i,j,k,t)/a)*tan(fi(j)))  
                end do
            end do
        end do
    end do
    open(20,file="/Users/wesleymoob/Desktop/vor.grd",form='unformatted',access='stream')
    write(20)vor
    close(20)
    !散度
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    if(i==1 .and. j==1)then
                        dux=(u(i+1,j,k,t)-u(i,j,k,t))/lat(j)
                        dvy=(v(i,j+1,k,t)-v(i,j,k,t))/lon
                    else if(i==nx.and.j==ny)then
                        dux=(u(i,j,k,t)-u(i-1,j,k,t))/lat(j)
                        dvy=(v(i,j,k,t)-v(i,j-1,k,t))/lon
                    else if(i==1 .and. j==ny)then
                        dux=(u(i+1,j,k,t)-u(i,j,k,t))/lat(j)
                        dvy=(v(i,j,k,t)-v(i,j-1,k,t))/lon
                    else if(i==nx.and.j==1)then
                        dux=(u(i,j,k,t)-u(i-1,j,k,t))/lat(j)
                        dvy=(v(i,j+1,k,t)-v(i,j,k,t))/lon
                    else if(i==1 .and. j/=1 .and. j/=ny)then
                        dux=(u(i+1,j,k,t)-u(i,j,k,t))/lat(j)
                        dvy=(v(i,j+1,k,t)-v(i,j-1,k,t))/(2.*lon)
                    else if(i==nx.and.j/=1.and.j/=ny)then
                        dux=(u(i,j,k,t)-u(i-1,j,k,t))/lat(j)
                        dvy=(v(i,j+1,k,t)-v(i,j-1,k,t))/(2.*lon)
                    else if(j==1.and.i/=1.and.i/=nx)then
                        dux=(u(i+1,j,k,t)-u(i-1,j,k,t))/(2.*lat(j))
                        dvy=(v(i,j+1,k,t)-v(i,j,k,t))/lon
                    else if(j==ny.and.i/=1.and.i/=nx)then
                        dux=(u(i+1,j,k,t)-u(i-1,j,k,t))/(2.*lat(j))
                        dvy=(v(i,j,k,t)-v(i,j-1,k,t))/lon
                    else
                        dux=(u(i+1,j,k,t)-u(i-1,j,k,t))/(2.*lat(j))
                        dvy=(v(i,j+1,k,t)-v(i,j-1,k,t))/(2.*lon)
                    end if
                    div(i,j,k,t)= dux+dvy-((v(i,j,k,t)/a)*tan(fi(j)))
                end do
            end do
        end do
    end do
    open(30,file="/Users/wesleymoob/Desktop/div.grd",form='unformatted',access='stream')
    write(30)div
    close(30)
    !垂直速度(m/s)
    do t=1,nt
        omg(:,:,1,t)=0
    end do
    dp=(/75,75,150,200,100,100,50,50,50,50/)
    M=(nz+1)*nz*0.5
    do t=1,nt
        do k=2,nz
            omg(:,:,k,t)=omg(:,:,k-1,t)+((div(:,:,k-1,t)+div(:,:,k,t))*dp(k))/2.
        end do
    end do
    !垂直速度订正(m/s)
    M=0.5*(nz+1)*nz
    do t=1,nt
        do k=1,nz
            omgg(:,:,k,t)=omg(:,:,k,t)-(((k+1)*k)/real(2*M))*(omg(:,:,nz,t)-0)
        end do
    end do
    open(40,file="/Users/wesleymoob/Desktop/omg.grd",form='unformatted',access='stream')
    write(40)omgg
    close(40)
    !读取高度场
    do t=1,nt
        do k=1,nz
            open(11,file='/Users/wesleymoob/Desktop/height/'//trim(levelfile(k))//'/'//trim(timefile(t)))
            do l=1,4
                read(11,*)
            enddo
            do j=ny,1,-1
                read(11,*)(ht(i,j,k,t),i=1,nx)        
            end do
            close(11)
        end do
    end do
    !读取温度露点差(摄氏度)
    do t=1,nt
        do k=1,nz
            open(50,file='/Users/wesleymoob/Desktop/t-td/'//trim(levelfile(k))//'/'//trim(timefile(t)))
            do l=1,4
                read(50,*)
            enddo
            do j=ny,1,-1
                read(50,*)(ttd(i,j,k,t),i=1,nx)        
            end do
            close(50)
        end do
    end do
    open(60,file="/Users/wesleymoob/Desktop/ttd.grd",form='unformatted',access='stream')
    write(60)ttd
    close(60)
    !读取温度(摄氏度)
    do t=1,nt
        do k=1,nz
            open(70,file='/Users/wesleymoob/Desktop/temper/'//trim(levelfile(k))//'/'//trim(timefile(t)))
            do l=1,4
                read(70,*)
            enddo
            do j=ny,1,-1
                read(70,*)(temp(i,j,k,t),i=1,nx)        
            end do
            close(70)
        end do
    end do
    !计算位温(K)
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    zeta(i,j,k,t)=(temp(i,j,k,t)+273.15)*(1000./p(k))**0.286
                end do
            end do
        end do
    end do
    open(80,file="/Users/wesleymoob/Desktop/zeta.grd",form='unformatted',access='stream')
    write(80)zeta
    close(80)
    !计算露点温度
    td=temp-ttd
    !计算水汽压
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    e(i,j,k,t)=6.11*10**((7.5*td(i,j,k,t))/(273.3+td(i,j,k,t)))
                end do
            end do
        end do
    end do
    !计算比湿(g/g)
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    q(i,j,k,t)=(0.622*e(i,j,k,t))/(p(k)-(0.378*e(i,j,k,t)))
                end do
            end do
        end do
    end do
    !计算Cp
    cp=1004*(1+0.85*q)
    !计算抬升凝结高度以及抬升凝结温度
    !用1000hPa高度上的个参量计算气团q的真值，假设凝结高度为300hPa

    do t=1,nt
        do k=1,nz
            do i=1,nx
                do j=1,ny
                    pl=1.0
                    ql(i,j,k,t)=0
                    do while(ql(i,j,k,t)-q(i,j,k,t)<0)
                        tl(i,j,k,t)=((zeta(i,j,k,t)*(pl/1000.)**0.286)-273.15)
                        el(i,j,k,t)=6.11*10**((7.5*tl(i,j,k,t))/(273.3+tl(i,j,k,t)))
                        ql(i,j,k,t)=(0.622*el(i,j,k,t))/(pl-(0.378*el(i,j,k,t)))
                        pl=pl+1
                    end do
                    zl(i,j,k,t)=(9.8*ht(i,j,k,t)+cp(i,j,k,t)*(temp(i,j,k,t)+273.15)-cp(i,j,k,t)*(tl(i,j,k,t)+273.15))/9.8
                end do
            end do
        end do
    end do
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    se(i,j,k,t)=(temp(i,j,k,t)+273.15)*(1000./(p(k)-e(i,j,k,t)))**0.286&
                    *exp(zl(i,j,k,t)*q(i,j,k,t)/(cp(i,j,k,t)*(tl(i,j,k,t)+273.15)))
                end do
            end do
        end do
    end do
    open(12,file="/Users/wesleymoob/Desktop/se.grd",form='unformatted',access='stream')
    write(12)se
    close(12)
    !计算涡度、温度平流
    do t=1,nt
        do k=1,nz
            do j=1,ny
                do i=1,nx
                    if(i==1 .and. j==1)then
                        dty=(temp(i,j+1,k,t)-temp(i,j,k,t))/lon
                        dtx=(temp(i+1,j,k,t)-temp(i,j,k,t))/lat(j)
                        dwy=(vor(i,j+1,k,t)-vor(i,j,k,t))/lon
                        dwx=(vor(i+1,j,k,t)-vor(i,j,k,t))/lat(j)
                    else if(i==nx.and.j==ny)then
                        dty=(temp(i,j,k,t)-temp(i,j-1,k,t))/lon
                        dtx=(temp(i,j,k,t)-temp(i-1,j,k,t))/lat(j)
                        dwy=(vor(i,j,k,t)-vor(i,j-1,k,t))/lon
                        dwx=(vor(i,j,k,t)-vor(i-1,j,k,t))/lat(j)
                    else if(i==1 .and. j==ny)then
                        dty=(temp(i,j,k,t)-temp(i,j-1,k,t))/lon
                        dtx=(temp(i+1,j,k,t)-temp(i,j,k,t))/lat(j)
                        dwy=(vor(i,j,k,t)-vor(i,j-1,k,t))/lon
                        dwx=(vor(i+1,j,k,t)-vor(i,j,k,t))/lat(j)
                    else if(i==nx.and.j==1)then
                        dty=(temp(i,j+1,k,t)-temp(i,j,k,t))/lon
                        dtx=(temp(i,j,k,t)-temp(i-1,j,k,t))/lat(j)
                        dwy=(vor(i,j+1,k,t)-vor(i,j,k,t))/lon
                        dwx=(vor(i,j,k,t)-vor(i-1,j,k,t))/lat(j)
                    else if(i==1 .and. j/=1 .and. j/=ny)then
                        dty=(temp(i,j+1,k,t)-temp(i,j-1,k,t))/(2*lon)
                        dtx=(temp(i+1,j,k,t)-temp(i,j,k,t))/lat(j)
                        dwy=(vor(i,j+1,k,t)-vor(i,j-1,k,t))/(2*lon)
                        dwx=(vor(i+1,j,k,t)-vor(i,j,k,t))/lat(j)
                    else if(i==nx.and.j/=1.and.j/=ny)then
                        dty=(temp(i,j+1,k,t)-temp(i,j-1,k,t))/(2*lon)
                        dtx=(temp(i,j,k,t)-temp(i-1,j,k,t))/lat(j)
                        dwy=(vor(i,j+1,k,t)-vor(i,j-1,k,t))/(2*lon)
                        dwx=(vor(i,j,k,t)-vor(i-1,j,k,t))/lat(j)
                    else if(j==1.and.i/=1.and.i/=nx)then
                        dty=(temp(i,j+1,k,t)-temp(i,j,k,t))/lon
                        dtx=(temp(i+1,j,k,t)-temp(i-1,j,k,t))/(2.*lat(j))
                        dwy=(vor(i,j+1,k,t)-vor(i,j,k,t))/lon
                        dwx=(vor(i+1,j,k,t)-vor(i-1,j,k,t))/(2.*lat(j))
                    else if(j==ny.and.i/=1.and.i/=nx)then
                        dty=(temp(i,j,k,t)-temp(i,j-1,k,t))/lon
                        dtx=(temp(i+1,j,k,t)-temp(i-1,j,k,t))/(2.*lat(j))
                        dwy=(vor(i,j,k,t)-vor(i,j-1,k,t))/lon
                        dwx=(vor(i+1,j,k,t)-vor(i-1,j,k,t))/(2.*lat(j))
                    else
                        dty=(temp(i,j+1,k,t)-temp(i,j-1,k,t))/(2.*lon)
                        dtx=(temp(i+1,j,k,t)-temp(i-1,j,k,t))/(2.*lat(j))
                        dwy=(vor(i,j+1,k,t)-vor(i,j-1,k,t))/(2.*lon)
                        dwx=(vor(i+1,j,k,t)-vor(i-1,j,k,t))/(2.*lat(j))
                    end if
                    tpl(i,j,k,t)=-(u(i,j,k,t)*dtx(i,j,k,t)+v(i,j,k,t)*dty(i,j,k,t))
                    wdpl(i,j,k,t)=-(u(i,j,k,t)*dwx(i,j,k,t)+v(i,j,k,t)*dwy(i,j,k,t))
                end do
            end do
        end do
    end do
    open(90,file="/Users/wesleymoob/Desktop/tpl.grd",form='unformatted',access='stream')
    write(90)tpl
    close(90)
    open(100,file="/Users/wesleymoob/Desktop/wdpl.grd",form='unformatted',access='stream')
    write(100)wdpl
    close(100)
end program name

主要采用了有限差分等方法，其中较为复杂的部分在于计算相当位温。

运算结果使用grads制图，如上所示。

luna2222

想请问一下数据来源！谢谢！（ps:如果方便的话能否发到我的邮箱1217759768.qq.com里 谢谢！！！

WesleyMoob

luna2222 发表于 2022-11-18 09:08
想请问一下数据来源！谢谢！（ps:如果方便的话能否发到我的邮箱1217759768.qq.com里 谢谢！！！

数据来源于学校专业课程的实习，可能不太方便发诶

luna2222

WesleyMoob 发表于 2022-11-20 10:21
数据来源于学校专业课程的实习，可能不太方便发诶

好叭谢谢 没事嘿嘿

234kk

感谢大佬分享