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

tianmeng

      随着研究深入，我们常常需要直接上手修改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

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

tianmeng

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

nightfalls

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

tianmeng

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

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

谷梁陌尘

厉害，收藏学习

Houckin

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

zhenbuxiangshuo

谢谢分享，太厉害了

BOBO.jas

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

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

nightfalls

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

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

BOBO.jas

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

好滴，谢谢