cesm 边界层高度计算方法改进对理查森数判定对感热修正对f修正-CSDN博客

本文链接：https://blog.csdn.net/mayubins/article/details/126745988

本文探讨了如何修改理查森数的判定标准，引入感热修正，并优化f值处理，以提高大气湍流模型的准确性。着重于1月1日至5日的边界层高度模拟，新方法与旧算法进行了对比分析。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

修改1对理查森数的判定

首先介绍一下理查森数，

$N^2 = \frac{g}{\bar{\theta}}{}\frac{\partial \theta}{\partial z}$

Ri表示浮力项与流剪切项的比值的无量纲数。在物理学上，理查森数用来表示势能和动能的比值,在物理海洋学中，理查森数被用来研究海洋湍流，海洋混合。在大气上，理查森数表示大气静力稳定度与垂直风切变的比值。

现在根据RI来判别大气是不是稳定

一般来说，当RI>0.25时，大气是动力稳定的

经过试验，发现一般来说，大气的上部ri都是大于零的，高度越高这个数值越大，

在低层有时会出现ri<0的情况吗，这就是低层大气不稳定。

我们通过奥布霍夫和理查森数一起进行判断

当奥布霍夫长度大于零，理查森数大于0.25 ，则大气是稳定的，使用稳定边界层公式。

修改2 对感热修正

从上面的分析，我们知道当F*变小时，h会变大

F*是感热通量 $F* = \rho *Cp* \overline{ w' \theta' }$

如果除去 Cp则会将温度通量代入，这样会显著增加h

但是之前的分析认为，代入温度通量的边界层高度h 太高了，这里就回到原先的公式，直接将感热通量F*代入计算。

修改3 对f修正

由于科氏项在分母位，在低纬地区，由于科氏参数很小，则会容易造成计算不稳定。这需要对函数形式进行变形，防止在赤道附近出现计算不稳定，或者计算结果远远偏离实际。

那么最简单的方法就是当纬度 $\varphi$ 低于一定的纬度 $\varphi_0$ 时，则认为是这个维度 $\varphi_0$ .

当然南半球的纬度时负的这个需要注意到。我这里采用的方法是当低于15 度时则认为是15度

版本控制v3.1

 subroutine calc_pbl_h_vector_no_iterate(ri,h_ncol_pver,pbl_h_old,obklen,latvals,pver,n,T_c,tau_star,hsb_star,g_in,n2,PBL_H)
        !--------------------------------------!
        !purpose:calculate the height of PBL   !
        !mothod from the zili(2005)  paper     !
        !                                      !
        !mayubin 2022-5-19                     !
        !--------------------------------------!

        !=================!
        !  input variable !
        !=================!
        
        integer, intent(in)  :: n
        integer, intent(in)  :: pver
        real(r8),intent(in)  :: pbl_h_old(n)
        real(r8),intent(in)  :: obklen(n)
        real(r8),intent(in)  :: g_in  
        real(r8),intent(in)  :: T_c(n)
        real(r8),intent(in)  :: tau_star(n) 
        real(r8),intent(in)  :: hsb_star(n)      
        real(r8),intent(in)  :: n2(n,pver),ri(n,pver)
        real(r8),intent(in)  :: latvals(n)     !lat in radians (rad)
        real(r8),intent(in)  :: h_ncol_pver(n,pver)
        !=================!
        !  local variable !
        !=================!

        real(r8)             :: C_R =0.6_r8
        real(r8)             :: C_CN =1.36_r8
        real(r8)             :: C_NS =0.51_r8
        real(r8)             :: NN=0.01_r8     !ping lv
        real(r8)             :: f(n)           !defult 0.00005_r8 
        real(r8)             :: h_1(n)
        real(r8)             :: h_2(n)
        real(r8)             :: h_3(n)
        real(r8)             :: beta_b(n)
        real(r8)             :: h_pver_sum 
        integer              :: i
        integer              :: j
        integer              :: k
        real(r8)             :: n2_vertical_mean(n)
        real(r8)             :: pbl_h_stable(n)
        integer              :: diedaiindex =2
        integer,parameter    :: zuidadiedaicishu = 10 !die dai zui da ci shu 
        real(r8)             :: pbl_h_array(n,zuidadiedaicishu)
        integer              :: pingjuncengshu = 4
        real(r8)             :: h_pver(pver)
        integer              :: deltamin,deltamin_index
        real(r8)             :: delta(pver) 
        real(r8)             :: hsb_star_nocp(n)
        real(r8)             :: h_2_temp
        real(r8)             :: pbl_h_init(n)
        real(r8)             :: tau_star_2(n)
        real(r8)             :: f_min
        !==================!
        !  output variable !
        !==================!

        real(r8),intent(out) :: PBL_H(n)

        !===================!
        ! ji suan guo cheng !
        !===================!

        !!!!!====zheli xian ceshi yixia nayi ceng suanchulaide zuizhengchang ==
        ! brunt vaisaila frequency chuizhifangxiang pingjun    
        do i=1,n
                n2_vertical_mean(i) = 0
                do k =24,28
                        n2_vertical_mean(i) = n2_vertical_mean(i)+n2(i,k)
                end do 
                n2_vertical_mean(i) = n2_vertical_mean(i)/5
        end do 
do i = 1,n

        if(latvals(i)<0.and.latvals(i)>-15*3.14/180)then
                f(i) = -15*3.14/180
        end if
        if(latvals(i) >= 0 .and. latvals(i) <15* 3.14/180)then
                f(i) = 15*3.14/180
        end if

end do


        f=2_r8*7.292_r8*0.00001_r8*sin(latvals)
        

        beta_b = g_in/T_c
  
        ! hsb_star  ke neng duo cheng le yi ge C_p ,yao bu yao chu diao ,zheyang hui
        ! xian zhu zeng jia bian jie ceng gao du  
        !tau_star he u_star
        hsb_star_nocp  = hsb_star 

        h_1 = (f**2)/((C_R**2)*tau_star)
        h_2 = (sqrt(abs(n2_vertical_mean))*abs(f))/(tau_star*(C_CN**2))
        h_3 = abs(f*hsb_star_nocp*beta_b)/(C_CN**2*tau_star**2)
        pbl_h_stable = sqrt(1/(h_1+h_2+h_3))

!dui xinfangan de maxmin h jinxing xianzhi 

        do i = 1,n
                if(pbl_h_stable(i)>3000)then
                        pbl_h_stable(i) = 3000
                end if
                if(pbl_h_stable(i)<10)then 
                        pbl_h_stable = 10
                end if  
        end do



        !dui wen ding du jin xing pan duan 

        do  i = 1, n
                
                if ( obklen(i) > 0 .and. ri(i,31)>0.25_r8) then

                        PBL_H(i) = pbl_h_stable(i) 
                else
                        PBL_H(i) = pbl_h_old(i)
                end if
        end do

      end subroutine calc_pbl_h_vector_no_iterate