WRF内部程序学习，第一弹：中间变量的输出

tianmeng · 发表于 2023-2-12 18:13:52

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

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

x

   随着研究深入，我们常常需要直接上手修改WRF子程序的一部分，为了评估修改效果，首先涉及输出一些WRF程序内部变量（local var）作为检验。    因此WRF内部程序学习的第一部分，应该是WRF中间变量的输出，尤其是wrfout默认输出中没有的中间变量的输出。

第一部分.如果是Registry.EM_COMMON注册表中已有但没有输出的变量，情况是非常简单的，只要设置该变量属性为h，即可输出即可：
      具体可参考气象家园贴：https://bbs.06climate.com/forum.php?mod=viewthread&tid=2811

   由于该方法需要重新编译WRF程序，一个更简单的方法是，在wrf.exe的运行目录下增加一个myfile.txt文件，写入注册表中已有但是没有输出的变量，这里给个例子，如果我想输出湍流动能，长短波辐射率和湍流交换系数：
   +:h:0:tke,RTHRATLW,RTHRATEN,RTHRATSW,EXCH_M,EXCH_H
   保存后修改namelist.input，在&time_control后加入：
   iofields_filename = "myfile.txt",
   ignore_iofields_warning = .true.,
   无需编译，即可在wrfout文件中新增变量，用同样方法，修改myfile.txt的“+”为“-”号，也可轻松去掉不想要的变量，从而缩小wrfout输出文件大小。

第二部分. 如果是Registry.EM_COMMON注册表中没有的变量（但是WRF内部程序有）或者自己新增的变量（内部程序也没有），情况会略复杂，但也能轻松搞定，参考WRF官网提供的Annual_Report_WRF_PHYSICS.pdf，可以找到关于WRF框架的介绍，即WRF程序包含三层结构：求解器（dyn_em/solver）、驱动程序（phys/driver）和独立方案（phys/module），这里简单介绍一下它们的功能：
（1）第一层：求解器solve_rk.F，WRF动力和物理模块之间桥梁，主要用于变量传递和并行。
（2）第二层：驱动程序（microphysics_driver, cumulus_driver, pbl_driver, and radiation_driver等），求解器不是直接调用不同的物理方案，而是调用物理驱动程序（driver），因此driver是求解器和单个物理方案之间的接口。除了子网格涡流扩散（sub-grid eddy diffusion），每个物理参数化方案都有自己的driver程序或共享driver程序，例如pbl_driver包括表面层、次表面层和边界层的驱动，radiation_driver包括长波和短波辐射的驱动。这两个模块是我们需要输出中间变量的重灾区。在每个driver中，CASE SELECT结构用于确定不同物理参数化方案的选择（对应namelis.input），这里我们也可以插入自己编写的参数化物理包（略，详见Annual_Report_WRF_PHYSICS）。
（3）第三层: 独立方案（module_*_*.F），这就没什么可介绍的了，找对自己需要使用和修改的方案就可以了，Annual_Report_WRF_PHYSICS有介绍命名规则。
   可以想见，如果我们需要输出一个注册表中不存在的变量，我们需要通过这三层结构，才能传递给注册表进行输出，这里我们从独立方案（3）开始改起，然后是驱动程序（2）和求解器（1），最后传递到注册表Registry.EM_COMMON，  如果我们非常熟悉整个程序，只需要两步，A确定需要输出的变量，B增加Fortran声明，自下而上三层挨个增加声明，最后加入注册表重新编译。

   如果不想看复杂的程序，但只要能确定需要输出的变量名，我们跟着一个模块可输出的变量，手慢的10分钟也可以轻松搞定，举个例子：

   假设我对YSU边界层方案计算流程一窍不通，而我想输出YSU方案的一个中间量体积理查孙数bulk Richardson number，因为它的计算能够确定边界层高度。如果后续我要利用梯度观测来修正这个数，那我首先得已知它在程序里叫（brup），在输出过程中，我首先随便跟着一个叫EXCH_H（热力湍流扩散系数）的量来写，因为我知道它在注册表里（见第一部分例子），而且它也在YSU方案里，并且它在独立方案（module_bl_ysu.F）中和我要输出的变量是同维度的：
   （1）那我先找到二维变量exch_hx，它的前身是xkzh，这里xkzh还是一个二维的计算维度（its:ite, kts:kte），当赋值给exch_hx的时候，exch_hx的声明就变成存储维度（ims:ime, kms:kme）。当把各层的brup循环赋值给brup_tm2D后，我们跟着exch_hx，在下一行仿照exch_hx进行引用和声明，而它们同在一个叫subroutine ysu2d的子程序中：
   形参新增：
      subroutine ysu2d(j,ux,vx,tx,qx,p2d,p2di,pi2d,                            &
                     ...
                              exch_hx,                                                    &
                                brup_tm2D,                                                 &       !TM added

   声明新增：
         intent(inout) ::                                     exch_hx
         real, dimension( ims:ime, kms:kme )                                  , &
         intent(inout) ::                                     brup_tm2D    !TM added

   （2）修改完成后，我们发现subroutine ysu这个程序调用了subroutine ysu2d这个子程序，按照exch_hx我们捋到call ysu2d中有这么一句二维转三维的赋值语句，于是我们跟着它对brup_tm2D这个变量也进行了赋值：
         ,exch_hx=exch_h(ims,kms,j)                                     &
         ,brup_tm2D=brup_tm(ims,kms,j)                               & !TM added
      然后继续照猫画虎，在subroutine ysu中进行引用和声明：
      形参新增：
         subroutine ysu(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d,             &
         .....
         exch_h,                                                    &
         brup_tm,                                                    & ! TM added

      声明新增：
         real,    dimension( ims:ime, kms:kme, jms:jme )                         , &
         intent(inout) ::                                     exch_h
         real,    dimension( ims:ime, kms:kme, jms:jme )                         , &
         intent(inout) ::                                     brup_tm ! TM added

      由于subroutine ysu就是这个模块的主程序，因此我们发现，到一步，我们已经完成了几乎50%的工作。

   （3）向上我们就需要修改第二层驱动程序，由于我们修改的是边界层参数化方案，所以我们需要编写module_pbl_driver.F，其他参数化方案就打开对应驱动就可以了，那我们还是跟着exch_h走，在 CASE (YSUSCHEME)中，我们补充这个程序在调用了YSU参数化方案时：
            CALL ysu(                                           &
            U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy    &
            ......
            ,EXCH_H=exch_h,REGIME=regime                         &
         ,BRUP_TM=brup_tm                                  & !TM added
照葫芦画瓢，在SUBROUTINE pbl_driver程序中新增：
形参：
SUBROUTINE pbl_driver(                                        &
                  itimestep,dt,u_frame,v_frame                   &
                  ......
                  ,exch_h,exch_m,akhs,akms                         &
                 ,brup_tm                                        & !TM added

声明：
REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),          &
            INTENT(INOUT) ::                      RUBLTEN, &
                                                               RVBLTEN, &
                                                               RTHBLTEN, &
                                       EXCH_H,EXCH_M,TKE_PBL, &
                                                   BRUP_TM, & ! TM added

      这一步修改就结束了，是不是超简单。其实在跟随exch_h的脚步中，我们发现我们也是在跟随不同的子程序：先是subroutine ysu2d，然后我们发现subroutine ysu调用了它，紧接着SUBROUTINE pbl_driver调用了subroutine ysu，随着程序一步步向上追寻，我们也即将到达最后的目标，最上层的求解器。

      （4）在dyn_em文件夹中grep -wr pbl_driver，我们发现是哪个小可爱在更上层调用了这个驱动呢？module_first_rk_step_part1.F
打开这个程序，我们立刻又遇到了老熟人，exch_h，于是我们继续跟着它写：
               CALL pbl_driver(                                                             &
   &       AKHS=grid%akhs       ,AKMS=grid%akms                            &
               ......
   &       ,EXCH_H=grid%exch_h    ,EXCH_M=grid%exch_m                         &
    &       ,BRUP_TM=grid%brup_tm                                              & !TM added
至此，所有的修改工作已经完成，我们只需要在注册表中加入一句：state real BRUP_TM    ikj misc       1       Z    rh       "BRUP_TM"          "bulk Richardson number" " ,然后重新 ./compile em_real >& compile.log，即可。

   肝论文总有一种掏空自己的感觉，而没有功利性地阅读这些基础的程序、基础的文献，就是又一次的快乐充电。
   忘不了那个炎热的夏天，学长发我《WRF安装运行入门指南》的pdf，那种从安装到运行的简单快乐。有时候走的太远了，可以回来看看初心，阳后这一周的病假，放着文章不改，看着熟悉的代码，真正得到了一次身心的放松spa。

   祝大家科研生活愉快。

tbag · 发表于 2023-11-14 15:47:35

好东西，支持一下。
其实WRF自定义的东西很多很多，我之前问过一个气象公司，他们跑的都是超级自定义的自己的版本。

tianmeng · 发表于 2023-2-12 20:48:34

本帖最后由 tianmeng 于 2023-2-12 22:13 编辑

WRF内部程序学习第二弹写什么还没想好，有感兴趣的点吗？没人指导，纯自学，所以有问题还请大家随时指正。其实我想写一写WRF内部程序的荒川C网格设计，因为我在本科的时候就有一个很大的疑问，动量场变量总是比质量场变量多一层，那么大家画垂直廓线图的时候怎么对应高度呢？

nightfalls · 发表于 2023-2-13 17:00:36

太厉害了。我只会print出来在log文件里看

。

tianmeng · 发表于 2023-2-14 14:31:31

nightfalls 发表于 2023-2-13 17:00
太厉害了。我只会print出来在log文件里看。

其实找到计算程序之后，有些中间l量也可以通过复制公式，利用直接输出变量进行计算，和这样一层一层输出的结果是一致的，这样就节约了编译时间

谷梁陌尘 · 发表于 2023-2-17 09:15:02

厉害，收藏学习

Houckin · 发表于 2023-2-17 11:17:37

有兴趣，让我先把表层的先搞懂先

zhenbuxiangshuo · 发表于 2023-2-17 16:22:37

谢谢分享，太厉害了

BOBO.jas · 发表于 2023-3-7 20:07:49

nightfalls 发表于 2023-2-13 17:00
太厉害了。我只会print出来在log文件里看。

你好哇，这个log文件是指compile以后生成的compile.log文件吗，还是其他的啥文件哇

nightfalls · 发表于 2023-3-8 16:03:14

BOBO.jas 发表于 2023-3-7 20:07
你好哇，这个log文件是指compile以后生成的compile.log文件吗，还是其他的啥文件哇

不是。是运行中产生的rsl.out。

BOBO.jas · 发表于 2023-3-8 20:34:10

nightfalls 发表于 2023-3-8 16:03
不是。是运行中产生的rsl.out。

好滴，谢谢

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

WRF内部程序学习，第一弹：中间变量的输出

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

点评