はじめに

無次元化や座標のとり方はこのあたりと同じです。
https://qiita.com/ishigaki/items/866df91ea8f32136fa5a
https://qiita.com/ishigaki/items/60bb3c7fc9a66d278cfa
PISO法の概要についてはVersteeg "An Introduction to Computational Fluid Dynamics" http://ftp.demec.ufpr.br/disciplinas/TM702/Versteeg_Malalasekera_2ed.pdf の教科書を参考にしました。PISO法のベースとなるSIMPLE法については岡本「数値計算による流体力学- ポテンシャル流,層流,そして乱流へ -」を参考にしました。
コードのここおかしくないとかあったらコメントください。
コードの概要

Paramで計算につかうパラメータを設定してます。
mainで変数の宣言をして，初期化init, 境界条件BCを呼び出す。
solveUA, solveVAで1段目の予測速度u*, v*を計算。
solvePCで1段目の圧力修正量p'の計算をします。
updateで2段めの予測速度u**, v**を計算します。
solvePC2で2段めの圧力修正量p''を計算します。
update2で3段目の予測速度u***, v***を計算します。
solveTAで温度を計算します。
mainのfor n = 1:niterのループ中で速度が収束したら，ループを抜けます。次の時間ステップへ。
計算リスタートするときは2個目のmainを呼びます。
対流項は2次精度の中心差分（線形補間），拡散項は2次精度の中心差分で離散化。時間微分はEuler陰解法で離散化します。
下面が高温，上面が低温，側面は断熱。速度については滑りなし。
コード

using Parameters
using HDF5

@with_kw struct Param{T1,T2}
    dt::T1 = 0.01

    # Rayleigh number, Prandtl number, Grashoff number
    Ra::T1 = 1e4
    Pr::T1 = 1.0
    Gr::T1 = Ra/Pr
    nu::T1 = 1.0/sqrt(Gr)
    alpha::T1 = 1.0/Pr/sqrt(Gr)

    # boundary temp
    Th::T1 = 1.0
    Tc::T1 = 0.0

    nstep::T2 = 200
    niter::T2 = 500
    ipter::T2 = 100
    tol::T1 = 1e-5
    
    # mesh size
    nx::T2 = 40
    ny::T2 = 40
    
    # length
    Lx::T1 = 1.0
    Ly::T1 = 1.0
    dx::T1 = Lx/nx
    dy::T1 = Ly/ny

    ndata::T2 = 10

    outfile::String = "Data-piso/uv" 
end

function main(para)
    @unpack nx, ny, nstep, niter, ndata, tol = para
    
    ua = zeros(Float64, nx+2, ny+2)
    ua2 = zeros(Float64, nx+2, ny+2)
    ua3 = zeros(Float64, nx+2, ny+2)
    utemp = zeros(Float64, nx+2, ny+2)
    uold = zeros(Float64, nx+2, ny+2)
    apua = zeros(Float64, nx+2, ny+2)
    aeua = zeros(Float64, nx+2, ny+2)
    awua  = zeros(Float64, nx+2, ny+2)
    anua = zeros(Float64, nx+2, ny+2)
    asua = zeros(Float64, nx+2, ny+2)
    bua  = zeros(Float64, nx+2, ny+2)
    
    va = zeros(Float64, nx+2, ny+2)
    va2 = zeros(Float64, nx+2, ny+2)
    va3 = zeros(Float64, nx+2, ny+2)
    vtemp = zeros(Float64, nx+2, ny+2)
    vold = zeros(Float64, nx+2, ny+2)
    apva  = zeros(Float64, nx+2, ny+2)
    aeva = zeros(Float64, nx+2, ny+2)
    awva = zeros(Float64, nx+2, ny+2)
    anva = zeros(Float64, nx+2, ny+2)
    asva = zeros(Float64, nx+2, ny+2)
    bva = zeros(Float64, nx+2, ny+2)
    
    # temperature
    Ta = zeros(Float64, nx+2, ny+2)
    Told = zeros(Float64, nx+2, ny+2)
    apTa = zeros(Float64, nx+2, ny+2)
    aeTa = zeros(Float64, nx+2, ny+2)
    awTa = zeros(Float64, nx+2, ny+2)
    anTa = zeros(Float64, nx+2, ny+2)
    asTa = zeros(Float64, nx+2, ny+2)
    bTa = zeros(Float64, nx+2, ny+2)
    Ttemp = zeros(Float64, nx+2, ny+2)

    ue = zeros(Float64, nx+2, ny+2)
    uw = zeros(Float64, nx+2, ny+2)
    vn = zeros(Float64, nx+2, ny+2)
    vs = zeros(Float64, nx+2, ny+2)

    # pressure, corrected pressure
    pa = zeros(Float64, nx+2, ny+2)
    pc = zeros(Float64, nx+2, ny+2)
    ptemp = zeros(Float64, nx+2, ny+2)
    pc2 = zeros(Float64, nx+2, ny+2)
    appc = zeros(Float64, nx+2, ny+2)
    aepc = zeros(Float64, nx+2, ny+2)
    awpc = zeros(Float64, nx+2, ny+2)
    anpc = zeros(Float64, nx+2, ny+2)
    aspc = zeros(Float64, nx+2, ny+2)
    bpc = zeros(Float64, nx+2, ny+2)

    # coordinate
    x = zeros(Float64, nx+1, ny+1)
    y = zeros(Float64, nx+1, ny+1)
    
    rest = 0.0
    init(para, x, y)

    for l = 1:nstep
        BC(para, ua, va, Ta, ua2, va2)
        
        for n = 1:niter
            solveUA(para, ua, va, pa, apua, aeua, awua, anua, asua, bua, uold, utemp)
            solveVA(para, ua, va, Ta, pa, apva, aeva, awva, anva, asva, bva, vold, vtemp)
            solvePC(para, ua, va, pc, ptemp, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc)
            update(para, ua, va, pc, apua, apva, ua2, va2)
            solvePC2(para, ua, ua2, va, va2, pc, pc2, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc, ue, uw, vn, vs, ptemp)
            update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apua, aeua, awua, anua, asua, apva, aeva, awva, anva, asva)
            solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
            
            rest = converg(para, ua, va, utemp, vtemp)
            #println("iter, redisual ", n,", ",rest)
            
            if rest <=tol && n >=2
                if l%10 == 0
                    println(l, " th step, ", n, " th iteration, residual=", rest)
                end
                break
            end
            println(rest)
        end
        for j=2:ny+1
            for i=2:nx
                uold[i,j] = ua[i,j]
            end
        end
    
        for j=2:ny
            for i=2:nx+1
                vold[i,j] = va[i,j]
            end
        end
        
        for j=2:ny+1
            for i=2:nx+1
                Told[i,j] = Ta[i,j]
            end
        end
        
        if l%ndata == 0
            println("residual=", rest)
            output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
        end
    end
    
end

function main(para, restartFile)
    @unpack nx, ny, nstep, niter, ndata, tol = para
    # auxiliary velocity
    ua = zeros(Float64, nx+2, ny+2)
    ua2 = zeros(Float64, nx+2, ny+2)
    ua3 = zeros(Float64, nx+2, ny+2)
    utemp = zeros(Float64, nx+2, ny+2)
    uold = zeros(Float64, nx+2, ny+2)
    apua = zeros(Float64, nx+2, ny+2)
    aeua = zeros(Float64, nx+2, ny+2)
    awua  = zeros(Float64, nx+2, ny+2)
    anua = zeros(Float64, nx+2, ny+2)
    asua = zeros(Float64, nx+2, ny+2)
    bua  = zeros(Float64, nx+2, ny+2)
    
    va = zeros(Float64, nx+2, ny+2)
    va2 = zeros(Float64, nx+2, ny+2)
    va3 = zeros(Float64, nx+2, ny+2)
    vtemp = zeros(Float64, nx+2, ny+2)
    vold = zeros(Float64, nx+2, ny+2)
    apva  = zeros(Float64, nx+2, ny+2)
    aeva = zeros(Float64, nx+2, ny+2)
    awva = zeros(Float64, nx+2, ny+2)
    anva = zeros(Float64, nx+2, ny+2)
    asva = zeros(Float64, nx+2, ny+2)
    bva = zeros(Float64, nx+2, ny+2)
    
    # temperature
    Ta = zeros(Float64, nx+2, ny+2)
    Told = zeros(Float64, nx+2, ny+2)
    apTa = zeros(Float64, nx+2, ny+2)
    aeTa = zeros(Float64, nx+2, ny+2)
    awTa = zeros(Float64, nx+2, ny+2)
    anTa = zeros(Float64, nx+2, ny+2)
    asTa = zeros(Float64, nx+2, ny+2)
    bTa = zeros(Float64, nx+2, ny+2)
    Ttemp = zeros(Float64, nx+2, ny+2)

    ue = zeros(Float64, nx+2, ny+2)
    uw = zeros(Float64, nx+2, ny+2)
    vn = zeros(Float64, nx+2, ny+2)
    vs = zeros(Float64, nx+2, ny+2)

    # pressure, corrected pressure
    pa = zeros(Float64, nx+2, ny+2)
    pc = zeros(Float64, nx+2, ny+2)
    ptemp = zeros(Float64, nx+2, ny+2)
    pc2 = zeros(Float64, nx+2, ny+2)
    appc = zeros(Float64, nx+2, ny+2)
    aepc = zeros(Float64, nx+2, ny+2)
    awpc = zeros(Float64, nx+2, ny+2)
    anpc = zeros(Float64, nx+2, ny+2)
    aspc = zeros(Float64, nx+2, ny+2)
    bpc = zeros(Float64, nx+2, ny+2)

    # coordinate
    x = zeros(Float64, nx+1, ny+1)
    y = zeros(Float64, nx+1, ny+1)
    
    rest = 0.0
    
    nstart = restart(para, restartFile, x, y, ua, va, pa, Ta, uold, vold, Told)
    println("nstart=",nstart)

    for l = nstart+1:nstart+nstep
        BC(para, ua, va, Ta, ua2, va2)
        
        for n = 1:niter
            solveUA(para, ua, va, pa, apua, aeua, awua, anua, asua, bua, uold, utemp)
            solveVA(para, ua, va, Ta, pa, apva, aeva, awva, anva, asva, bva, vold, vtemp)
            solvePC(para, ua, va, pc, ptemp, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc)
            update(para, ua, va, pc, apua, apva, ua2, va2)
            solvePC2(para, ua, ua2, va, va2, pc, pc2, apua, aeua, awua, anua, asua, bua, apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc, ue, uw, vn, vs, ptemp)
            update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apua, aeua, awua, anua, asua, apva, aeva, awva, anva, asva)
            solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
            
            rest = converg(para, ua, va, utemp, vtemp)
            #println("iter, redisual ", n,", ",rest)
            
            if rest <=tol && n >=2
                if l%10 == 0
                    println(l, " th step, ", n, " th iteration, residual=", rest)
                end
                break
            end
            println(rest)
        end
        for j=2:ny+1
            for i=2:nx
                uold[i,j] = ua[i,j]
            end
        end
    
        for j=2:ny
            for i=2:nx+1
                vold[i,j] = va[i,j]
            end
        end
        
        for j=2:ny+1
            for i=2:nx+1
                Told[i,j] = Ta[i,j]
            end
        end
        
        if l%ndata == 0
            println("residual=", rest)
            output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
        end
    end
    
end

function init(para, x, y)
    @unpack nx, ny, dx, dy = para
    
    for j in 2:ny+1
        for i in 2:nx+1
            x[i,j] = (i-2+0.5)*dx
        end
    end
    for j in 2:ny+1
        for i in 2:nx+1
            y[i,j] = (j-2+0.5)*dy
        end
    end
end

function BC(para, ua, va, Ta, ua2, va2)
    @unpack nx, ny, Th, Tc = para
    for i=1:nx+2
        ua[i, ny+2] = 0.0
        va[i, ny+1] = 0.0
        ua[i, 1] = 0.0
        va[i, 1] = 0.0
        Ta[i, ny+2] = Tc
        Ta[i, 1] = Th
        # ua2[i, ny+2]をupdate2でつかうが，updateではここの値は更新されないので，BCで与える。
        ua2[i, ny+2] = 0.0
        va2[i, ny+1] = 0.0
        ua2[i, 1] = 0.0
        va2[i, 1] = 0.0
    end
    
    for j=1:ny+2
        ua[nx+1, j] = 0.0
        va[nx+2, j] = 0.0
        ua[1, j] = 0.0
        va[1, j] = 0.0
        Ta[nx+2, j] = Ta[nx+1, j]
        Ta[1, j] = Ta[2,j]
        ua2[nx+1, j] = 0.0
        va2[nx+2, j] = 0.0
        ua2[1, j] = 0.0
        va2[1, j] = 0.0
    end
end

function solveUA(para, ua, va, pa, apua, aeua, awua, anua, asua, bua, uold, utemp)
    @unpack nx, ny, dx, dy, dt, nu,  = para
    
    for j=2:ny+1
        for i=1:nx+1
            if i==nx+1
                aeua[i,j] = nu*dy/dx
            else
                aeua[i,j] = nu*dy/dx - 0.25*dy*(ua[i+1, j] + ua[i,j])
            end
            if i==1
                awua[i,j] = nu*dy/dx
            else
                awua[i,j] = nu*dy/dx + 0.25*dy*(ua[i, j] + ua[i-1,j])
            end

            if j==2
                asua[i,j] = 2nu*dx/dy
            else
                asua[i,j] = nu*dx/dy + 0.25dx*(va[i+1,j-1] + va[i,j-1])
            end
            if j == ny+1
                anua[i,j] = 2nu*dx/dy
            else
                anua[i,j] = nu*dx/dy - 0.25dx*(va[i+1,j] + va[i,j])
            end
            
            apua[i,j] = aeua[i,j] + awua[i,j] + anua[i,j] + asua[i,j] + dx*dy/dt
            bua[i,j] = -dy*(pa[i+1,j] - pa[i,j]) +dx*dy/dt*uold[i,j]
            
            if i!=1 && i!=nx+1
                utemp[i,j] = (aeua[i,j]*ua[i+1,j] + awua[i,j]*ua[i-1,j] + anua[i,j]*ua[i,j+1] + asua[i,j]*ua[i,j-1] + bua[i,j])/apua[i,j]
            end
        end
    end

    for j=2:ny+1
        for i=2:nx
            ua[i,j] = utemp[i,j]
        end
    end
end

function solveVA(para, ua, va, Ta, pa, apva, aeva, awva, anva, asva, bva, vold, vtemp)
    @unpack nx, ny, dx, dy, dt, nu,  = para
    
    for j=1:ny+1
        for i=2:nx+1
            if i==2
                awva[i,j] = 2nu*dy/dx
            else
                awva[i,j] = nu*dy/dx + 0.25dy*(ua[i-1,j+1] + ua[i-1,j])
            end
            if i == nx+1
                aeva[i,j] = 2nu*dy/dx
            else
                aeva[i,j] = nu*dy/dx - 0.25dy*(ua[i,j+1] + ua[i,j])
            end

            if j==ny+1
                anva[i,j] = nu*dx/dy
            else
                anva[i,j] = nu*dx/dy - 0.25*dx*(va[i, j+1] + va[i,j])
            end
            if j==1
                asva[i,j] = nu*dx/dy
            else
                asva[i,j] = nu*dx/dy + 0.25*dx*(va[i, j] + va[i,j-1])
            end
            
            apva[i,j] = aeva[i,j] + awva[i,j] + anva[i,j] + asva[i,j] + dx*dy/dt
            bva[i,j] = -dx*(pa[i,j+1] - pa[i,j]) + dx*dy/dt*vold[i,j] + dx*dy*0.5*(Ta[i, j]+Ta[i, j+1])
            
            if j!=1 && j!=ny+1
                vtemp[i,j] = (aeva[i,j]*va[i+1,j] + awva[i,j]*va[i-1,j] + anva[i,j]*va[i,j+1] + asva[i,j]*va[i,j-1] + bva[i,j])/apva[i,j]
            end
        end
    end
    for j=2:ny
        for i=2:nx+1
            va[i,j] = vtemp[i,j]
        end
    end
end

function solvePC(para, ua, va, pc, ptemp, apua, aeua, awua, anua, asua, bua, 
        apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc)
    @unpack nx, ny, dx, dy, nu, ipter = para
    
    for j=2:ny+1
        for i=2:nx+1
            pc[i,j] = 0.0
            bpc[i,j] = 0.0
        end
    end
    
    for j=2:ny+1
        for i=2:nx+1
            if i == 2
                awpc[i,j] = 0.0
            else
                awpc[i,j] = dy*dy/apua[i-1,j]
            end
            
            if i == nx+1
                aepc[i,j] = 0.0
            else          
                aepc[i,j] = dy*dy/apua[i,j]
            end
            
            if j == 2
                aspc[i,j] = 0.0
            else
                aspc[i,j] = dx*dx/apva[i,j-1]
            end
                    
            if j == ny+1
                anpc[i,j] = 0.0
            else
                anpc[i,j] = dx*dx/apva[i,j]
            end
                        
            appc[i,j] = aepc[i,j] + awpc[i,j] + anpc[i,j] + aspc[i,j]
            bpc[i,j] = -dy*(ua[i,j] - ua[i-1,j]) - dx*(va[i,j] - va[i,j-1])
        end
    end
    
    # SOR法で圧力の解を求めてる
    for ip=1:ipter
        for j=2:ny+1
            for i=2:nx+1
                ptemp[i,j] = pc[i,j]
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                pc[i,j] = (aepc[i,j]*ptemp[i+1,j] + awpc[i,j]*ptemp[i-1,j] + anpc[i,j]*ptemp[i,j+1] + aspc[i,j]*ptemp[i,j-1] + bpc[i,j])/appc[i,j]
            end
        end
        rest=convergP(para, pc, ptemp)
        
        if rest < 1e-5 && ip >2
            #println("rest PC ", rest, ", ip=", ip)
            break
        end
        
    end
end    

function update(para, ua, va, pc, apua, apva, ua2, va2)
    @unpack nx, ny, dx, dy,  = para
    
    # u** = u* + u'
    for j=2:ny+1
        for i=2:nx
            ua2[i,j] = ua[i,j] -dy*(pc[i+1,j] - pc[i,j])/apua[i,j] 
        end
    end
    
    for j=2:ny
        for i=2:nx+1
            va2[i,j] = va[i,j] -dx*(pc[i,j+1] - pc[i,j])/apva[i,j]
        end
    end
end

function solvePC2(para, ua, ua2, va, va2, pc, pc2, apua, aeua, awua, anua, asua, bua, 
    apva, aeva, awva, anva, asva, bva, appc, aepc, awpc, anpc, aspc, bpc, ue, uw, vn, vs, ptemp)
    @unpack nx, ny, dx, dy, nu, ipter = para

    for j=2:ny+1
        for i=2:nx+1
            pc2[i,j] = 0.0
            bpc[i,j] = 0.0
        end
    end

    for j=2:ny+1
        for i=2:nx+1

            if i==nx+1
                ue[i,j] = awua[i,j]*(ua2[i-1,j]-ua[i-1,j])                          
            else
                ue[i,j] = aeua[i,j]*(ua2[i+1,j]-ua[i+1,j]) +awua[i,j]*(ua2[i-1,j]-ua[i-1,j]) +anua[i,j]*(ua2[i,j+1]-ua[i,j+1]) +asua[i,j]*(ua2[i,j-1]-ua[i,j-1])
            end

            if i==2
                uw[i,j] = aeua[i-1,j]*(ua2[i,j]-ua[i,j]) 
            else
                uw[i,j] = aeua[i-1,j]*(ua2[i,j]-ua[i,j]) +awua[i-1,j]*(ua2[i-2,j]-ua[i-2,j]) +anua[i-1,j]*(ua2[i-1,j+1]-ua[i-1,j+1]) +asua[i-1,j]*(ua2[i-1,j-1]-ua[i-1,j-1])
            end

            if j==ny+1
                vn[i,j] = asva[i,j]*(va2[i,j-1]-va[i,j-1])
            else
                vn[i,j] = aeva[i,j]*(va2[i+1,j]-va[i+1,j]) +awva[i,j]*(va2[i-1,j]-va[i-1,j]) +anva[i,j]*(va2[i,j+1]-va[i,j+1]) +asva[i,j]*(va2[i,j-1]-va[i,j-1])
            end

            if j==2
                vs[i,j] = anva[i,j-1]*(va2[i,j]-va[i,j])
            else
                vs[i,j] = aeva[i,j-1]*(va2[i+1,j-1]-va[i+1,j-1]) +awva[i,j-1]*(va2[i-1,j-1]-va[i-1,j-1]) +anva[i,j-1]*(va2[i,j]-va[i,j]) +asva[i,j-1]*(va2[i,j-2]-va[i,j-2])
            end
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            bpc[i,j] = -dy/apua[i,j]*ue[i,j] +dy/apua[i-1,j]*uw[i,j] -dx/apva[i,j]*vn[i,j] +dx/apva[i,j-1]*vs[i,j]           
        end
    end

    # SOR法で圧力の解を求めてる
    for ip=1:ipter
        for j=2:ny+1
            for i=2:nx+1
                ptemp[i,j] = pc2[i,j]
            end
        end
        for j=2:ny+1
            for i=2:nx+1
                pc2[i,j] = (aepc[i,j]*ptemp[i+1,j] + awpc[i,j]*ptemp[i-1,j] + anpc[i,j]*ptemp[i,j+1] + aspc[i,j]*ptemp[i,j-1] + bpc[i,j])/appc[i,j]
            end
        end
        
        rest=convergP(para, pc2, ptemp)
        if rest < 1e-5
            #println("rest PC2 ", rest, ", ip=", ip)
            break
        end
    end
end

function update2(para, ua, ua2, ua3, va, va2, va3, pa, pc, pc2, apua, aeua, awua, anua, asua, apva, aeva, awva, anva, asva)
    @unpack nx, ny, dx, dy,  = para

    # u** = u* + u'
    for j=2:ny+1
        for i=2:nx
            ua3[i,j] = ua2[i,j] +
            (aeua[i,j]*(ua2[i+1,j]-ua[i+1,j]) + awua[i,j]*(ua2[i-1,j]-ua[i-1,j]) + anua[i,j]*(ua2[i,j+1]-ua[i,j+1]) + asua[i,j]*(ua2[i,j-1]-ua[i,j-1]) -
            dy*(pc2[i+1,j] - pc2[i,j]))/apua[i,j] 
        end
    end

    for j=2:ny+1
        for i=2:nx
            ua[i,j] = ua3[i,j]
        end
    end
    
    for j=2:ny
        for i=2:nx+1
            va3[i,j] = va2[i,j] +
            (aeva[i,j]*(va2[i+1,j]-va[i+1,j]) + awva[i,j]*(va2[i-1,j]-va[i-1,j]) + anva[i,j]*(va2[i,j+1]-va[i,j+1]) + asva[i,j]*(va2[i,j-1]-va[i,j-1]) -
            dx*(pc2[i,j+1] - pc2[i,j]))/apva[i,j]
        end
    end

    

    for j=2:ny
        for i=2:nx+1
            va[i,j] = va3[i,j]
        end
    end

    for j=2:ny+1
        for i=2:nx+1
            pa[i,j] = pa[i,j] + pc[i,j] + pc2[i,j]
        end
    end
    pref = pa[2,2]
    for j=2:ny+1
        for i=2:nx+1
            pa[i,j] -= pref
        end
    end
end

function solveTA(para, ua, va, Ta, apTa, aeTa, awTa, anTa, asTa, bTa, Told, Ttemp)
    @unpack nx, ny, dx, dy, dt, alpha,  = para
    
    for j=2:ny+1
        for i=2:nx+1
            if i==2
                awTa[i,j] = 2alpha*dy/dx
            else
                awTa[i,j] = alpha*dy/dx + 0.5dy*ua[i-1,j]
            end
            if i==nx+1
                aeTa[i,j]=2alpha*dy/dx
            else
                aeTa[i,j] = alpha*dy/dx - 0.5dy*ua[i,j]
            end

            if j == ny+1
                anTa[i,j] = 2alpha*dx/dy
            else
                anTa[i,j] = alpha*dx/dy - 0.5*dx*va[i,j]
            end
            if j == 2
                asTa[i,j] = 2alpha*dx/dy
            else
                asTa[i,j] = alpha*dx/dy + 0.5*dx*va[i,j-1]
            end
            
            apTa[i,j] = aeTa[i,j] + awTa[i,j] + anTa[i,j] + asTa[i,j] + dx*dy/dt
            bTa[i,j] = dx*dy/dt*Told[i,j]
            
            Ttemp[i,j] = (aeTa[i,j]*Ta[i+1,j] + awTa[i,j]*Ta[i-1,j] + anTa[i,j]*Ta[i,j+1] + asTa[i,j]*Ta[i,j-1] + bTa[i,j])/apTa[i,j]
        end
    end
    for j=2:ny+1
        for i=2:nx+1
            Ta[i,j] = Ttemp[i,j]
        end
    end 
end

function converg(para, ua, va, utemp, vtemp)
    @unpack nx, ny = para
    
    restua = 0.0
    restva = 0.0
    
    for j=2:ny+1
        for i=2:nx
            restua = restua + abs(ua[i,j] - utemp[i,j])
            
        end
    end
    for j=2:ny
        for i=2:nx+1
            restva = restva + abs(va[i,j] - vtemp[i,j])
        end
    end

    rest = (restua + restva)/2.0
    return rest/(nx*ny)
end

function convergP(para, pc, ptemp)
    @unpack nx, ny = para
    
    rest = 0.0
    
    for j=2:ny+1
        for i=2:nx+1
            rest = rest + abs(pc[i,j] - ptemp[i,j])
        end
    end

    
    return rest/(nx*ny)
end

function output(para, ua, va, pa, Ta, uold, vold, Told, x, y, l)
    @unpack outfile, dt, nx, ny = para
    #println("time=", time)
    filename = outfile*string(l)*".h5"
    h5open(filename, "w") do file
        write(file, "U", (ua[2:nx+1, 2:ny+1]+ua[1:nx, 2:ny+1])/2)
        write(file, "V", (va[2:nx+1, 2:ny+1]+va[2:nx+1, 1:ny])/2)
        write(file, "T", Ta[2:nx+1, 2:ny+1])
        write(file, "P", pa[2:nx+1, 2:ny+1])
        write(file, "X", x[2:nx+1, 2:ny+1])
        write(file, "Y", y[2:nx+1, 2:ny+1])
        write(file, "ua_latest", ua)
        write(file, "va_latest", va)
        write(file, "pa_latest", pa)
        write(file, "Ta_latest", Ta)
        write(file, "ua_old", uold)
        write(file, "va_old", vold)
        write(file, "Ta_old", Told)
        write(file, "n", l)
    end
end

function restart(para, restartFile, x, y, ua, va, pa, Ta, uold, vold, Told)
    @unpack nx, ny, dx, dy = para
    for j in 2:ny+1
        for i in 2:nx+1
            x[i,j] = (i-2+0.5)*dx
        end
    end
    for j in 2:ny+1
        for i in 2:nx+1
            y[i,j] = (j-2+0.5)*dy
        end
    end
    
    file = h5open(restartFile, "r") 
    ua .= read(file, "ua_latest")
    va .= read(file, "va_latest")
    pa .= read(file, "pa_latest")
    Ta .= read(file, "Ta_latest")
    uold .= read(file, "ua_old")
    vold .= read(file, "va_old")
    Told .= read(file, "Ta_old")
    nstart = read(file, "n")
    
    return nstart
end
計算結果

Ra=1e8, Pr=5.3, nx=ny=512, dt=2e-3の計算結果
SIMPLE法，PISO法の理屈

有限体積法による離散化
SIMPLE法
非定常SIMPLE
PISO
もしくは