Physics modules and interface in WRF:


Page reproduced from the WRF original docs page at mmm.ucar.edu: http://www.mmm.ucar.edu/wrf/users/docs/wrf-phy.html

WRF PHYSICS----
History of this documentation
Original : Annual report: WRF physics (Chen and Dudhia 2000)
Modified: 1. WRF 1.2 release (Changes , 04/17/02)


I. INTRODUCTION
II. PHYSICS SCHEME
III. PHYSICS INTERFACE - A THREE-LEVEL PHYSICS STRUCTURE
IV. IMPLEMENTING A NEW PHYSICS SCHEME
V REFERENCES



I. INTRODUCTION
To simulate real weather and to do simulations with coarse resolutions, a minimum set of physics components is required, namely radiation, boundary layer and land-surface parameterization, convective parameterization, subgrid eddy diffusion, and microphysics. Since the model is developed for both research and operational groups, sophisticated physics schemes and simple physics schemes are needed in the model. The objectives of the WRF physics development are to implement a basic set of physics into the WRF model and to design a user friendly physics interface. Since the WRF model is targeted at resolutions of 1-10 km, some of physics schemes might not work properly in this high resolution (e.g. cumulus parameterization). However, at this early stage of model development, only existing physics schemes are implemented, and most of them are taken from current mesoscale and cloud models. In the future, new physics schemes designed for resolutions of 1-10 km should be developed and implemented.
II. PHYSICS SCHEMES
Currently, several physics components have been included in WRF: microphysics, cumulus parameterization, long wave radiation, short wave radiation, boundary layer turbulence (PBL), surface layer, land-surface parameterization, and subgrid scale diffusion. The rest of this section will briefly introduce the schemes, which have been implemented.
2.1 Microphysics
(a) Kessler (kesslerscheme) [module_mp_kessler.F]Kessler scheme (Kessler 1969), which was taken from COMMAS, is a simple warm cloud scheme, and it includes water vapor, cloud water, and rain. The microphysical processes included are: the production, fall and evaporation of rain, the accretion and autoconversion of cloud water, and the production of cloud water from condensation.


(b) Purdue Lin scheme (linscheme) [module_mp_lin.F]
Six classes of hydrometeos are included : water vapor, cloud water, rain, cloud ice, snow, and graupel. All parameterization production terms are based on Lin et al. (1983) and Rutledge and Hobbs (1984) with some modifications, including saturation adjustment following Tao (1989) and ice sedimentation. This is a relatively sophisticated microphysics scheme in WRF, and it is more suitable for use in research studies. The scheme is taken from Purdue cloud model and the details can be found in Chen and Sun (2002).


(c) NCEP simple ice (ncepcloud3) [module_mp_ncloud3.F]
This scheme follows Hong et al. (1998) with some modification, now including the ice sedimentation effect. Three categories of hydrometers are included: vapor, cloud water/ice, and rain/snow. The cloud ice and cloud water are counted as the same category, and they are distinguished by temperature, namely cloud ice can only exist when the temperature is less than or equal to freezing point, otherwise cloud water can exist. The same condition is applied to rain and snow. Though the ice phase is included, it is considered simple enough for using in operational models.


(d) NCEP mixed phase (ncepcloud5) [module_mp_ncloud5.F]
This scheme is similar to NCEP simple ice scheme. However, rain and snow (cloud ice and cloud water) are in different categories. It allows supercooled water to exist, and gradual melting of snow as it falls. Ice sedimentation is now included. Details can be found in Hong et al. (1998).


(e) Eta microphysics (etampscheme) [module_mp_eta.F]
The scheme explicitly predicts the cloud water/ice mixing ratio as described in Zhao and Carr (1997). Liquid and frozen precipitation are derived diagnostically from the cloud mixing ratio and are assumed to fall to the ground in a single time step. This scheme is implemented to conduct tests with old Eta model. (This paragraph was provided by Tom Black.)

(f)
Eta Grid-scale Cloud and Precipitation scheme of 2001 (current name "etampnew", future name "EGCP01") [module_mp_etanew.F]
The scheme predicts changes in water vapor and total condensate that are advected in the model. Local storage arrays retain first-guess information that extract contributions of cloud water, rain, cloud ice, and precipitation ice of variable density in the form of snow, graupel, or sleet. The density of precipitation ice is estimated from a local array that stores information on the total growth of ice by vapor deposition and accretion of liquid water. Sedimentation is treated by partitioning the time-averaged flux of precipitation into a grid box between local storage in the box and fall out through the bottom of the box. This approach, together with modifications in the treatment of rapid microphysical processes, permits large time steps to be used with stable results. The mean size of precipitation ice is assumed to be a function of temperature following the observational results of Ryan (1996). Mixed-phase physics are considered at temperatures warmer than -10C, whereas ice saturation is assumed for cloudy conditions at colder temperatures. Further description of the scheme can be found in Sec. 3.1 of the November 2001 Technical Procedures Bulletin (TPB) at http://www.emc.ncep.noaa.gov/mmb/mmbpll/eta12tpb/ and on the COMET page at http://meted.ucar.edu/nwp/pcu2/etapcp1.htm. (This paragraph was provided by Brad Ferrier)
Currently, the sedimentation process, which could be combined with the calculation of vertical velocities in the future, is counted inside microphysics, and a smaller time step is allowed to calculate the vertical flux of precipitation to prevent instability. The saturation adjustment is also included inside the microphysics, and it might be separated as an individual subroutine in the future.


2.2 Convective schemes

(a) New Kain-Fritsch (kfetascheme) [module_cu_kfeta.F]
The modified version of the Kain-Fritsch scheme (KF-Eta) is based on Kain and Fritsch (1990, 1993), but has been modified based on testing within the Eta model. As with the original KF scheme, it utilizes a simple cloud model with moist updrafts and downdrafts, including the effects of detrainment, entrainment, and relatively crude microphysics. It differs from the original KF scheme in the following ways:



  • A minimum entrainment rate is imposed to suppress widespread convection in marginally unstable, relatively dry environments
  • Shallow (non precipitating) convection is allowed for any updraft that does not reach minimum cloud depth for precipitating clouds; this minimum depth varies as a function of cloud-base temperature.
  • Entrainment rate is allowed to vary as a function of low-level convergence.
  • Downdraft changes:

  • source layer is entire 150 - 200 mb deep layer just above cloud base.
  • mass flux is specified as a fraction of updraft mass flux at cloud base; fraction is a function of source layer RH rather than wind shear or other parameters, i.e., old precipitation efficiency relationship not used.
  • detrainment is specified to occur in updraft source layer and below.


(This paragraph is provided by Jack Kain.)
(b) Bett-Miller-Janjic (bmjscheme) [module_cu_bmj.F]
The scheme is based on the Betts-Miller convective adjustment scheme (Betts 1986; Betts and Miller 1986). Primary modifications were made by Janjic (1990, 1994, 2000) including the introduction of "cloud efficiency" to provide an additional degree of freedom in determining target profiles of heat and moisture. Shallow convective adjustment is an important part of the parameterization. (This paragraph was provided by Tom Black.)
(c) Kain-Fritsch (kfscheme) [module_cu_kf.F]The Kain-Fritsch scheme (Kain and Fritsch 1990, 1993) is taken from the MM5 model, and a simple cloud model, which includes detrainment, entrainment, updraft, and downdraft, is incorporated into this scheme. This version only includes deep convection. This scheme is implemented to conduct tests with old Eta model.
2.3 Long wave radiation
(a) RRTM (rrtmscheme) [module_ra_rrtm.FThis RRTM, which is taken from MM5, is based on Mlawer et al. (1997), and it is a spectral-band scheme using the correlated-k method. It uses pre-set tables to accurately represent longwave processes due to water vapor, ozone, CO2, and trace gases (if present) as well as accounting for cloud optical depth.


(b) ETA GFDL long wave (gfdllwscheme) [module_ra_gfdleta.F]
This longwave radiation is from GFDL. It follows the simplified exchange method of Fels and Schwarzkopf (1975) and Schwarzkopf and Fels (1991), with calculation over spectral bands associated with carbon dioxide, water vapor, and ozone. Schwarzkopf and Fels (1985) transmission coefficients for carbon dioxide, a Roberts et al. (1976) water vapor continuum, and the effects of water vapor-carbon dioxide overlap and of a Voigt line-shape correction are included. The Rodgers (1968) formulation is adopted for ozone absorption. Clouds are randomly overlapped. (This paragraph was provided by Kenneth Campana.) This scheme is implemented to conduct tests with old Eta model.

2.4 Short wave radiation
(a) Simple short wave (swscheme) [module_ra_sw.F]The scheme is base on Dudhia (1989) and taken from MM5. It is a simple downward integration of solar flux, accounting for clear-air scattering, water vapor absorption (Lacis and Hansen 1974), and cloud albedo and absorption. It uses look-up tables for clouds from Stephens (1978).


(b) Goddard short wave (gsfcswscheme) [module_ra_gsfcsw.F]
This scheme is base on Chou and Suarez (1994).


(c) ETA GFDL short wave (gfdlswscheme) [module_ra_gfdleta.F]
This shortwave radiation is a GFDL version of the Lacis and Hansen (1974) parameterization. Effects of atmospheric water vapor, ozone (both from Lacis and Hansen, 1974) and carbon dioxide (Sasamori, et al., 1972) are employed. Clouds are randomly overlapped. SW calculations are made using a daylight-mean cosine solar zenith angle over the time interval (given in namelist.input). The diurnal cycle is approximated during the hour by weighting SW heating rates and fluxes with the actual cosine zenith angle at each model time-step. (This paragraph was provided by Kenneth Campana.) This scheme is implemented to conduct tests with old Eta model.


2.5 Surface layer
(a) Similarity theory (sfclayscheme) [module_bl_sfclay.F]Uses stability functions from Paulson (1970), Dyer and Hicks (1970), and Webb (1970) to compute surface exchange coefficients for heat, moisture, and momentum.


(b) MYJ surface scheme (myjsfcscheme) [module_bl_myjsfc.F]
The surface layer scheme (Janjic 1996a, 2002) is based on the similarity theory (Monin and Obukhov, 1954). The scheme includes parameterizations of viscous sublayer. Over water surfaces the viscous sublayer is parameterized explicitly following Janjic (1994). Over land, the effects of the viscous sublayer are taken into account through variable roughness height for temperature and humidity as proposed by Zilitinkevitch (1995). The Beljaars (1994) correction is applied in order to avoid singularities in the case of unstable surface layer and vanishing wind speed. The surface fluxes are computed by an iterative method. (This paragraph was provided by Zavisa Janjic)

2.6 Land-surface
(a) Thermal diffusion (slabscheme) [module_bl_slab.F]This is based on the MM5 5-layer soil temperature model. Layers are 1, 2, 4, 8, and 16 cm thick. Below these the temperature is fixed at diurnal average. Energy budget includes radiation, sensible, and latent heat flux. It also allows for snow-cover flag.


(b) OSU/MM5 scheme (lsmscheme) [module_bl_lsm.F]
Described by Chen and Dudhia (2000). This is a 4-layer soil temperature and moisture model with canopy moisture and snow cover prediction. It includes root zone, and other vegetation effects, drainage, and runoff. Provides sensible and latent heat fluxes to boundary-layer scheme.

2.7 Boundary layer parameterization
(a) MRF (mrfscheme) [module_bl_mrf.F]The scheme is described by Hong and Pan (1996). This uses a so-called countergradient flux for heat and moisture in unstable conditions. It uses enhanced vertical flux coefficients in the PBL, and the PBL height is determined from a critical bulk Richardson number. It handles vertical diffusion with an implicit local scheme, and it is based on local Ri in the free atmosphere.


(b) MYJ (myjpblscheme) [module_bl_myjpbl.F]
The parameterization of turbulence in the PBL and in the free atmosphere (Janjic, 1990, 1996b, 2002) represents a nonsingular implementation of the Mellor-Yamada Level 2.5 turbulence closure model (Mellor and Yamada, 1982) through the full range of atmospheric turbulent regimes. In this implementation, an upper limit is imposed on the master length scale. This upper limit depends on TKE and the buoyancy and shear of the driving flow. In the unstable range the functional form of the upper limit is derived from the requirement that the TKE production be nonsingular in the case of growing turbulence. In the stable range the upper limit is derived from the requirement that the ratio of the variance of the vertical velocity deviation and TKE cannot be smaller than that corresponding to the regime of vanishing turbulence. The TKE production/dissipation differential equation is solved iteratively The empirical constants have been revised as well (Janjic, 1996b, 2002). (This paragraph was provided by Zavisa Janjic)

2.8 Subgrid eddy diffusion
(a) Simple diffusion (diff_opt=1)A simple diffusion is built using K theory, where K is constant. This option comes with the choice of diff_opt =1 in the namelist.input file. The momentum eddy coefficients are khdif for the horizontal and kvdif for the vertical. The heat eddy coefficients are three times of momentum eddy coefficient.


(b) Stress/deformation form (diff_opt=2) [dyn_eh or dyn_em/module_diffusion.F ]
The subgrid diffusion equations with stress formulae are derived to match the WRF governing equations (Chen and Dudhia 2000). Three eddy diffusion coefficients are available: constant, function of turbulence kinetic energy, and function of deformation, and the option flag (km_opt) is in the namelist.input file. The constant coefficients (km_opt =1) takes the values of khdif and kvdif as in (a). The eddy coefficients in terms of TKE (km_opt=2) or Smagorinsky (km_opt=3) , as well as the TKE equations, mainly follow ARPS model, Moeng (1984), and Moeng and Wyngaard (1989). This code is newly written for the WRF model. Note that the vertical subgrid diffusion is only called if there is no boundary layer scheme.

III. PHYSICS INTERFACE - A THREE-LEVEL PHYSICS STRUCTURE
The WRF model does not allow the use of common blocks and therefore, all variables have to be passed into subroutines through argument lists. The use of the modules (a feature of FORTRAN 90) makes the design of physics interface easier. A three-level structure (solver, driver, and individual scheme) of physics (Figure 1) has been constructed. The idea of three-level structure is so that users can easily implement their schemes into the WRF with little need to understand other parts of WRF code.external image wrf_phy_fig1.gif

Figure 1: A three-level physics code structure.

3.1 Level 1 - Solver
The first level, solver (solve_eh.F for height model and solve_em.F for mass model), is the main subroutine for calling dynamics and physics as the schematic diagram shows in Figure 2. The solver is a bridge, which connects dynamics and physics. In addition, the parallel processes regarding the model part, such as multi threading and message passing, are also dealt with in this subroutine.
external image wrf_phy_fig2.gif

Figure 2: A schematic diagram of the WRF model. 'wrf' is the main program and solve_em is the major subroutine for mass model solver, which calls dynamics and physics. Boxes in dash-line and double-line (individual schemes) represent physics codes.

3.1 Level 1 - Solver
3.2 Level 2 - Driver
Instead of calling different physics schemes directly, the solver calls the physics drivers, which are the interfaces between the solver and individual physics schemes, and shared by the two dynamic cores. The mechanism that allows different dynamic cores to use the same physics driver is the physics preparation calls in each solver (phy_prep and moist_physics_prep). This subroutine calculates variables used in all physics routines. Except subgrid eddy diffusion, every physics component has its own driver, such as the microphysics_driver, cumulus_driver, pbl_driver, and radiation_driver. The pbl_driver includes the surface layer, land/sub-surface, and boundary layer, and the radiation_driver includes long wave and short wave radiations. In each driver, CASE SELECT statements are used to branch to the choice of scheme. The flags for choosing different physics schemes are located in the namelist.input file. The design of the driver interface is to provide a plug-compatible environment for users to plug in their own physics packages.
3.3 Level 3 - Individual physics scheme.
The prohibition of using common blocks separates users' codes from the WRF code quite well, and it allows users to easily implement their packages into WRF, though some minor changes of the WRF code might be required. One module comes with each physics scheme, and all the subroutines related to the scheme should be included into this module. The details are given in Section 4.

IV. IMPLEMENTING A NEW PHYSICS SCHEME
This section provides some information about the coding rules of WRF physics, and the procedure of implementing a new physic scheme in detail.
4.1 Naming rules
Except for subgrid scale diffusion, each physics scheme has its own module, and the module is named as module_yy_xxx.F. Here yy represents different physics components and it is defined as follows:



ra is for radiation parameterization,bl is for PBL parameterization,cu is for cumulus parameterization, andmp is for microphysics.

and xxx is named to match each individual scheme. For example, module_cu_kf.F is a cumulus parameterization module, and it includes all subroutines related to the old Kain-Fritsch scheme.
Some variables are also named corresponding to different physics components, such as physics tendencies and the frequencies of calling different physics components. As seen in the WRF code, physics tendencies are named as zzzyyTEN (e.g. RTHRATEN, RUBLTEN, RQCCUTEN) where zzz is RU, RV, RTH, RQV, RQC, RQR, etc. STEPyy (e.g. STEPCU, STEPRA, STEPCL) is the frequency, in terms of model time steps, of calling physics.
(Note: the physics tendency arrays will remain constant, and add to the other tendencies, at each model timestep until the physics is called again after STEPxx timesteps. Microphysics and subgrid eddy diffusion have no corresponding STEPyy since they are called every time step. To change the frequencies of calling physics such as radiation, cumulus, and PBL, please go to file namelist.input and modify the values of RADT (for radiation), BLDT (for PBL), and CUDT (for cumulus). They are in the unit of minutes.)


4.2 Coding rules
(a)
As mentioned earlier, no common block is allowed in the WRF model if the scheme is to be considered for the official version. All variables used in subroutines have to pass through argument lists.

(b)
FORTRAN 90 is the official language used in the model part and it has some C++ features, such as pointers, user defined data type, and modules.

(c)
Before calculating physics, all variables or decoupled variables (e.g. u and v) have to be horizontally interpolated to the column of mass fields (A grid) if it is necessary, but the vertical ones (e.g. w) can be still staggered. This also means that some calculated physics tendencies (e.g. RUBLTEN and RVBLTEN) need to be interpolated back to horizontal staggered grid later on. Another attention has to paid to the vertical index order between your code and the WRF model. In WRF, kms (smallest number) is the bottom level and kme (largest number) is the top level. In a use's scheme, if 1 is at the top level, then you have to reverse the order in the k direction. The following is the vertical structure of the WRF. Please also see 4.4.a for more information about indices.


kme - H (no data at this level)kme ----- Fkme-1 - Hkme-1 ----- F . . .kms+2 - Hkms+2 ----- Fkms+1 - Hkms+1 ----- Fkms - H (half level, u, v, w, p, T, q, etc..)kms ----- F (full level, w)

4.3 Implementing a new scheme
Assume that we are interested in implementing a new cumulus parameterization scheme such as the Arakawa scheme, into the WRF model. The procedure is as follows: (Note: In the following procedure, the names of modules, files and the pieces of WRF code are italic and in light blue color. In addition, new code regarding the new scheme is bolded.)

(a)
Create a new cumulus module module_cu_arakawa.F(The pseudo codes of different physics components are available in Appendix B in Chen and Dudhia 2000). All subroutines related to the new scheme should be included inside this module, such as the scheme itself and its initialization. The scheme-related local constants, which won't be used in other schemes, have to be placed at the top of the module.

(b)
Declare a new package in the file Registry, which is in the directory Registry. For example,

package kfscheme cu_physics==1 - -package arakawascheme cu_physics==4 - -

Package is a new defined data type. The first line is the declaration of an existing Kain-Fritsch scheme, and the second bolded line is the declaration of the new scheme package. In the model, arakawascheme corresponds to the number of 4 and the string arakawascheme should be used to present the scheme itself (examples are given later) instead of using the number 4. Currently, there are seven different categories of physics packages defined in the model:
cu_physics
cumulus parameterization
mp_physics
microphysics
ra_lw_physics
long wave radiation
ra_sw_physics
short wave radiation
bl_sfclay_physics
surface layer
bl_surface_physics
land-surface layer
bl_pbl_physics
boundary layer

The declaration of the mp_physics is slightly different from the others. For example, the Kessler scheme package in the Registry is defined as:

package kesslerscheme mp_physics==1 - moist:qv,qc,qr

The fifth column (moist:qv,qc,qr) contains the moisture fields required in this scheme. This also implies that only these moisture fields are active (memory is allocated) and are available for use in the other physics parameterization schemes, such as cumulus parameterization, radiation, and pbl parameterization. ||

(c)
Reassign cu_physics = 4 in the file namelist.input in order to choose the new scheme. After compiling the model, the schemes corresponding to different numbers can be found in the file module_state_description.F. The choices of different schemes for different physics components are also in the namelist.input such as, mp_physics, ra_lw_physics, ra_sw_physics, bl_sfclay_physics, bl_surface_physics, bl_pbl_physics, cu_physics. Their definitions are the same as those in the file Registry.

(d)
Initialize initial values for the new scheme inside module_physics_init.F if necessary. Otherwise, skip this step.
a. Go to subroutine cu_init in module_physics_init.Fb. Add a new case statement in SELECT CASE, such as:


cps_loop: SELECT CASE(config_flags%cu_physics)
CASE (KFSCHEME)
CALL kfinit(...)......
CASE (ARAKAWASCHEME)
CALL arakawainit(...) include subroutine arakawainit inside the module of the new scheme (module_cu_arakawa.F)
CASE DEFAULT
ENDSELECT cps_loop

The subroutine arakawainit should be included inside the module of the new scheme (module_cu_arakawa.F). As mentioned earlier, the name "ARAKAWASCHEME" has to be used (with case free) instead of the number, and it spells exactly the same as the one in the Registry (the second column of the package declaration).

(e)
Make the change in the driver cumulus_driver.F (microphsics_driver.F, for microphysics, pbl_driver.F for PBL, and radiation_driver.F for radiation).
(i) Include the new module into the driver such as:


USE module_cu_kfUSE module_cu_arakawa <--- new scheme

(ii) Check the available variables passed through argument list. The meanings and units of variables are explained inside each driver. If there is any variable, which is not available in this driver, you need to go one level up (solve_eh.F/solve_em.F) and pass the variable through argument list. Meanwhile, you also need to make the same change to the interface file physics_driver.int. If the new required variable is not available in the model, you have to declare a new variable in the Registry (it will magically show up in the file solve_eh.F/solve_em.F) or just use any inactive working space/variable and then, calculate it in the subroutine phy_prep, which is called by solve_eh.F/solve_em.F, and pass it to driver.

(iii) Check all global constants defined in the module_model_constants.F. If there is any global constant, which is not defined yet, please add the new one into the file module_model_constants.F and it will show up in the driver (cumulus_driver.F) subroutines since the module module_model_constants is included inside the drivers. Only use this for constants that might be needed in more than one physics package. Otherwise, define the constants at the top of the new physics module.
(iv) Create a new SELECT CASE in cumulus_driver.F


cps_loop: SELECT CASE(config_flags%cu_physics)
CASE (KFSCHEME)
CALL KFCPS(...)
CASE (ARAKAWASCHEME) <--- new scheme
CALL arakawa(...)

include subroutine arakawa inside the module of the new scheme (module_cu_arakawa.F)

CASE DEFAULT
END SELECT cps_loop

Again, the name "ARAKAWASCHEME" should match the one in the Registry (second column) you have defined.

(f)
Check physics tendencies currently exist in the file Registry. Create new tendency variables in the Registry if it is necessary. They will show up in the file solve_eh.F/solve_em.F.

(g)
Update physics tendencies to diffusion tendencies (zzz_tendf)

If there is any new created or calculated tendencies, these tendencies have to be passed into subroutine update_phy_ten (called by solve_eh.F/solve_em.F), and passed down to the bottom of the calling tree of subroutines such as phy_cu_ten for cumulus parameterization (phy_ra_ten for radiation and phy_bl_ten for PBL). Then, create a new case statement in SELECT CASE, such as:

SELECT CASE (config_flags%cu_physics)
CASE (KFSCHEME)
CALL add_a2a(rt_tendf,RTHCUTEN, ...)CALL add_a2a(moist_tendf(ims,kms,jms,P_QV),RQVCUTEN, ...)IF (P_Qx .ge. PARAM_FIRST_SCALAR) &
CALL add_a2a(moist_tendf(ims,kms,jms,P_Qx),RQxCUTEN, ...)
CASE(ARAKAWASCHEME)
CALL add_a2a(rt_tendf,RTHCUTEN, ...)CALL add_a2a(moist_tendf(ims,kms,jms,P_QV),RQVCUTEN, ...)
IF (P_Qx .ge. PARAM_FIRST_SCALAR) &
CALL add_a2a(moist_tendf(ims,kms,jms,P_Qx),RQxCUTEN, ...)
(Note: P_Qx is P_QC (cloud water), P_QR (rain), P_QI (ice), P_QS (snow), or P_QG (graupel))
CASE DEFAULT
END SELECT
It is important to test the existence of any microphysical variables (e.g. IF (P_Qx .ge. PARAM_FIRST_SCALAR)) before using that particular microphysical variable and its tendency.
Every physics tendency in the new scheme has to accumulate into the non-advective tendencies such as rt_tendf, moist_tendf, and so on, by means of a series of calls (CALL add_a2a, CAll add_a2c_u, or CALL add_a2c_v). The subroutine add_a2a is used to add the physics tendencies into the non-advective tendencies on the same grid points such as moisture, temperature, and vertical velocity fields. The subroutine add_a2c_u (add_a2c_v) is used to add the u (v) physics tendencies, which is on A grid, to u (v) non-advective tendencies, which is on C grid. Please also read 4.4d and 4.4g for other details.
It is also important to know that the P_QV, P_QC, P_QR, P_QI, P_QS, or P_QG have be used instead of numbers (2,3,...). In any case, a loop over microphysical variables is required for any physics scheme. The loop has to run from PARAM_FIRST_SCALAR (not 1) to number of microphysical variables. Again, the num_3d_m (the name could be changed after passing into physics driver) instead of a number has to be used for the number of microphysical variables.
(Note: Currently microphysics is treated differently from other physics. The other physics drivers are called before the sound-wave short steps to allow diabatic terms to be used for theta on the short steps. The microphysics driver is called at the end of the time step after all the variables have been updated. This means that the microphysics directly updates the variables and has no associated tendencies. This is done at the end to allow an accurate saturation adjustment, and the diabatic term is saved for the next time-step's short steps as an approximation. [In the future we could separate the saturation adjustment from the rest of the microphysics so that we can have the microphysics treated like the other physics.])



(h)
Add the new module into the Makefile which is in the directory src and modify the dependencies.

4.4 General notes
(a)
The example codes show three sets of indices such as

...ids,ide, jds,jde, kds,kde, &ims,ime, jms,jme, kms,kme, &its,ite, jts,jte, kts,kte )...

The first set (ids, ...) refers to domain starting and ending indices, and should only be used for boundary checks (not often needed in the physics). The second set (ims, ...) refers to memory dimensions and are distinct from the first set to enable distributed memory application. Arrays passed in from the solver are dimensioned by these indices. The third set (its, ...) refers to tile dimensions, and differs from the domain dimensions when there is multiple threading. Local arrays inside the physics generally need only to be dimensioned by the tile size, e. g. (its: ite, kts: kte), if the physics routine is called inside a j loop (do j= jts, jte). All loops should use these limits for i, j, and k, and oper-ations should not affect array elements outside these limits.

(b)
The index order for WRF is IKJ. Physics code can conform to this by calling the lower level routines inside a J loop and having local variables dimensioned (i, k), and then having an inner i-loop allowing vectorization (see XXX2D example for pbl scheme). Column physics is also allowed but discouraged because of its usually poor vector performance and non-compliance with the WRF storage and loop-order standards. It requires local (k) arrays in the physics routines (XXX1D examples). Array syntax is used in only a few limited situations.

(c)
Note that module_ model_ constants contains many of the basic physical constants. The physics driver uses this module and passes selected constants to the physics routines. It is recommended to use these constants as much as possible, and not to define alternative values inside the physics module. This also applies to the saturation vapor pressure calculation.

(d)
The Registry can also be used to add variables required at the solver level. For instance, if you have a cumulus scheme with momentum transport, you will have to add RUCUTEN and RVCUTEN to the Registry, and make additions to the routine phy_ cu_ ten. Recall also that momentum tendencies in the physics are calculated on the A-grid, and so they need to be interpolated back to the C-grid (e. g. see phy_ bl_ ten) afterwards.

(e)
The Registry can also be used to add species for moisture and chemistry. Moisture variables are carried in WRF as a 4d array (i, k, j, im) where im denotes species number. This is general enough to add any number of species. Note that any new additions in moisture variable will be advected and diffused and will contribute to total density. The 4d array exists down to the physics driver level and is separated into 3d moisture arrays such as moist( ims, kms, jms, P_ QV) by use of the P_ Qx constants in the CALL from the driver. P_ QV, P_ QC, P_ QR, P_ QI, P_ QS, P_ QG indicate water vapor, cloud water, rain, cloud ice, snow, and graupel, respectively.

(f)
The 4d moisture array is not suitable for number concentration or "moment" variables, but these can possibly use the chemical scalar array, or we may add a general scalar array. Turbulent kinetic energy and tracers would also be examples of scalars that do not contribute to mass.

(g)
Horizontally dependent physics such as diffusion has some special considerations. All operations should be on variables within the tile. They can use values outside the tile, but it should be ensured that those values are accounted for in the distributed-memory exchanges of "halo" zones around each patch. This will need expert consultation for specific cases that go beyond the halo zones WRF currently has.
V. REFERENCES
· Beljaars, A.C.M., 1994: The parameterization of surface fluxes in large-scale models under free convection. Quart. J. Roy. Meteor. Soc., 121, 255-270.
· Betts, A. K., 1986: A new convective adjustment scheme. Part I: Observational and theoretical basis. Quart. J. Roy. Meteor. Soc., 112, 677-691.
· Betts, A. K., and M. J. Miller, 1986: A new convective adjustment scheme. Part II: Single column tests using GATE wave, BOMEX, and arctic air-mass data sets. Quart. J. Roy. Meteor. Soc., 112, 693-709.
· Chen, F., and J. Dudia, 2000: Coupling an advanced land-surface/ hydrology model with the Penn State/ NCAR MM5 modeling system. Part I: Model description and implementation. Mon. Wea. Rev., in press.
· Chen, S.-H., and W.-Y. Sun, 2002: A one-dimensional time dependent cloud model. J. Meteor. Soc. Japan , 80, 99-118.
· Chen, S.-H., and J. Dudhia, 2000: Annual report: WRF physics, Air Force Weather Agency, 38pp.
· Chou M.-D., and M. J. Suarez, 1994: An efficient thermal infrared radiation parameterization for use in general circulation models. NASA Tech. Memo. 104606, 3, 85pp.
· Dudia, J., 1989: Numerical study of convection observed during the winter monsoon experiment using a mesoscale two-dimensional model. J. Atmos. Sci., 46, 3077-3107.
· Dyer, A. J., and B. B. Hicks, 1970: Flux-gradient relationships in the constant flux layer. Quart. J. Roy. Meteor. Soc., 96, 715-721.
· Hong, S.-Y., and H.-L. Pan, 1996: Nonlocal boundary layer vertical diffusion in a medium-range forecast model. Mon. Wea. Rev., 124, 2322-2339.
· Hong, S.-Y., H.-M. H. Juang, and Q. Zhao, 1998: Implementation of prognostic cloud scheme for a regional spectral model, Mon. Wea. Rev., 126, 2621-2639.
· Janjic, Z. I., 1990: The step-mountain coordinate: physical package. Mon. Wea. Rev., 118, 1429-1443.
· Janjic, Z. I., 1994: The step-mountain eta coordinate model: further developments of the convection, viscous sublayer and turbulence closure
schemes. Mon. Wea. Rev., 122, 927-945.
· Janjic, Z. I., 1996a: The surface layer in the NCEP Eta Model. Eleventh Conference on Numerical Weather Prediction, Norfolk, VA, 19-23 August 1996; Amer. Meteor. Soc., Boston, MA, 354-355.
· Janjic, Z. I., 2000: Comments on "Development and Evaluation of a Convection Scheme for Use in Climate Models." Journal of the Atmospheric Sciences, Vol. 57, p. 3686.
· Janjic, Z. I., 2002: Nonsingular Implementation of the Mellor-Yamada Level 2.5 Scheme in the NCEP Meso model. NCEP Office Note No. 437, 61 pp.
· Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric circula-tion. Meteor. Monogr., No. 32, Amer. Meteor. Soc., 84pp.
· Kain, J. S., and J. M. Fritsch, 1990: A one-dimensional entraining/ detraining plume model and its application in convective parameterization. J. Atmos. Sci., 47, 2784-2802.
. Kain, J. S., and J. M. Fritsch, 1993: Convective parameterization for mesoscale models: The Kain-Fritcsh scheme. The representation of cumulus convection in numerical models, K. A. Emanuel and D.J. Raymond, Eds., Amer. Meteor. Soc., 246 pp.
·
Lacis, A. A., and J. E. Hansen, 1974: A parameterization for the absorption of solar radiation in the earth's atmosphere. J. Atmos. Sci., 31, 118-133.
· Lin, Y.-L., R. D. Farley, and H. D. Orville, 1983: Bulk parameterization of the snow field in a cloud model. J. Climate Appl. Meteor., 22, 1065-1092.
· Meng, C.-H., 1984: A large-eddy-simulation model for the study of planetary boundary-layer turbulence. J. Atmos. Sci., 41, 2052-2062.
· Meng, C.-H., and J. C. Wyngaard, 1989: Evaluation of turbulence and dissipation closures in second-oeder modeling. J. Atmos. Sci., 46, 23112330.
·
Mlawer, E. J., S. J. Taubman, P. D. Brown, M. J. Iacono, and S. A. Clough, 1997: Radiative trans-fer for inhomogeneous atmosphere: RRTM, a validated correlated-k model for the long-wave. J. Geophys. Res., 102( D14), 16663-16682.
· Monin, A.S. and A.M. Obukhov, 1954: Basic laws of turbulent mixing in the surface layer of the atmosphere. Contrib. Geophys. Inst. Acad. Sci. USSR, (151), 163-187 (in Russian).
·
Paulson, C. A., 1970: The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J. Appl. Meteor., 9, 857-861.
· Rutledge, S. A., and P. V. Hobbs, 1984: The mesoscale and microscale structure and organization of clouds and precipitation in midlatitude cyclones. XII: A diagnostic modeling study of precipitation development in narrow cloud-frontal rainbands. J. Atmos. Sci., 20, 2949-
2972.
· Stephens, G. L., 1978: Radiation profiles in extended water clouds. Part II: Parameterization schemes. J. Atmos. Sci., 35, 2123-2132.
· Tao, W.-K., 1989: An ice-water saturation adjustment. Mon. Wea. Rev., 117, 231-235. Webb, E. K., 1970: Profile relationships: The log-linear range, and extension to strong stability. Quart. J. Roy. Meteor. Soc., 96, 67-90.
· Zilitinkevich, S.S., 1995: Non-local turbulent transport: pollution dispersion aspects of coherent structure of convective flows. In: Air
Pollution III - Volume I. Air Pollution Theory and Simulation (Eds. H. Power, N. Moussiopoulos and C.A. Brebbia). Computational Mechanics Publications, Southampton Boston, 53-60.