关于cesm中寻找浮力频率N的记录

前几天模拟器调试，一直是登不上的。今天终于可以登上了，我来寻找N和f了，把这两个加上我的h公式计算就完整了。

首先我们看到pbl_utils.f90 文件的前面一个计算子程序austausch_atm

咦我们看到什么好东西了？

这不是理查森数和风剪切吗？

那么浮力频率不就可以计算出来了？

那么我们看看那些程序调用了它

确实找到了两个程序

分别在不同的文件里

我们进去看看是怎么调用的

 

我们看到，这两个变量分别就是ri 和s2

看看定义 

没有找到，应该是写在模块里了 

这样，我在hb_diff,F90文件里计算一下ri 和s2 ,然后在vertical_diffusion.F90里面调用一下就可以 

实际上我可能考虑麻烦了

vertical_diffusion.F90里面有ri ,再看看有没有s2甚至N.

真的没有

我正准备在hb_diff里写一个计算n2的子程序

突然

我发现他本来就有，这样我就直接调用就好了。 

既然vertical_diffusion.f90有s2，ri ,那肯定有N^2了

然后我再改一下我的子程序，

然后编译一下看看有没有问题

出现了一些低级错误

 

少打一个逗号，

加上再试一下

我估计还是会有一些错误，因为数组维度不一样

 

 又一个低级错误，没定义n2

定义一下重新

但是我隐约觉得不对

因为二维数组呀，传到一位数组了

我觉得我的整个子程序都得改了

一旦涉及N和f 就是二维平面的问题了，我裂开。

果然，似乎说我的数组不匹配了。

在n2 后面加（：ncol）这样不知道行不行。不行的话真是难搞呀

实际上n2 是一个二维数组

那么我需要对垂直方向取平均才可以变其为一位数组，而且我本来也是需要这么计算的

 我们先把数组维数对上，再来考虑是垂直的几层进行平均

我直接到hb_diff 里改trbintd函数

 

增加了一个对垂直方向平均的n2

我就怕其他文件应用的时候出问题。

 这里有两个地方

分别改一下吧

我还是觉得这样改不妥，计算平均的可以单独放一个程序里就好了

这样增加了调用出错的可能性

最终，我是输入一个二维数组到子程序里，由子程序计算一下浮力频率的垂直平均。

 subroutine calc_pbl_h_vector(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     !
        !                                      !
        !--------------------------------------!
        integer, intent(in)  :: n
        integer, intent(in)  :: pver
        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)
        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 !pinglv
        real(r8)  :: f=0.00005_r8 
        real(r8)  :: h_1(n)
        real(r8)  :: h_2(n)
        real(r8)  :: h_3(n)
        real(r8)  :: beta_b(n)
        integer   :: i
        integer   :: k
!zhihou  hui tian jia f he N de jisuan 

        real(r8)  :: n2_vertical_mean(n) ! brunt vaisaila frequency chuizhifangxiang pingjun    
!
        real(r8),intent(out)  :: PBL_H(n)
              !mayubin 2022-3-25 
         do i=1,n
                n2_vertical_mean(i) = 0
                do k =1,pver
                   n2_vertical_mean(i) = n2_vertical_mean(i)+n2(i,k)
                end do 
                n2_vertical_mean(i) = n2_vertical_mean(i)/pver
        end do 

        beta_b = g_in/T_c

        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*beta_b)/(C_CN**2*tau_star**2)

        PBL_H = sqrt(1/(h_1+h_2+h_3))
       
      end subroutine calc_pbl_h_vector

这样计算出的pbl_h 我觉得还是存在一些问题，有些地区极值很大，而且梯度很大，就好像很多奇异点，这画出来的图我都看着头皮发麻

 我觉得和这三项之间的计算转换存在一定的关系。关于这个就是计算方案的问题了，而不是技术上的问题了。

针对这个问题，我把第三项去掉

 我们看到，边界层高度变得很大，这可能和系数有关系，我应该分别计算三项分别的高度。

就在我考察这个公式的时候我突然发现一个问题，那就是这个方程的第三项的量纲似乎是不对的，这让我感到非常高兴，因为我又可以修正问题了。

之后我会将这个问题解决然后提出我的新计算公式重新带入到模式里面进行计算。

今天是周六，我成功将N代入到计算公式里面，又意外发现了文章里公式的重大错误，今天是重要的一天。

但是f的作用还没加进去，这个会在后面加进去。

对于全球来说，f（y）的变化是非常明显的，这样边界层高度便会随着纬度而变化。