本内容相关视频讲解：WRF物理参数化方案的修改
背景——模型耦合
什么是模型耦合(Couple)？
经过搜索，可以看到如下定义：
Coupling refers to communication and interchange of information between models.
直译而言，模式的耦合可以理解为模式间的信息交换与相互作用，那么，又怎样理解模式间的信息交换与作用？模式间的信息又指什么？
按笔者个人的经验，模式的“信息”可以理解为：模式的输入与输出，而模式间的耦合，可以理解为两个及以上的模型之间输入与输出变量间的相互连接。
比如说，WRF-Chem便是一个典型的在线耦合模式，中尺度天气预测模式WRF提供了气象边界场与气象条件，这些输出的气象条件作为一个输入，用来计算大气内各种化学物质的传输、沉降过程，同时，化学模型计算所输出的各类化学物质的源汇，又可以作为一个输入，去计算其当前气象条件下，气溶胶粒子对气象条件的影响，这种模式间不同输入与输出的连接与迭代，便是”耦合“。
而最简单的理解就是：一个模式输出的结果，作为另一个模式的输入，便可以称之为”耦合“。
这类定义随着学科的差异而有着不同的定义，在地学领域内，模型的耦合可分为在线（Online） 与**离线（Offline）**耦合。
离线耦合通常仅仅只是将一个模型的输出结果提供给另一个模型，而并不考虑二者间的反馈作用，比方说WRF/SLEUTH的耦合，将SLEUTH输出的土地利用类型变化输入至WRF中，以此来研究土地利用类型对气象天气的影响，不考虑气象条件对土地利用类型的影响。离线耦合的模式彼此间较为独立，运行也是分开运行。
在线耦合则是考虑了模式间的反馈作用，WRF-Chem便是在线耦合模式的代表，在线耦合模式模块间存在着内部间的相互迭代，考虑两个模块间的相互作用，在线耦合也被称作”fully coupled“
那么，当我们想实现一个在线的模式耦合，我们应当如何做？我们如何应用一个新的物理方案到WRF中呢？
关于增加一个新的物理参数化方案到WRF中的步骤，已有博客Implementing a new physics scheme进行了极好的概括，不过，由于年代较远且细节不足等原因，仅仅依赖此篇博客，或许仍然会存在一些问题，在此，我将以自己所做的工作为实例，讲述耦合新的方案到WRF中的流程。
技术路线——WRF结构
在正式开始修改源码前，我们首先应当对WRF各个模块与其间的关系有一个初步的了解。WRF模式，尤其是模拟真实大气的部分结构较为复杂，会考虑真实三维大气从地面到高空的全部过程。
总体而言，WRF代码的框架可以概括为以下部分：
- 用来读取WPS前处理文件、完成初始化与求解的dyn_em
- 对于网格与模拟域的处理：frame
- 读取变量、变量处理、插值：share
- 物理参数化方案: phys
由于WRF由多个模块组成，模块与模块间存在相互多次调用的情况，在修改WRF源码时，应当注意修改部分的对应。
物理参数方案
WRF中的物理参数化方案包括长短波辐射方案、积云方案、云微物理方案、边界层参数化方案、陆面方案等。
这些方案的代码全都在WRF/phys的文件夹下，通常，除了各个方案的代码以外，还包括物理常数模块、以及驱动物理参数化的模块，如，地面方案的surface_driver模块。
对物理源码进行修改时，除了变量的增添删改，还应当模块的使用与对应模块的修改，如，当你给CLM模块中的函数增加了一个变量，除了在函数内与模块内声明外，还应当在使用了CLM木块的surface_driver增加变量的声明。
初始化
输入WRF中的任何变量，都需要在module_initialize_real.F中进行初始化，即，在wrfinput文件中生成，以及 first_rk_step模块根据选择的参数化方案分配处理变量。
输入变量
在模式初始化时，如果需要增加新的变量并进行初始化的话，需要在相应的模块optional_input
进行增加，否则不会对新增变量进行处理。
变量描述——Registry
WRF会根据你的namelist设置以及Registry文件中对于变量以及pacakge的描述与定义，自动生成一段源代码，该代码会控制wrfout结果的输入与输出。
在增加变量描述时，一定要确保变量的描述特征与所需求一致，比如维度以及变量的驶入输出以及插值覆盖，这些属性都会影响运行时WRF对新增变量的处理。
关于Registry的具体运算，建议参考：How to use WRF Registry
并行运算
WRF在运行时，有时需要读取外部数据，比较典型的代表便是运行陆面方案。
由于WRF运行往往需要进行并行运算，因此，在读取外部数据后，还需要将其传递给其他进程以便并行运算。
通常，WRF中读取外部数据的会使用命令
IF ( wrf_dm_on_monitor() ) THEN
监测WRF是否在进程中，如果WRF在运行，则Open Read Close读取数据。
随后，使用，CALL wrf_dm_bcast_real 函数命令，将读取的数据传播至其他并行进程。
可参考：How to parallelize a part of reading an aerosol file in module_initialize_real.F in dyn_em directory
模块编译
WRF作为一个中尺度模式，其各个过程都存储与相应独立的模块（module）中。每个moudle相互独立，具有不同的功能，又彼此联系耦合，当一个module想使用另一个module时，则需要在module开头 USE moudle module name，比如说，在CLM模块中，经常使用：
use shr_kind_mod, only: r8 => shr_kind_r8
use module_cam_support, only: endrun
使用shr_kind模块中的r8提高变量的精度（双浮点），使用cam_support模块用来DEBUG，并停止模式运行。
流程总结
总结一下WRF修改物理参数化方案的流程。
针对只修改原方案，不增加方案的：
- 修改对应phys/文件夹下的代码，包括计算模块与驱动模块
- 如果不涉及变量的改变，则直接编译即可，如果涉及变量的增添，在修改初始化部分代码与变量输入模块代码。
- 修改Registry文件，增加变量描述、
针对需要增加新方案的，参考博客Implementing a new physics scheme，相比起简单修改方案，需要多出以下步骤：
- 修改Makefile文件，编译自己的新方案模块。
- 分配模块空间，确保在namelist.input中可以选择。
- 修改驱动文件，创建一个新Case，同时在Regeistry中增加pacakge。
如果涉及到变量的增加，则需同样的步骤。
上述修改完并编译成功后，可能会出现变量未正确输入或运算的问题，建议在代码中加上write函数输出确认变量的输入输出过程。
此外，可能存在即使模式正确耦合运行，却运行速度降低的问题，一般是由于并行部分处理有问题，需要检查并行部分的处理。
实例
基本信息
对于不同的研究方向与问题而言，对于WRF的修改必然都会有所不同，但是，有些东西是共通的，不同的只有其中的物理过程差异，而WRF架构的变量输入、输出、并行、模块这一部分则是每一个对于WRF 的modification需要了解的。
我本次修改想实现的目的是：
- 将CLM 陆面模式中的SNICAR积雪反照率方案耦合至Noah-mp模式中。
- 增加一个输入变量，MBC，表示雪中BC浓度，此变量可通过MET_EM文件输入至wrfinput中。
- 增加的变量MBC将参与新耦合的雪反照率方案中。
根据我的需求，我对原本的SNICAR方案进行了简化：
- 不考虑雪变质、老化等过程，忽略雪粒有效半径的计算，将雪有效半径固定。
- 不考虑BC 的粒径与coated作用，仅考虑浓度。
- 将晴天与非晴天单独分开处理，即将散射与直射分开，只考虑晴天于非晴天两种情况。
下面开始讲我的实际操作：
函数添加
这一部分主要在于不同变量间的声明，注意，耦合的关键在于新函数计算所需的变量可以由Noah-MP提供，在不同函数间，声明的变量具有不同的名称与相同的物理含义，需要辨别。
创建一个新的函数，反照率函数：
SUBROUTINE SNOWALB_SNICAR_DIR (parameters,NBAND,NSNOW,ISNOW,COSZ,BCS,UALB,SNICE,SNLIQ,SNEQV,ALBSND,ALBSNI)
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! INput
type (noahmp_parameters), INTENT(IN) :: parameters
INTEGER,INTENT(IN) :: NBAND !number of waveband classes=numrad
INTEGER,INTENT(IN) :: ISNOW !actual no. of snow layers
INTEGER,INTENT(IN) :: NSNOW !maximum no. of snow layers
REAL,INTENT(IN) :: COSZ !cosINe solar zenith angle
REAL(r8),INTENT(IN) :: BCS !bc IN snow (ng/g)
REAL,INTENT(IN) :: UALB !Albedo of underlyINg surface
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNICE !snow ice mass (kg/m2)
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNLIQ !snow liq mass (kg/m2)
REAL, INTENT(IN) :: SNEQV !snow water eqv. [mm]
! REAL(r8), INTENT(IN) :: h2osno_liq(-NSNOW+1:0) ! liquid water content (col,lyr) [kg/m2]
! REAL(r8), INTENT(IN) :: h2osno_ice(-NSNOW+1:0) ! ice content (col,lyr) [kg/m2]
! INTEGER, INTENT(IN) :: snw_rds(-NSNOW+1:0) ! snow effective radius (col,lyr) [microns, m^-6]
! REAL(r8), INTENT(IN) :: mss_cnc_aer_in(-NSNOW+1:0) ! mass concentration of all aerosol species (col,lyr,aer) [kg/kg],此处只考虑一种
! REAL(r8), INTENT(IN) :: albsfc(numrad) ! albedo of surface underlyINg snow (col,bnd) [frc]
! OUTput
REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSND !snow albedo for direct(1=vis, 2=nir)
REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSNI !snow albedo for diffuse
! ---------------------------------------------------------------------------------------------
! ------------------------ local variables ----------------------------------------------------
!-----------------------------------------------------------------------
!
! variables for snow radiative transfer calculations
! Local variables representINg sINgle-column values of arrays:
INTEGER :: snl_lcl ! negative number of snow layers [nbr]
INTEGER :: snw_rds_lcl(-NSNOW+1:0) ! snow effective radius [m^-6]
REAL(r8):: flx_slrd_lcl(1:numrad_snw) ! direct beam INcident irradiance [W/m2] (set to 1)
REAL(r8):: flx_slri_lcl(1:numrad_snw) ! diffuse INcident irradiance [W/m2] (set to 1)
REAL(r8):: mss_cnc_aer_lcl(-NSNOW+1:0) ! aerosol mass concentration (lyr,aer_nbr) [kg/kg]
REAL(r8):: h2osno_lcl ! total column snow mass [kg/m2]
REAL(r8):: h2osno_liq_lcl(-NSNOW+1:0) ! liquid water mass [kg/m2]
REAL(r8):: h2osno_ice_lcl(-NSNOW+1:0) ! ice mass [kg/m2]
REAL(r8):: albsfc_lcl(1:numrad_snw) ! albedo of underlyINg surface [frc]
REAL(r8):: ss_alb_snw_lcl(-NSNOW+1:0) ! sINgle-scatter albedo of ice graINs (lyr) [frc]
REAL(r8):: asm_prm_snw_lcl(-NSNOW+1:0) ! asymmetry parameter of ice graINs (lyr) [frc]
REAL(r8):: ext_cff_mss_snw_lcl(-NSNOW+1:0) ! mass extINction coefficient of ice graINs (lyr) [m2/kg]
REAL(r8):: ss_alb_aer_lcl ! sINgle-scatter albedo of aerosol species (aer_nbr) [frc]
REAL(r8):: asm_prm_aer_lcl ! asymmetry parameter of aerosol species (aer_nbr) [frc]
REAL(r8):: ext_cff_mss_aer_lcl ! mass extINction coefficient of aerosol species (aer_nbr) [m2/kg]
REAL(r8) :: albout(1:NBAND) ! snow albedo, averaged INto 2 bands (=0 if no sun or no snow) (col,bnd) [frc]
REAL(r8) :: flx_abs(-NSNOW+1:1,1:NBAND) ! absorbed flux IN each layer per unit flux INcident on top of snowpack (col,lyr,bnd) [frc]
! Other local variables
INTEGER :: APRX_TYP ! two-stream approximation type (1=EddINgton, 2=Quadrature, 3=Hemispheric Mean) [nbr]
INTEGER :: DELTA ! flag to use Delta approximation (Joseph, 1976) (1= use, 0= don't use)
REAL(r8):: flx_wgt(1:numrad_snw) ! weights applied to spectral bands, specific to direct and diffuse cases (bnd) [frc]
INTEGER :: flg_nosnl ! flag: =1 if there is snow, but zero snow layers, =0 if at least 1 snow layer [flg]
INTEGER :: trip ! flag: =1 to redo RT calculation if result is unREAListic
INTEGER :: flg_dover ! defINes conditions for RT redo (explaINed below)
REAL(r8):: albedo_1 ! temporary snow albedo [frc]
REAL(r8):: flx_sum ! temporary summation variable for NIR weightINg
REAL(r8):: albout_lcl(1:numrad_snw) ! snow albedo by band [frc]
REAL(r8):: flx_abs_lcl(-NSNOW+1:1,1:numrad_snw)! absorbed flux per unit INcident flux at top of snowpack (lyr,bnd) [frc]
REAL(r8):: L_snw(-NSNOW+1:0) ! h2o mass (liquid+solid) IN snow layer (lyr) [kg/m2]
REAL(r8):: tau_snw(-NSNOW+1:0) ! snow optical depth (lyr) [unitless]
REAL(r8):: L_aer(-NSNOW+1:0) ! aerosol mass IN snow layer (lyr,nbr_aer) [kg/m2]
REAL(r8):: tau_aer(-NSNOW+1:0) ! aerosol optical depth (lyr,nbr_aer) [unitless]
REAL(r8):: tau_sum ! cumulative (snow+aerosol) optical depth [unitless]
REAL(r8):: tau_clm(-NSNOW+1:0) ! column optical depth from layer bottom to snowpack top (lyr) [unitless]
REAL(r8):: omega_sum ! temporary summation of sINgle-scatter albedo of all aerosols [frc]
REAL(r8):: g_sum ! temporary summation of asymmetry parameter of all aerosols [frc]
REAL(r8):: tau(-NSNOW+1:0) ! weighted optical depth of snow+aerosol layer (lyr) [unitless]
REAL(r8):: omega(-NSNOW+1:0) ! weighted sINgle-scatter albedo of snow+aerosol layer (lyr) [frc]
REAL(r8):: g(-NSNOW+1:0) ! weighted asymmetry parameter of snow+aerosol layer (lyr) [frc]
REAL(r8):: tau_star(-NSNOW+1:0) ! transformed (i.e. Delta-EddINgton) optical depth of snow+aerosol layer (lyr) [unitless]
REAL(r8):: omega_star(-NSNOW+1:0) ! transformed (i.e. Delta-EddINgton) SSA of snow+aerosol layer (lyr) [frc]
REAL(r8):: g_star(-NSNOW+1:0) ! transformed (i.e. Delta-EddINgton) asymmetry paramater of snow+aerosol layer (lyr) [frc]
INTEGER :: g_idx, c_idx, l_idx ! gridcell, column, and landunit INdices [idx]
INTEGER :: bnd_idx ! spectral band INdex (1 <= bnd_idx <= numrad_snw) [idx]
INTEGER :: rds_idx ! snow effective radius INdex for retrievINg Mie parameters from lookup table [idx]
INTEGER :: snl_btm ! INdex of bottom snow layer (0) [idx]
INTEGER :: snl_top ! INdex of top snow layer (-4 to 0) [idx]
INTEGER :: fc ! column filter INdex
INTEGER :: i ! layer INdex [idx]
INTEGER :: j ! aerosol number INdex [idx]
INTEGER :: n ! tridiagonal matrix INdex [idx]
INTEGER :: m ! secondary layer INdex [idx]
INTEGER :: ix,k ! an INdex
REAL(r8):: F_direct(-NSNOW+1:0) ! direct-beam radiation at bottom of layer INterface (lyr) [W/m^2]
REAL(r8):: F_net(-NSNOW+1:0) ! net radiative flux at bottom of layer INterface (lyr) [W/m^2]
REAL(r8):: F_abs(-NSNOW+1:0) ! net absorbed radiative energy (lyr) [W/m^2]
REAL(r8):: F_abs_sum ! total absorbed energy IN column [W/m^2]
REAL(r8):: F_sfc_pls ! upward radiative flux at snowpack top [W/m^2]
REAL(r8):: F_btm_net ! net flux at bottom of snowpack [W/m^2]
REAL(r8):: F_sfc_net ! net flux at top of snowpack [W/m^2]
REAL(r8):: energy_sum ! sum of all energy terms; should be 0.0 [W/m^2]
REAL(r8):: F_direct_btm ! direct-beam radiation at bottom of snowpack [W/m^2]
REAL(r8):: mu_not ! cosINe of solar zenith angle (used locally) [frc]
INTEGER :: err_idx ! counter for number of times through error loop [nbr]
REAL(r8):: lat_coord ! gridcell latitude (debuggINg only)
REAL(r8):: lon_coord ! gridcell longitude (debuggINg only)
INTEGER :: sfctype ! underlyINg surface type (debuggINg only)
REAL(r8):: pi ! 3.1415...
LOGICAL :: opened
INTEGER :: lu_unit=10
logical, external :: wrf_dm_on_monitor
! INtermediate variables for radiative transfer approximation:
REAL(r8):: gamma1(-NSNOW+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: gamma2(-NSNOW+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: gamma3(-NSNOW+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: gamma4(-NSNOW+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: lambda(-NSNOW+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: GAMMA(-NSNOW+1:0) ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: mu_one ! two-stream coefficient from Toon et al. (lyr) [unitless]
REAL(r8):: e1(-NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (lyr)
REAL(r8):: e2(-NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (lyr)
REAL(r8):: e3(-NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (lyr)
REAL(r8):: e4(-NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (lyr)
REAL(r8):: C_pls_btm(-NSNOW+1:0) ! INtermediate variable: upward flux at bottom INterface (lyr) [W/m2]
REAL(r8):: C_mns_btm(-NSNOW+1:0) ! INtermediate variable: downward flux at bottom INterface (lyr) [W/m2]
REAL(r8):: C_pls_top(-NSNOW+1:0) ! INtermediate variable: upward flux at top INterface (lyr) [W/m2]
REAL(r8):: C_mns_top(-NSNOW+1:0) ! INtermediate variable: downward flux at top INterface (lyr) [W/m2]
REAL(r8):: A(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: B(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: D(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: E(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: AS(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: DS(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: X(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
REAL(r8):: Y(-2*NSNOW+1:0) ! tri-diag INtermediate variable from Toon et al. (2*lyr)
! ---------------------------------------------------------------------------------------------
! zero albedos for all poINts
ALBSND(1: NBAND) = 0.
ALBSNI(1: NBAND) = 0.
! DefINe constants
pi = SHR_CONST_PI
! always use Delta approximation for snow
DELTA = 1
! Zero absorbed radiative fluxes:
do i=-NSNOW+1,1,1
flx_abs_lcl(:,:) = 0._r8
flx_abs(i,:) = 0._r8
enddo
h2osno_lcl = SNEQV
if (h2osno_lcl > min_snw) then
! Assign local (single-column) variables to global values
! If there is snow, but zero snow layers, we must create a layer locally.
! This layer is presumed to have the fresh snow effective radius.
if ( ISNOW > -1) then
flg_nosnl = 1
snl_lcl = -1
h2osno_ice_lcl(0) = h2osno_lcl
h2osno_liq_lcl(0) = 0._r8
snw_rds_lcl(0) = nint(snw_rds_min)
else
flg_nosnl = 0
snl_lcl = ISNOW
h2osno_liq_lcl(:) = SNLIQ(:)
h2osno_ice_lcl(:) = SNICE(:)
snw_rds_lcl(:) = 150
endif
! Set local aerosol array
snl_btm = 0
snl_top = snl_lcl+1
! Set local aerosol array
mss_cnc_aer_lcl(:) = BCS*1.0_r8
! Set spectral underlying surface albedos to their corresponding VIS or NIR albedos
albsfc_lcl(1) = UALB
albsfc_lcl(nir_bnd_bgn:nir_bnd_end) = UALB
! Incident flux weighting parameters
! - sum of all VIS bands must equal 1
! - sum of all NIR bands must equal 1
!
! Spectral bands (5-band case)
! Band 1: 0.3-0.7um (VIS)
! Band 2: 0.7-1.0um (NIR)
! Band 3: 1.0-1.2um (NIR)
! Band 4: 1.2-1.5um (NIR)
! Band 5: 1.5-5.0um (NIR)
! 5-band weights DIRECT
flx_wgt(1) = 1._r8
flx_wgt(2) = 0.49352158521175_r8
flx_wgt(3) = 0.18099494230665_r8
flx_wgt(4) = 0.12094898498813_r8
flx_wgt(5) = 0.20453448749347_r8
! Loop over snow spectral bands
do bnd_idx = 1,numrad_snw
mu_not = COSZ ! must set here, because of error handling
flg_dover = 1 ! default is to redo
err_idx = 0 ! number of times through loop
do while (flg_dover > 0)
! DEFAULT APPROXIMATIONS:
! VIS: Delta-Eddington
! NIR (all): Delta-Hemispheric Mean
! WARNING: DO NOT USE DELTA-EDDINGTON FOR NIR DIFFUSE - this sometimes results in negative albedo
!
! ERROR CONDITIONS:
! Conditions which cause "trip", resulting in redo of RT approximation:
! 1. negative absorbed flux
! 2. total absorbed flux greater than incident flux
! 3. negative albedo
! NOTE: These errors have only been encountered in spectral bands 4 and 5
!
! ERROR HANDLING
! 1st error (flg_dover=2): switch approximation (Edd->HM or HM->Edd)
! 2nd error (flg_dover=3): change zenith angle by 0.02 (this happens about 1 in 10^6 cases)
! 3rd error (flg_dover=4): switch approximation with new zenith
! Subsequent errors: repeatedly change zenith and approximations...
if ( bnd_idx == 1) then
if (flg_dover == 2) then
APRX_TYP = 3
elseif (flg_dover == 3) then
APRX_TYP = 1
if (COSZ > 0.5_r8) then
mu_not = mu_not - 0.02_r8
else
mu_not = mu_not + 0.02_r8
endif
elseif (flg_dover == 4) then
APRX_TYP = 3
else
APRX_TYP = 1
endif
else
if (flg_dover == 2) then
APRX_TYP = 1
elseif (flg_dover == 3) then
APRX_TYP = 3
if (COSZ > 0.5_r8) then
mu_not = mu_not - 0.02_r8
else
mu_not = mu_not + 0.02_r8
endif
elseif (flg_dover == 4) then
APRX_TYP = 1
else
APRX_TYP = 3
endif
endif
! Set direct or diffuse incident irradiance to 1
! (This has to be within the bnd loop because mu_not is adjusted in rare cases)
flx_slrd_lcl(bnd_idx) = 1._r8/(mu_not*pi) ! this corresponds to incident irradiance of 1.0
flx_slri_lcl(bnd_idx) = 0._r8
! Pre-emptive error handling: aerosols can reap havoc on these absorptive bands.
! Since extremely high soot concentrations have a negligible effect on these bands, zero them.
if ( (numrad_snw == 5).and.((bnd_idx == 5).or.(bnd_idx == 4)) ) then
mss_cnc_aer_lcl = 0._r8
endif
! Define local Mie parameters based on snow grain size and aerosol species,
! retrieved from a lookup table.
do i=snl_top,snl_btm,1
rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
! snow optical properties (direct radiation)
ss_alb_snw_lcl(i) = ss_alb_snw_drc(rds_idx,bnd_idx)
asm_prm_snw_lcl(i) = asm_prm_snw_drc(rds_idx,bnd_idx)
ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_drc(rds_idx,bnd_idx)
enddo
! aerosol species optical properties only BC nocoated
ss_alb_aer_lcl = ss_alb_bc1(1,bnd_idx)
asm_prm_aer_lcl = asm_prm_bc1(1,bnd_idx)
ext_cff_mss_aer_lcl = ext_cff_mss_bc1(1,bnd_idx)
! 1. snow and aerosol layer column mass (L_snw, L_aer [kg/m^2])
! 2. optical Depths (tau_snw, tau_aer)
! 3. weighted Mie properties (tau, omega, g)
! Weighted Mie parameters of each layer
do i=snl_top,snl_btm,1
L_snw(i) = h2osno_ice_lcl(i)+h2osno_liq_lcl(i)
tau_snw(i) = L_snw(i)*ext_cff_mss_snw_lcl(i)
L_aer(i) = L_snw(i)*mss_cnc_aer_lcl(i)
tau_aer(i) = L_aer(i)*ext_cff_mss_aer_lcl
tau_sum = 0._r8
omega_sum = 0._r8
g_sum = 0._r8
tau_sum = tau_sum + tau_aer(i)
omega_sum = omega_sum + (tau_aer(i)*ss_alb_aer_lcl)
g_sum = g_sum + (tau_aer(i)*ss_alb_aer_lcl*asm_prm_aer_lcl)
tau(i) = tau_sum + tau_snw(i)
if(tau(i) == 0) then
write(6,*) 'FATAL ERROR in SNICAR RT, tau(',i,') is the denominatoer can not equal to ',tau(i)
call endrun()
end if
omega(i) = (1/tau(i))*(omega_sum+(ss_alb_snw_lcl(i)*tau_snw(i)))
g(i) = (1/(tau(i)*omega(i)))*(g_sum+ (asm_prm_snw_lcl(i)*ss_alb_snw_lcl(i)*tau_snw(i)))
enddo
! DELTA transformations, if requested
if (DELTA == 1) then
do i=snl_top,snl_btm,1
g_star(i) = g(i)/(1+g(i))
omega_star(i) = ((1-(g(i)**2))*omega(i)) / (1-(omega(i)*(g(i)**2)))
tau_star(i) = (1-(omega(i)*(g(i)**2)))*tau(i)
enddo
else
do i=snl_top,snl_btm,1
g_star(i) = g(i)
omega_star(i) = omega(i)
tau_star(i) = tau(i)
enddo
endif
! Total column optical depth:
! tau_clm(i) = total optical depth above the bottom of layer i
tau_clm(snl_top) = 0._r8
do i=snl_top+1,snl_btm,1
tau_clm(i) = tau_clm(i-1)+tau_star(i-1)
enddo
! Direct radiation at bottom of snowpack:
F_direct_btm = albsfc_lcl(bnd_idx)*mu_not*exp(-(tau_clm(snl_btm)+tau_star(snl_btm))/mu_not)*pi*flx_slrd_lcl(bnd_idx)
! Intermediates
! Gamma values are approximation-specific.
! Eddington
if (APRX_TYP==1) then
do i=snl_top,snl_btm,1
gamma1(i) = (7-(omega_star(i)*(4+(3*g_star(i)))))/4
gamma2(i) = -(1-(omega_star(i)*(4-(3*g_star(i)))))/4
gamma3(i) = (2-(3*g_star(i)*mu_not))/4
gamma4(i) = 1-gamma3(i)
mu_one = 0.5
enddo
! Quadrature
elseif (APRX_TYP==2) then
do i=snl_top,snl_btm,1
gamma1(i) = (3**0.5)*(2-(omega_star(i)*(1+g_star(i))))/2
gamma2(i) = omega_star(i)*(3**0.5)*(1-g_star(i))/2
gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2
gamma4(i) = 1-gamma3(i)
mu_one = 1/(3**0.5)
enddo
! Hemispheric Mean
elseif (APRX_TYP==3) then
do i=snl_top,snl_btm,1
gamma1(i) = 2 - (omega_star(i)*(1+g_star(i)))
gamma2(i) = omega_star(i)*(1-g_star(i))
gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2
gamma4(i) = 1-gamma3(i)
mu_one = 0.5
enddo
endif
! Intermediates for tri-diagonal solution
do i=snl_top,snl_btm,1
lambda(i) = sqrt(abs((gamma1(i)**2) - (gamma2(i)**2)))
GAMMA(i) = gamma2(i)/(gamma1(i)+lambda(i))
e1(i) = 1+(GAMMA(i)*exp(-lambda(i)*tau_star(i)))
e2(i) = 1-(GAMMA(i)*exp(-lambda(i)*tau_star(i)))
e3(i) = GAMMA(i) + exp(-lambda(i)*tau_star(i))
e4(i) = GAMMA(i) - exp(-lambda(i)*tau_star(i))
enddo !enddo over snow layers
! Intermediates for tri-diagonal solution
do i=snl_top,snl_btm,1
C_pls_btm(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-(tau_clm(i)+tau_star(i))/mu_not)* &
(((gamma1(i)-(1/mu_not))*gamma3(i))+ &
(gamma4(i)*gamma2(i))))/((lambda(i)**2)-(1/(mu_not**2)))
C_mns_btm(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-(tau_clm(i)+tau_star(i))/mu_not)* &
(((gamma1(i)+(1/mu_not))*gamma4(i))+ &
(gamma2(i)*gamma3(i))))/((lambda(i)**2)-(1/(mu_not**2)))
C_pls_top(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-tau_clm(i)/mu_not)*(((gamma1(i)-(1/mu_not))* &
gamma3(i))+(gamma4(i)*gamma2(i))))/((lambda(i)**2)-(1/(mu_not**2)))
C_mns_top(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
exp(-tau_clm(i)/mu_not)*(((gamma1(i)+(1/mu_not))* &
gamma4(i))+(gamma2(i)*gamma3(i))))/((lambda(i)**2)-(1/(mu_not**2)))
enddo
! Coefficients for tridiaganol matrix solution
do i=2*snl_lcl+1,0,1
!Boundary values for i=1 and i=2*snl_lcl, specifics for i=odd and i=even
if (i==(2*snl_lcl+1)) then
A(i) = 0
B(i) = e1(snl_top)
D(i) = -e2(snl_top)
E(i) = flx_slri_lcl(bnd_idx)-C_mns_top(snl_top)
elseif(i==0) then
A(i) = e1(snl_btm)-(albsfc_lcl(bnd_idx)*e3(snl_btm))
B(i) = e2(snl_btm)-(albsfc_lcl(bnd_idx)*e4(snl_btm))
D(i) = 0
E(i) = F_direct_btm-C_pls_btm(snl_btm)+(albsfc_lcl(bnd_idx)*C_mns_btm(snl_btm))
elseif(mod(i,2)==-1) then ! If odd and i>=3 (n=1 for i=3)
n=floor(i/2.0)
A(i) = (e2(n)*e3(n))-(e4(n)*e1(n))
B(i) = (e1(n)*e1(n+1))-(e3(n)*e3(n+1))
D(i) = (e3(n)*e4(n+1))-(e1(n)*e2(n+1))
E(i) = (e3(n)*(C_pls_top(n+1)-C_pls_btm(n)))+(e1(n)*(C_mns_btm(n)-C_mns_top(n+1)))
elseif(mod(i,2)==0) then ! If even and i<=2*snl_lcl
n=(i/2)
A(i) = (e2(n+1)*e1(n))-(e3(n)*e4(n+1))
B(i) = (e2(n)*e2(n+1))-(e4(n)*e4(n+1))
D(i) = (e1(n+1)*e4(n+1))-(e2(n+1)*e3(n+1))
E(i) = (e2(n+1)*(C_pls_top(n+1)-C_pls_btm(n)))+(e4(n+1)*(C_mns_top(n+1)-C_mns_btm(n)))
endif
enddo
AS(0) = A(0)/B(0)
DS(0) = E(0)/B(0)
do i=-1,(2*snl_lcl+1),-1
X(i) = 1/(B(i)-(D(i)*AS(i+1)))
AS(i) = A(i)*X(i)
DS(i) = (E(i)-(D(i)*DS(i+1)))*X(i)
enddo
Y(2*snl_lcl+1) = DS(2*snl_lcl+1)
do i=(2*snl_lcl+2),0,1
Y(i) = DS(i)-(AS(i)*Y(i-1))
enddo
! Downward direct-beam and net flux (F_net) at the base of each layer:
do i=snl_top,snl_btm,1
F_direct(i) = mu_not*pi*flx_slrd_lcl(bnd_idx)*exp(-(tau_clm(i)+tau_star(i))/mu_not)
F_net(i) = (Y(2*i-1)*(e1(i)-e3(i))) + (Y(2*i)*(e2(i)-e4(i))) + &
C_pls_btm(i) - C_mns_btm(i) - F_direct(i)
enddo
! Upward flux at snowpack top:
F_sfc_pls = (Y(2*snl_lcl+1)*(exp(-lambda(snl_top)*tau_star(snl_top))+ &
GAMMA(snl_top))) + (Y(2*snl_lcl+2)*(exp(-lambda(snl_top)* &
tau_star(snl_top))-GAMMA(snl_top))) + C_pls_top(snl_top)
! Net flux at bottom = absorbed radiation by underlying surface:
F_btm_net = -F_net(snl_btm)
! Bulk column albedo and surface net flux
albedo_1 = F_sfc_pls/((mu_not*pi*flx_slrd_lcl(bnd_idx))+flx_slri_lcl(bnd_idx))
F_sfc_net = F_sfc_pls - ((mu_not*pi*flx_slrd_lcl(bnd_idx))+flx_slri_lcl(bnd_idx))
trip = 0
! Absorbed flux in each layer
do i=snl_top,snl_btm,1
if(i==snl_top) then
F_abs(i) = F_net(i)-F_sfc_net
else
F_abs(i) = F_net(i)-F_net(i-1)
endif
flx_abs_lcl(i,bnd_idx) = F_abs(i)
! ERROR check: negative absorption
if (flx_abs_lcl(i,bnd_idx) < -0.00001) then
trip = 1
endif
enddo
flx_abs_lcl(1,bnd_idx) = F_btm_net
if (flg_nosnl == 1) then
! If there are no snow layers (but still snow), all absorbed energy must be in top soil layer
!flx_abs_lcl(:,bnd_idx) = 0._r8
!flx_abs_lcl(1,bnd_idx) = F_abs(0) + F_btm_net
! changed on 20070408:
! OK to put absorbed energy in the fictitous snow layer because routine SurfaceRadiation
! handles the case of no snow layers. Then, if a snow layer is addded between now and
! SurfaceRadiation (called in Hydrology1), absorbed energy will be properly distributed.
flx_abs_lcl(0,bnd_idx) = F_abs(0)
flx_abs_lcl(1,bnd_idx) = F_btm_net
endif
!Underflow check (we've already tripped the error condition above)
do i=snl_top,1,1
if (flx_abs_lcl(i,bnd_idx) < 0._r8) then
flx_abs_lcl(i,bnd_idx) = 0._r8
endif
enddo
F_abs_sum = 0._r8
do i=snl_top,snl_btm,1
F_abs_sum = F_abs_sum + F_abs(i)
enddo
!ERROR check: absorption greater than incident flux
! (should make condition more generic than "1._r8")
if (F_abs_sum > 1._r8) then
trip = 1
endif
!ERROR check:
if ((albedo_1 < 0._r8).and.(trip==0)) then
write(6,*) 'ERROR: albedo_1 <0 = ', albedo_1
trip = 1
endif
! Set conditions for redoing RT calculation
if ((trip == 1).and.(flg_dover == 1)) then
flg_dover = 2
elseif ((trip == 1).and.(flg_dover == 2)) then
flg_dover = 3
elseif ((trip == 1).and.(flg_dover == 3)) then
flg_dover = 4
elseif((trip == 1).and.(flg_dover == 4).and.(err_idx < 20)) then
flg_dover = 3
err_idx = err_idx + 1
write(6,*) "SNICAR WARNING: Both approximations failed with new zenith angle :(. Zenith= ", mu_not
elseif((trip == 1).and.(flg_dover == 4).and.(err_idx >= 20)) then
flg_dover = 0
write(6,*) "SNICAR ERROR: FOUND A WORMHOLE. STUCK IN INFINITE LOOP! Called from: "
write(6,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
write(6,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
write(6,*) "SNICAR STATS: soot1(0)= ", mss_cnc_aer_lcl
call endrun()
else
flg_dover = 0
endif
enddo !enddo while (flg_dover > 0)
! Energy conservation check:
! Incident direct+diffuse radiation equals (absorbed+bulk_transmitted+bulk_reflected)
energy_sum = (mu_not*pi*flx_slrd_lcl(bnd_idx)) + flx_slri_lcl(bnd_idx) - (F_abs_sum + F_btm_net + F_sfc_pls)
if (abs(energy_sum) > 0.00001_r8) then
write (6,"(a,e12.6,a,i6,a,i6)") "SNICAR ERROR: Energy conservation error of : ", energy_sum
! call endrun()
endif
albout_lcl(bnd_idx) = albedo_1
! Check that albedo is less than 1
if (albout_lcl(bnd_idx) > 1.0) then
write (6,*) "SNICAR ERROR: Albedo > 1.0 "
write (6,*) "SNICAR STATS: bnd_idx= ",bnd_idx
write (6,*) "SNICAR STATS: albout_lcl(bnd_idx)= ",albout_lcl(bnd_idx), " albsfc_lcl(bnd_idx)= ",albsfc_lcl(bnd_idx)
write (6,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
write (6,*) "SNICAR STATS: coszen= ", COSZ, " flg_slr= ", 1
write (6,*) "SNICAR STATS: soot(0)= ", mss_cnc_aer_lcl
write (6,*) "SNICAR STATS: L_snw(-4)= ", L_snw(-4)
write (6,*) "SNICAR STATS: L_snw(-3)= ", L_snw(-3)
write (6,*) "SNICAR STATS: L_snw(-2)= ", L_snw(-2)
write (6,*) "SNICAR STATS: L_snw(-1)= ", L_snw(-1)
write (6,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
call endrun()
endif
enddo ! loop over wvl bands
! Weight output NIR albedo appropriately
albout(1) = albout_lcl(1)
flx_sum = 0._r8
do bnd_idx= nir_bnd_bgn,nir_bnd_end
flx_sum = flx_sum + flx_wgt(bnd_idx)*albout_lcl(bnd_idx)
end do
albout(2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
! Weight output NIR absorbed layer fluxes (flx_abs) appropriately
flx_abs(:,1) = flx_abs_lcl(:,1)
do i=snl_top,1,1
flx_sum = 0._r8
do bnd_idx= nir_bnd_bgn,nir_bnd_end
flx_sum = flx_sum + flx_wgt(bnd_idx)*flx_abs_lcl(i,bnd_idx)
enddo
flx_abs(i,2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
end do
elseif ((h2osno_lcl < min_snw) .and. (h2osno_lcl > 0._r8)) then
albout(1) = UALB
albout(2) = UALB
! There is either zero snow, or no sun
else
albout(1) = 0._r8
albout(2) = 0._r8
endif ! if column has snow and coszen > 0
ALBSND(1)=albout(1)
ALBSNI(1)=albout(1)
ALBSND(2)=albout(2)
ALBSNI(2)=albout(2)
write (8,*) "check the ALBSND", ALBSND(1)
write (8,*) "check the ALBSNI", ALBSNI(2)
END SUBROUTINE SNOWALB_SNICAR_DIR
新函数定义好后，在对应的地方使用：
private :: SNOWALB_SNICAR_DIR
IF(OPT_ALB == 3) THEN
CALL SNOWALB_SNICAR_DIR (parameters,NBAND,NSNOW,ISNOW,COSZ,BCS,UALB,SNICE,SNLIQ,SNEQV,ALBSND,ALBSNI)
END IF
变量添加
变量添加比较麻烦，由于Fortran需要对每一个变量都声明并分配的特点，我们要找出所有与该过程相关的函数，并将对应的变量声明添加上去。
在NoahMPLSM模块中，我主要新增了两个变量 UALB与BCS.。
其中，UALB从背景反照率ALBBCK中获取，BCS则来源于自己添加的变量MBC。
为此，我们需要增加相应的变量声明：
REAL(r8),INTENT(IN) :: BCS !bc IN snow (ng/g)
REAL,INTENT(IN) :: UALB !Albedo of underlyINg surface
注意，该声明涉及所有使用学反照率的函数与模块，我们都需要进行添加。
修改完所有以后，我们在Registry中加入：
state real mbc ij misc 1 - i012rh "MBC" "BC in snow" "ng/g"
这样，MBC就可以通过grid%mbc直接赋值了。
在对应的dyn_em文件夹下，修改module_first_rk_step_part1.F加入mbc变量。
修改phys文件夹所有与noahmp模块相关的模块，包括：
module_sf_noahmplsm.F module_sf_noahmpdrv.F module_surface_driver.F，如果修改noahmp的初始化函数的话，还需修改moudle_phys_init.F。
大致的顺序是：WRF使用first_rk_step读取处理初始化变量，first_rk_step使用surface_driver模块对地面变量初始化，real完成后，WRF会根据选择的物理方案，调用phys_init针对方案对不同方案的各个变量进行初始化，在使用Noah-MP方案时，phys_init会调用Noahmpdrv模块里的NOAHMPINIT函数，noahmp-drv模块则需要使用noahlsm模块中的NOAHMPSHFLX这个函数进行运算，不同模块间都存在相互耦合的关系。
除了上述以外，为了确保MBC被正确读取以及之后进行并行运算，最好修改share下的 module_wps_io_arw.F:
VarName='FLAG_MBC'
call retrieve_index(index,VarName,varname_all,nrecs,iret)
call mpi_file_read_at(iunit,file_offset(index)+5*4, &
igarb,1,mpi_integer4, &
mpi_status_ignore, ierr)
flag_mbc=igarb
以及moule_optional_input：
CALL optional_mbc ( grid , fid , &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
SUBROUTINE optional_mbc ( grid , fid , &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
USE module_io_wrf
USE module_domain , ONLY : domain
USE module_configure , ONLY : grid_config_rec_type
USE module_io_domain
IMPLICIT NONE
TYPE ( domain ) :: grid
INTEGER , INTENT(IN) :: fid
INTEGER :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte
INTEGER :: itmp , icnt , ierr
flag_name = ' '
flag_mbc = 0
flag_name(1:8) = 'MBC '
CALL wrf_get_dom_ti_integer ( fid, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
IF ( ierr .EQ. 0 ) THEN
flag_mbc = itmp
END IF
END SUBROUTINE optional_mbc
来保证初始化的并行运算。
致谢
感谢NCAR的 He cenlin教授、UCLA的薛永康教授以及本所冯锦明研究员的回复与解答。
3423
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
