# Presentation: 592551

## Tutorial Notes: WRF Software

- John Michalakes, Dave Gill
- NCAR
- WRF Software Architecture Working Group

## Outline

- Introduction
- Computing Overview
- Software Overview
- Examples

## Resources

- WRF project home page
- http://www.wrf-model.org
- WRF users page (linked from above)
- http://www.mmm.ucar.edu/wrf/users
- On line documentation (also from above)
- http://www.mmm.ucar.edu/wrf/WG2/software_v2
- WRF users help desk
- wrfhelp@ucar.edu

## Introduction

- Characteristics of WRF Software
- Developed from scratch beginning around 1998
- Requirements emphasize flexibility over a range of platforms, applications, users; performance
- WRF develops rapidly. First released Dec 2000; Current Release WRF v3.0  (March 2008)
- Supported by flexible efficient architecture and implementation called the WRF Software Framework

## WRF Software Framework Overview

- Implementation of  WRF Architecture
- Hierarchical organization
- Multiple dynamical cores
- Plug compatible physics
- Abstract interfaces (APIs) to external packages
- Performance-portable
- Designed from beginning to be adaptable to today’s computing environment for NWP

## WRF Supported Platforms

- Petascale precursor systems

**Notes:**
This slide shows the systems WRF has been ported to and, in most cases is running on out in the user community.
A few recent developments: Experimental implementation on IBM BG/L
Also Apple G5 workstation ports, but not yet distributed memory parallel
Many variations of UNIX, both by the major vendors and also by a number of Integrators and Turnkey systems, as well as the home grown research dept. varieties.
May have said this the last time I was here but Linux architecture is probably widest distribution in terms of number of CPUs running WRF and also probably the most vexing in terms of community model support.

## WRF Supported Platforms

- Petascale precursor systems

**University of São Paulo **
**“****Clothesline Computer”**

**Notes:**
This slide shows the systems WRF has been ported to and, in most cases is running on out in the user community.
A few recent developments: Experimental implementation on IBM BG/L
Also Apple G5 workstation ports, but not yet distributed memory parallel
Many variations of UNIX, both by the major vendors and also by a number of Integrators and Turnkey systems, as well as the home grown research dept. varieties.
May have said this the last time I was here but Linux architecture is probably widest distribution in terms of number of CPUs running WRF and also probably the most vexing in terms of community model support.

The U. Sao Paulo has a unique low-cost supercomputer installation. See attached. They have 2.2 GhZ IA32 mother boards and power supplies attached to 1'x4' wooden planks, 4 motherboards per plank, suspended with window sash cord from their space-limited computer room ceiling. They call it a "clothesline" computer (actually, the Brazilian word for clothesline). The planks are suspended from pullies and can be individually lowered and raised for maintenance and cleaning. The blue cables are GigE network over which they run MPICH. Very ingenious (assuming their machine room is beaver- and termite-proof).

## Computing Overview

- APPLICATION

## Computing Overview

- APPLICATION

- SYSTEM

## Computing Overview

- APPLICATION

- HARDWARE

- SYSTEM

## Hardware: The Computer

- The ‘N’ in NWP
- Components
- Processor
- A program counter
- Arithmetic unit(s)
- Some scratch space (registers)
- Circuitry to store/retrieve from memory device
- Memory
- Secondary storage
- Peripherals
- The implementation has been continually refined, but the basic idea hasn’t changed much

## Hardware has not changed much...

## Hardware has not changed much...

*~50,000 flop/s*
*48hr 12km WRF CONUS in 600 years*

## Hardware has not changed much...

*~50,000 flop/s*
*48hr 12km WRF CONUS in 600 years*

## Hardware has not changed much...

*~50,000 flop/s*
*48hr 12km WRF CONUS in 600 years*

*~5,000,000,000 flop/s*
*48 12km WRF CONUS in 52 Hours*

## ...how we use it has

- Fundamentally, processors haven’t changed much since 1960
- Quantitatively, they haven’t improved nearly enough
- 100,000x increase in peak speed
- > 4,000x increase in memory size
- These are too slow and too small for even a moderately large NWP run today
- We make up the difference with _parallelism_
- Ganging multiple processors together to achieve 1011-12 flop/second
- Aggregate available memories of 1011-12 bytes

*~100,000,000,000 flop/s*
*48 12km WRF CONUS in under 15 minutes*

## Parallel computing terms -- hardware

**Processor**:
- A device that reads and executes instructions in sequence to produce perform operations on data that it gets from a memory device producing results that are stored back onto the memory device
**Node**: One memory device connected to one or more processors.
- Multiple processors in a node are said to share-memory and this is “shared memory parallelism”
- They can work together because they can see each other’s memory
- The latency and bandwidth to memory affect performance
**Cluster**: Multiple nodes connected by a network
- The processors attached to the memory in one node can not see the memory for processors on another node
- For processors on different nodes to work together they must send messages between the nodes. This is  “distributed memory parallelism”
**Network: **
- Devices and wires for sending messages between nodes
- Bandwidth – a measure of the number of bytes that can be moved in a second
- Latency – the amount of time it takes before the first byte of a message arrives at its destination

**Notes:**
jjjj

## Parallel Computing Terms – System Software

- Process:
- A set of instructions to be executed on a processor
- Enough state information to allow process execution to stop on a processor and be picked up again later, possibly by another processor
- Processes may be lightweight or heavyweight
- Lightweight processes, e.g. shared-memory threads, store very little state; just enough to stop and then start the process
- Heavyweight processes, e.g. UNIX processes, store a lot more (basically the memory image of the job)

*“**The only thing one does directly with hardware is pay for it.”*

## Parallel Computing Terms – System Software

- Every job has at least one heavy-weight *process*.
- A job with more than one process is a distributed-memory parallel job
- Even on the same node, heavyweight processes do not share memory†
- Within a heavyweight process you may have some number of lightweight processes, called *threads.*
- Threads are shared-memory parallel; only threads in the same memory space can work together*. *
- A thread never exists by itself; it is always inside a heavy-weight process.
- Processes (heavy-weight) are the vehicles for distributed memory parallelism
- Threads are the vehicles for shared-memory parallelism

## Jobs, Processes, and Hardware

- MPI is used to start up and pass messages between multiple heavyweight processes
- The **mpirun **command** **controls the number of processes and how they are mapped onto nodes of the parallel machine
- Calls to MPI routines send and receive messages and control other interactions between processes
- http://www.mcs.anl.gov/mpi
- OpenMP is used to start up and control threads within each process
- Directives specify which parts of the program are multi-threaded
**OpenMP** environment variables determine the number of threads in each process
- http://www.openmp.org
- The number of processes (number of MPI processes times the number of threads in each process) usually corresponds to the number of processors

## Examples

- If the machine consists of 4 nodes, each with 4 processors, how many different ways can you run a job to use all 16?
- 4 MPI processes, each with 4 threads
**setenv OMP_NUM_THREADS 4**
**mpirun –np 4 wrf.exe**
- 8 MPI processes, each with 2 threads
**setenv OMP_NUM_THREADS 2**
**mpirun –np 8 wrf.exe**
- 16 MPI processes, each with 1 thread
**setenv OMP_NUM_THREADS 1**
**mpirun –np 16 wrf.exe**

## Examples (cont.)

- Note, since there are 4 nodes, we can never have fewer than 4 MPI processes because nodes do not share memory
- What happens on this same machine for the following?
**setenv OMP_NUM_THREADS 4**
**mpirun –np 32**

## Application:  WRF

- WRF uses *domain decomposition* to divide total amount of work over parallel processes
- Since the process model has two levels, the decomposition of the application over processes has two levels:
- The domain is first broken up into rectangular pieces that are assigned to heavy-weight processes. These pieces are called *patches*
- The *patches *may be further subdivided into smaller rectangular pieces that are called *tiles, *and these are assigned to *threads* within the process*.*

## Single version of code for efficient execution on:
Distributed-memory
Shared-memory 
Clusters of SMPs
Vector and microprocessors

**Model domains are decomposed for parallelism on two-levels**
***Patch:*** section of model domain  allocated to a distributed memory  node
***Tile:*** section of a patch allocated to a shared-memory processor within a node; this is also the scope of a model layer subroutine.
- Distributed memory parallelism is over patches; shared memory parallelism is over tiles within patches

**Single version of code for efficient execution on:**
- Distributed-memory
- Shared-memory
- Clusters of SMPs
- Vector and microprocessors

- Parallelism in WRF: Multi-level Decomposition

- Logical domain

- 1 Patch, divided into multiple tiles

- Inter-processor communication

## Distributed Memory Communications

**                      ****(module_diffusion.F )**
**SUBROUTINE horizontal_diffusion_s (tendency, rr, var, . . .**
**. . .**
**   ****DO j = jts,jte**
**   ****DO k = kts,ktf**
**   ****DO i = its,ite**
**      ****mrdx=msft(i,j)*rdx**
**      ****mrdy=msft(i,j)*rdy**
**      ****tendency(i,k,j)=tendency(i,k,j)-                          &**
**           ****(mrdx*0.5*((rr(****i+1****,k,j)+rr(i,k,j))*H1(****i+1****,k,j)-      &**
**                      ****(rr(****i-1****,k,j)+rr(i,k,j))*H1(i  ,k,j))+     &**
**            ****mrdy*0.5*((rr(i,k,****j+1****)+rr(i,k,j))*H2(i,k,****j+1****)-      &**
**                      ****(rr(i,k,****j-1****)+rr(i,k,j))*H2(i,k,j  ))-     &**
**            ****msft(i,j)*(H1avg(i,k+1,j)-H1avg(i,k,j)+             &**
**                       ****H2avg(i,k+1,j)-H2avg(i,k,j)              &**
**                                ****)/dzetaw(k)                     &**
**           ****)**
**   ****ENDDO**
**   ****ENDDO**
**   ****ENDDO**
** ****. . .**

- Example code fragment that requires communication between patches
- Note the tell-tale +1 and –1 expressions in indices for **rr**, **H1,** and **H2 ** arrays on right-hand side of assignment. These are *horizontal data dependencies *because the indexed operands may lie in the patch of a neighboring processor. That neighbor’s updates to that element of the array won’t be seen on this processor. We have to communicate.

## Halo updates

- Distributed Memory MPI Communications

- memory on one processor

- memory on neighboring processor

*****
**+ ***
*****

*****

## Halo updates

- Distributed Memory MPI Communications

- memory on one processor

- memory on neighboring processor

*****
**+ ***
*****

*****

- *

## Halo updates
Periodic boundary updates
Parallel transposes
Nesting scatters/gathers

- Periodic boundary updates
- Parallel transposes
- Nesting scatters/gathers

- Distributed Memory (MPI) Communications

## Review

- APPLICATION
- (WRF)

- HARDWARE
- (Processors, Memories, Wires)

- SYSTEM
- (UNIX, MPI, OpenMP)

| Domain | contains | Patches | contain | Tiles |
| --- | --- | --- | --- | --- |
| Job | contains | Processes | contain | Threads |
| Cluster | contains | Nodes | contain | Processors |

- Distributed
- Memory
- Parallel

- Shared
- Memory
- Parallel

## WRF Software Overview

- Architecture
- Directory structure
- Model Layer Interface
- Data Structures
- I/O
- Registry

## WRF Software Architecture

- Hierarchical software architecture
- Insulate scientists' code from parallelism and other architecture/implementation-specific details
- Well-defined interfaces between layers, and external packages for communications, I/O, and model coupling facilitates code reuse and exploiting of community infrastructure, e.g. ESMF.

- Registry

## WRF Software Architecture

- Driver Layer
- Allocates, stores, decomposes model domains, represented abstractly as single data objects
- Contains top-level time loop and algorithms for integration over nest hierarchy
- Contains the calls to I/O, nest forcing and feedback routines supplied by the Mediation Layer
- Provides top-level, non package-specific access to communications, I/O, etc.
- Provides some utilities, for example **module_wrf_error**, which is used for diagnostic prints and error stops

- Registry

## WRF Software Architecture

- Mediation Layer
- Provides to the Driver layer
- Solve solve routine, which takes a domain object and advances it one time step
- I/O routines that Driver calls when it is time to do some input or output operation on a domain
- Nest forcing and feedback routines
- The Mediation Layer and not the Driver knows the specifics of what needs to be done
- The sequence of calls to Model Layer routines for doing a time-step is known in Solve routine
- Responsible for dereferencing driver layer data objects so that individual fields can be passed to Model layer Subroutines
- Calls to message-passing are contained here as part of solve routine

- Registry

## WRF Software Architecture

- Model Layer
- Contains the information about the model itself, with machine architecture and implementation aspects abstracted out and moved into layers above
- Contains the actual WRF model routines that are written to perform some computation over an arbitrarily sized/shaped subdomain
- All state data objects are simple types, passed in through argument list
- Model Layer routines don’t know anything about communication or I/O; and they are designed to be executed safely on **one thread **– they _never_ contain a **PRINT**, **WRITE**, or **STOP** statement
- These are written to conform to the Model Layer Subroutine Interface (more later) which makes them “tile-callable”

- Registry

## WRF Software Architecture

- Registry: an “Active” data dictionary
- Tabular listing of model state and attributes
- Large sections of interface code generated automatically
- Scientists manipulate model state simply by modifying Registry, without further knowledge of code mechanics

- Registry

**Notes:**
Dimensions, I/O and model coupling, communications, nesting, etc.

Interfaces between layers, communications

Approximately 1⁄4 of wrf is autogen’d in this way

## Call Structure superimposed on Architecture

## Call Structure superimposed on Architecture

## Call Structure superimposed on Architecture

**solve_em (dyn_em/solve_em.F)**

## Call Structure superimposed on Architecture

**solve_em (dyn_em/solve_em.F)**

## WRF Model Layer Interface

- Mediation layer / Model Layer Interface
- All state arrays passed through argument list as simple (not derived) data types
- Domain, memory, and run dimensions passed unambiguously in three physical dimensions
- Model layer routines are called from mediation layer in loops over tiles, which are multi-threaded
- Restrictions on model layer subroutines
- No I/O, communication, no stops or aborts (use wrf_error_fatal in frame/module_wrf_error.F)
- No common/module storage of decomposed data (exception allowed for set-once/read-only tables)
- Spatial scope of a Model Layer call is one “tile”
- Temporal scope of a call is limited by coherency

## WRF Model Layer Interface

- Mediation layer / Model Layer Interface
- Model layer routines are called from mediation layer in loops over tiles, which are multi-threaded
- All state arrays passed through argument list as simple data types
- Domain, memory, and run dimensions passed unambiguously in three physical dimensions
- Restrictions on model layer subroutines
- No I/O, communication, no stops or aborts (use wrf_error_fatal in frame/module_wrf_error.F)
- No common/module storage of decomposed data (exception allowed for set-once/read-only tables)
- Spatial scope of a Model Layer call is one “tile”
- Temporal scope of a call is limited by coherency

- SUBROUTINE solve_xxx (
-     . . .
- !$OMP DO PARALLEL
-    DO ij = 1, numtiles
-       its = i_start(ij) ; ite = i_end(ij)
-       jts = j_start(ij) ; jte = j_end(ij)
-       CALL model_subroutine( arg1, arg2, . . .
-            ids , ide , jds , jde , kds , kde ,
-            ims , ime , jms , jme , kms , kme ,
-            its , ite , jts , jte , kts , kte )
-    END DO
-     . . .
-  END SUBROUTINE

*           *_*template for model layer subroutine       *_
-  SUBROUTINE model_subroutine ( &
-    arg1, arg2, arg3, ... , argn,   &
-    ids, ide, jds, jde, kds, kde, &  ! Domain dims
-    ims, ime, jms, jme, kms, kme, &  ! Memory dims
-    its, ite, jts, jte, kts, kte  )  ! Tile dims
-  IMPLICIT NONE
-  ! Define Arguments (S and I1) data
-  REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: arg1, . . .
-  REAL, DIMENSION (ims:ime,jms:jme)         :: arg7, . . .
-   . . .
-  ! Define Local Data (I2)
-  REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: loc1, . . .
-   . . .
-  ! Executable code; loops run over tile
-  ! dimensions
-  DO j = jts, jte
-    DO k = kts, kte
-      DO i = MAX(its,ids), MIN(ite,ide)
-        loc(i,k,j) = arg1(i,k,j) + ...
-      END DO
-    END DO
-  END DO

## Domain dimensions

*           *_*template for model layer subroutine       *_
-  SUBROUTINE model ( &
-    arg1, arg2, arg3, ... , argn,   &
-    ids, ide, jds, jde, kds, kde, &  ! Domain dims
-    ims, ime, jms, jme, kms, kme, &  ! Memory dims
-    its, ite, jts, jte, kts, kte  )  ! Tile dims
-  IMPLICIT NONE
-  ! Define Arguments (S and I1) data
-  REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: arg1, . . .
-  REAL, DIMENSION (ims:ime,jms:jme)         :: arg7, . . .
-   . . .
-  ! Define Local Data (I2)
-  REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: loc1, . . .
-   . . .
-  ! Executable code; loops run over tile
-  ! dimensions
-  DO j = jts, jte
-    DO k = kts, kte
-      DO i = MAX(its,ids), MIN(ite,ide)
-        loc(i,k,j) = arg1(i,k,j) + ...
-      END DO
-    END DO
-  END DO

- Size of logical domain
- Used for bdy tests, etc.

## Domain dimensions

*           *_*template for model layer subroutine       *_
-  SUBROUTINE model ( &
-    arg1, arg2, arg3, ... , argn,   &
-    ids, ide, jds, jde, kds, kde, &  ! Domain dims
-    ims, ime, jms, jme, kms, kme, &  ! Memory dims
-    its, ite, jts, jte, kts, kte  )  ! Tile dims
-  IMPLICIT NONE
-  ! Define Arguments (S and I1) data
-  REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: arg1, . . .
-  REAL, DIMENSION (ims:ime,jms:jme)         :: arg7, . . .
-   . . .
-  ! Define Local Data (I2)
-  REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: loc1, . . .
-   . . .
-  ! Executable code; loops run over tile
-  ! dimensions
-  DO j = jts, jte
-    DO k = kts, kte
-      DO i = MAX(its,ids), MIN(ite,ide)
-        loc(i,k,j) = arg1(i,k,j) + ...
-      END DO
-    END DO
-  END DO

- Size of logical domain
- Used for bdy tests, etc.
- Memory dimensions
- Used to dimension dummy arguments
- Do not use for local arrays

## Domain dimensions

*           *_*template for model layer subroutine       *_
-  SUBROUTINE model ( &
-    arg1, arg2, arg3, ... , argn,   &
-    ids, ide, jds, jde, kds, kde, &  ! Domain dims
-    ims, ime, jms, jme, kms, kme, &  ! Memory dims
-    its, ite, jts, jte, kts, kte  )  ! Tile dims
-  IMPLICIT NONE
-  ! Define Arguments (S and I1) data
-  REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: arg1, . . .
-  REAL, DIMENSION (ims:ime,jms:jme)         :: arg7, . . .
-   . . .
-  ! Define Local Data (I2)
-  REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: loc1, . . .
-   . . .
-  ! Executable code; loops run over tile
-  ! dimensions
-  DO j = jts, jte
-    DO k = kts, kte
-      DO i = MAX(its,ids), MIN(ite,ide)
-        loc(i,k,j) = arg1(i,k,j) + ...
-      END DO
-    END DO
-  END DO

- Size of logical domain
- Used for bdy tests, etc.
- Memory dimensions
- Used to dimension dummy arguments
- Do not use for local arrays
- Tile dimensions
- Local loop ranges
- Local array dimensions

## Domain dimensions

*           *_*template for model layer subroutine       *_
-  SUBROUTINE model ( &
-    arg1, arg2, arg3, ... , argn,   &
-    ids, ide, jds, jde, kds, kde, &  ! Domain dims
-    ims, ime, jms, jme, kms, kme, &  ! Memory dims
-    its, ite, jts, jte, kts, kte  )  ! Tile dims
-  IMPLICIT NONE
-  ! Define Arguments (S and I1) data
-  REAL, DIMENSION (ims:ime,kms:kme,jms:jme) :: arg1, . . .
-  REAL, DIMENSION (ims:ime,jms:jme)         :: arg7, . . .
-   . . .
-  ! Define Local Data (I2)
-  REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: loc1, . . .
-   . . .
-  ! Executable code; loops run over tile
-  ! dimensions
-  DO j = jts, jte
-    DO k = kts, kte
-      DO i = MAX(its,ids), MIN(ite,ide)
-        loc(i,k,j) = arg1(i,k,j) + ...
-      END DO
-    END DO
-  END DO

- Size of logical domain
- Used for bdy tests, etc.
- Memory dimensions
- Used to dimension dummy arguments
- Do not use for local arrays
- Tile dimensions
- Local loop ranges
- Local array dimensions

- Patch dimensions
- Start and end indices of local distributed memory subdomain
- Available from mediation layer (solve) and driver layer; not usually needed or used at model layer

## Data Structures

## Data structures

- How is model state represented and stored?
- What do data objects look like in WRF?
- How do you manipulate these?

## Grid Representation in Arrays 

- Increasing indices in WRF arrays run
- West to East   (X, or I-dimension)
- South to North (Y, or J-dimension)
- Bottom to Top (Z, or K-dimension)
- Storage order in WRF is IKJ but this is a WRF Model convention, not a restriction of the WRF Software Framework
- The extent of the logical or *domain* dimensions is always the "staggered" grid dimension. That is, from the point of view of a non-staggered dimension, there is always an extra cell on the end of the domain dimension.

## Grid Indices Mapped onto Array Indices (C-grid example)

|  |  |  |  |
| --- | --- | --- | --- |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |

- jds = 1

- jde = 4

- ids = 1

- ide = 4

- Computation over
- mass points runs
- only ids..ide-1
- and  jds..jde-1
- Likewise, vertical
- computation over unstaggered fields
- run kds..kde-1

## What does data look like in WRF?

- At the highest level (Driver Layer)
- Represented as dynamically allocated fields in “derived data types” (DDTs) that contain all state for a “domain”
- Linked together in a tree to represent nest hierarchy; root pointer is head_grid, defined in frame/module_domain.F
- Recursive depth-first traversal for nesting

## What does data look like in WRF?

- At the highest level (Driver Layer)
- Represented as dynamically allocated fields in “derived data types” (DDTs) that contain all state for a “domain”
- Linked together in a tree to represent nest hierarchy; root pointer is head_grid, defined in frame/module_domain.F
- Driver passes the DDT for a domain to the solve routine, whose job it is to advance the state of that domain by one time step
- In the solver (Mediation Layer)
- Solve routine takes DDT argument named “grid” and *dereferences* fields as necessary over the course of computing one time step

- SUBROUTINE solve_em(  grid  ... )
-    TYPE(domain) :: grid
-   ...
-    CALL model_layer ( grid%u , grid%v , ... )
-   ...
- END SUBROUTINE

## What does data look like in WRF?

- At the highest level (Driver Layer)
- Represented as dynamically allocated fields in “derived data types” (DDTs) that contain all state for a “domain”
- Linked together in a tree to represent nest hierarchy; root pointer is head_grid, defined in frame/module_domain.F
- Driver passes the DDT for a domain to the solve routine, whose job it is to advance the state of that domain by one time step
- In the solver (Mediation Layer)
- Solve routine takes DDT argument named “grid” and *dereferences* fields as necessary over the course of computing one time step
- At the lowest level (Model Layer)
- All state data is simple Fortran data types and arrays, passed through the argument list (see Model Layer Interface, above)

## How is data manipulated in WRF?

- Things to do with a state field
- Define it (or remove it)
- Specify I/O on it
- Specify communication on it
- These get done using the Registry (coming up)

## I/O

## WRF I/O

- Streams: pathways into and out of model
- History + 11 auxiliary output streams  (10 &11 reserved for nudging)
- Input + 11 auxiliary input streams (10 & 11 reserved for nudging)
- Restart and boundary
- Attributes of streams
- Variable set
- The set of WRF state variables that comprise one read or write on a stream
- Defined for a stream at compile time in Registry
- Format
- The format of the data outside the program (e.g. NetCDF)
- Specified for a stream at run time in the namelist
- Additional namelist-controlled attributes of streams
- Dataset name
- Time interval between I/O operations on stream
- Starting, ending times for I/O (specified as intervals from start of run)

## The Registry

## WRF Registry

- "Active data-dictionary” for managing WRF data structures
- Database describing attributes of model state, intermediate, and configuration data
- Dimensionality, number of time levels, staggering
- Association with physics
- I/O classification (history, initial, restart, boundary)
- Communication points and patterns
- Configuration lists (e.g. namelists)
- Program for auto-generating sections of WRF from database:
- _570_ Registry entries  _30-thousand_ lines of automatically generated WRF code
- Allocation statements for state data, I1 data
- Argument lists for driver layer/mediation layer interfaces
- Interprocessor communications: Halo and periodic boundary updates, transposes
- Code for defining and managing run-time configuration information
- Code for forcing, feedback and interpolation of nest data
- Why?
- Automates time consuming, repetitive, error-prone programming
- Insulates programmers and code from package dependencies
- Allow rapid development
- Documents the data
- Reference: **Description of WRF Registry, **http://www.mmm.ucar.edu/wrf/software_v2

## Registry Mechanics

- WRF source
- */*.F

- CPP
- ____________
- Fortran90

## Registry Data Base

- Currently implemented as a text file: Registry/Registry.EM
- Types of entry:
*Dimspec* – Describes dimensions that are used to define arrays in the model
*State* – Describes state variables and arrays in the domain structure
*I1* – Describes local variables and arrays in solve
*Typedef* – Describes derived types that are subtypes of the domain structure
*Rconfig* – Describes a configuration (e.g. namelist) variable or array
*Package* – Describes attributes of a package (e.g. physics)
*Halo* – Describes halo update interprocessor communications
*Period* – Describes communications for periodic boundary updates
*Xpose* – Describes communications for parallel matrix transposes

## Registry State Entry: ordinary State 

- Elements
*Entry*: The keyword “state”
*Type*: The type of the state variable or array (real, double, integer, logical, character, or derived)
*Sym*: The symbolic name of the variable or array
*Dims*: A string denoting the dimensionality of the array or a hyphen (-)
*Use*: A string denoting association  4D scalar array (only significant if *dims* contains **f**)
*NumTLev*: An integer indicating the number of time levels (for arrays) or hypen (for variables)
*Stagger*: String indicating staggered dimensions of variable  (X, Y, Z, or hyphen)
*IO*: String indicating whether and how the variable is subject to I/O and Nesting
*DName*: Metadata name for the variable
*Units*: Metadata units of the variable
*Descrip*: Metadata description of the variable
- Example

**#      Type Sym  Dims   Use     Tlev Stag IO       Dname      Descrip**
- state  real  u   ikj**b**   dyn_em   2   X    irhusdf   "U"    "X WIND COMPONENT“

## Registry State Entry: ordinary State 

- This single entry results in 130 lines automatically added to 43 different locations of the WRF code:
- Declaration and dynamic allocation of arrays in TYPE(domain)
- Two 3D state arrays corresponding to the 2 time levels of U
- u_1 ( ims:ime , kms:kme , jms:jme )
- u_2 ( ims:ime , kms:kme , jms:jme )
- Two sets of LBC arrays for boundary and boundary tendencies
- u_bxs( jms:jme, kms:kme, spec_bdy_width )   ! west boundary
- u_bxe( jms:jme, kms:kme, spec_bdy_width )   ! east boundary
- u_bys( ims:ime, kms:kme, spec_bdy_width )   ! south boundary
- u_bye( ims:ime, kms:kme, spec_bdy_width )   ! north boundary
- Nesting code to interpolate, force, feedback, and smooth **u**
- Addition of **u** to the input, restart, history, and LBC  I/O streams

**#      Type Sym  Dims  Use  Tlev Stag    IO    Dname      Descrip**
- state  real  u   ikj**b**   -     2   X    irh**usdf**   "U"    "X WIND COMPONENT“

## State Entry: Defining a variable-set for an I/O stream

- Fields are added to a variable-set on an I/O stream in the Registry

_*IO*_ is a string that specifies if the variable is to be subject to initial, restart, history, or boundary I/O.  The string may consist of 'h' (subject to history I/O), 'i' (initial dataset), 'r' (restart dataset), or 'b' (lateral boundary dataset).  The 'h', 'r', and 'i' specifiers may appear in any order or combination.
- The ‘h’ and ‘i’ specifiers may be followed by an optional integer string consisting of ‘0’, ‘1’, ‘2’, ‘3’, ‘4’, and/or ‘5’.  Zero denotes that the variable is part of the principal input or history I/O stream. The characters ‘1’ through ‘5’ denote one of five auxiliary input or history I/O streams.

- #      Type Sym  Dims   Use     Tlev Stag IO    Dname      Descrip
- state  real  u   ikjb   dyn_em   2   X    irh   "U"  "X WIND COMPONENT“

## State Entry: Defining Variable-set for an I/O stream

- irh -- The state variable will be included in the input, restart, and history I/O streams
- irh13  -- The state variable has been added to the first and third auxiliary history output streams; it has been removed from the principal history output stream, because zero is not among the integers in the integer string that follows the character 'h'
- rh01 -- The state variable has been added to the first auxiliary history output stream; it is also retained in the principal history output
- i205hr  -- Now the state variable is included in the principal input stream as well as auxiliary inputs 2 and 5.  Note that the order of the integers is unimportant. The variable is also in the principal history output stream
- ir12h  -- No effect; there is only 1 restart data stream and ru added to it.

## Rconfig entry

- This defines namelist entries
- Elements
*Entry*: the keyword “rconfig”
*Type*: the type of the namelist variable (integer, real, logical, string )
*Sym*: the name of the namelist variable or array
*How* *set*: indicates how the variable is set: e.g. namelist or derived, and if namelist, which block of the namelist it is set in
*Nentries*: specifies the dimensionality of the namelist variable or array. If 1 (one) it is a variable and applies to all domains; otherwise specify max_domains (which is an integer parameter defined in module_driver_constants.F).
*Default*: the default value of the variable to be used if none is specified in the namelist; hyphen (-) for no default
- Example

**#         Type       Sym              How set          Nentries  Default**
- rconfig   integer spec_bdy_width  namelist,bdy_control    1          1

## Rconfig entry

- Result of this Registry Entry:
- Define an namelist variable “spec_bdy_width” in the bdy_control section of namelist.input
- Type integer (others: real, logical, character)
- If this is first entry in that section, define “bdy_control” as a new section in the namelist.input file
- Specifies that bdy_control applies to all domains in the run
- if Nentries is “max_domains” then the entry in the namelist.input file is a comma-separate list, each element of which applies to a separate domain
- Specify a default value of “1” if nothing is specified in the namelist.input file
- In the case of a multi-process run, generate code to read in the bdy_control section of the namelist.input file on one process and broadcast the value to all other processes

**#         Type       Sym              How set          Nentries  Default**
- rconfig   integer spec_bdy_width  namelist,bdy_control    1          1

-  **--- File: namelist.input ---**
- &bdy_control
-  **spec_bdy_width       = 5**,
-  spec_zone            = 1,
-  relax_zone           = 4,
-     . . .
-  /

## Examples: working with WRF software

-    Add a new physics package with time varying input source to the model

## Example: Input periodic SSTs

- Problem: adapt WRF to input a time-varying lower boundary condition, e.g. SSTs, from an input file for a new surface scheme
- Given: Input file in WRF I/O format containing 12-hourly SST’s
- Modify WRF model to read these into a new state array and make available to WRF surface physics

## Example: Input periodic SSTs

- Steps
- Add a new state variable and definition of a new surface layer package that will use the variable to the Registry
- Add to variable stream for an unused Auxiliary Input stream
- Adapt physics interface to pass new state variable to physics
- Setup namelist to input the file at desired interval

## Example: Input periodic SSTs

- Add a new state variable to Registry/Registry.EM and put it in the variable set for input on AuxInput #3
- Also added to History and Restart
- Result:
- 2-D variable named **nsst** defined and available in solve_em
- Dimensions: ims:ime, jms:jme
- Input and output on the AuxInput #3 stream will include the variable under the name NEW_SST

**#     type  symbol dims use  tl stag  io     dname     description       units**
**state real  nsst   ij   misc  1  -    i3rh  "NEW_SST" "Time Varying SST" "K“**

## Example: Input periodic SSTs

- Pass new state variable to surface physics

**        ****--- File: dyn_em/solve_em.F ---**
**CALL surface_driver(                                                    &**
**           ****. . .**
**   ****&   ,**_**NSST=grid%nsst**_**                                                  & ! new**
**   ****&   ,CAPG=grid%capg, EMISS=grid%emiss, HOL=grid%hol,MOL=grid%mol     &**
**   ****&   ,RAINBL=grid%rainbl                                              &**
**   ****&   ,RAINNCV=grid%rainncv,REGIME=grid%regime,T2=grid%t2,THC=grid%thc & **
**           ****. . .**
**   ****&                                                              )**

## Example: Input periodic SSTs

- Add new variable **nsst **to Physics Driver in Mediation Layer
- By making this an “Optional” argument, we preserve the driver’s compatibility with other cores and with versions of WRF where this variable hasn’t been added.

**        ****--- File: phys/module_surface_driver.F ---**
**   ****SUBROUTINE surface_driver(                                         &**
**          ****. . .**
**     ****&          ,**_**nsst**_**                                                 &**
**     ****&          ,capg,emiss,hol,mol                                   &**
**     ****&          ,rainncv,rainbl,regime,t2,thc                         &**
**                                                                      ****&**
**                                                                      ****))**
**         ****. . .**
**   ****REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT)::   nsst**

## Example: Input periodic SSTs

- Add call to Model-Layer subroutine for new physics package to Surface Driver
- Note the PRESENT test to make sure new optional variable **nsst** is available

**        ****--- File: phys/module_surface_driver ---**
**     ****sfclay_select: SELECT CASE(sf_sfclay_physics)**
**       ****CASE (SFCLAYSCHEME)**
**          ****. . .**
**       ****CASE (NEWSFCSCHEME)  ! <- This is defined by the Registry “package” entry**
**         ****IF (**_**PRESENT(nsst)**_**)  THEN**
**            ****CALL NEWSFCCHEME(                                        &**
**                ****nsst,                                                &**
**                ****ids,ide, jds,jde, kds,kde,                           &**
**                ****ims,ime, jms,jme, kms,kme,                           &**
**                ****i_start(ij),i_end(ij), j_start(ij),j_end(ij), kts,kte    )**
**         ****ELSE**
**           ****CALL wrf_error_fatal('Missing argument for NEWSCHEME in surface driver')**
**         ****ENDIF**
**          ****. . .**
**     ****END SELECT sfclay_select**

## Example: Input periodic SSTs

- Add definition for new physics package NEWSCHEME as setting 4 for namelist variable sf_sfclay_physics
- This creates a defined constant NEWSFCSCHEME and represents selection of the new scheme when the namelist variable sf_sfclay_physics is set to ‘4’ in the namelist.input file
- Clean –a and recompile so code and Registry changes take effect

**rconfig   integer  sf_sfclay_physics   namelist,physics   max_domains    0**
**package   sfclayscheme   sf_sfclay_physics==1        -             -**
**package   myjsfcscheme   sf_sfclay_physics==2        -             -**
**package   gfssfcscheme   sf_sfclay_physics==3        -             -**
_**package   newsfcscheme   sf_sfclay_physics==4**_**        -             -**

## Example: Input periodic SSTs

- Setup namelist to input SSTs from the file at desired interval
- Run code with sst_input file in run-directory

**     ****--- File: namelist.input ---**
**&time_control**
**   ****. . .**
** ****auxinput3_inname      = "sst_input"**
** ****auxinput3_interval_mo = 0**
** ****auxinput3_interval_d  = 0**
** ****auxinput3_interval_h  = 12**
** ****auxinput3_interval_m  = 0**
** ****auxinput3_interval_s  = 0**
**   ****. . .**
**/**
**   ****. . .**
**&physics**
** ****sf_sfclay_physics  = 4, 4, 4**
**   ****. . .**
**/**

## Example: Input periodic SSTs

- A few notes...
- The read times and the time-stamps in the input file must match exactly
- We haven’t done anything about what happens if the file runs out of time periods (the last time period read will be used over and over again, though you’ll see some error messages in the output if you set debug_level to be 1 or greater in namelist.input)
- We haven’t said anything about what generates sst_input

## Example: Working with WRF Software

- Computing and outputting a Diagnostic

## Example: Compute a Diagnostic

- Problem: Output global average and global maximum and lat/lon location of maximum for 10 meter wind speed in WRF
- Steps:
- Modify solve to compute wind-speed and then compute the local sum and maxima at the end of each time step
- Use reduction operations built-in to WRF software to compute the global qualitities
- Output these on one process (process zero, the “monitor” process)

## Example: Compute a Diagnostic

- Compute local sum and local max and the local indices of the local maximum

-    **--- File: dyn_em/solve_em.F  (near the end) ---**
- ! Compute local maximum and sum of 10m wind-speed
-    sum_ws = 0.
-    max_ws = 0.
-    DO j = jps, jpe
-      DO i = ips, ipe
-        wind_vel = sqrt( grid%u10(i,j)*grid%u10(i,j) + grid%v10(i,j)*grid%v10(i,j) )
-        IF ( wind_vel .GT. max_ws ) THEN
-           **max_ws** = wind_vel
-           **idex** = i
-           **jdex** = j
-        ENDIF
-        **sum_ws** = sum_ws + wind_vel
-      ENDDO
-    ENDDO

## Example: Compute a Diagnostic

- Compute global sum,  global max, and indices of the global max

- ! Compute global sum
-    **sum_ws** = wrf_dm_sum_real ( sum_ws )
- ! Compute global maximum and associated i,j point
-    CALL wrf_dm_maxval_real ( **max_ws**, **idex, jdex** )

## Example: Compute a Diagnostic

- On the process that contains the maximum value, obtain the latitude and longitude of that point; on other processes set to an artificially low value.
- The use parallel reduction to store that result on every process

-    IF ( ips .LE. idex .AND. idex .LE. ipe .AND.  &
-         jps .LE. jdex .AND. jdex .LE. jpe ) THEN
-       glat = grid%xlat(idex,jdex)
-       glon = grid%xlong(idex,jdex)
-    ELSE
-       glat = -99999.
-       glon = -99999.
-    ENDIF
- ! Compute global maximum to find glat and glon
-    **glat** = wrf_dm_max_real ( glat )
-    **glon** = wrf_dm_max_real ( glon )

## Example: Compute a Diagnostic

- Output the value on process zero, the “monitor”
- Output from process zero of a 4 process run

- ! Print out the result on the monitor process
-    IF ( wrf_dm_on_monitor() ) THEN
-       WRITE(outstring,*)'Avg. ',sum_ws/((ide-ids*1)*(jde-jds+1))
-       CALL wrf_message ( TRIM(outstring) )
-       WRITE(outstring,*)'Max. ',max_ws,' Lat. ',glat,' Lon. ',glon
-       CALL wrf_message ( TRIM(outstring) )
-    ENDIF

-      **--- Output file: rsl.out.0000 ---**
-     . . .
-   Avg.    5.159380
-   Max.    15.09370     Lat.    37.25022     Lon.   -67.44571
- Timing for main: time 2000-01-24_12:03:00 on domain   1:    8.96500 elapsed seconds.
-   Avg.    5.166167
-   Max.    14.97418     Lat.    37.25022     Lon.   -67.44571
- Timing for main: time 2000-01-24_12:06:00 on domain   1:    4.89460 elapsed seconds.
-   Avg.    5.205693
-   Max.    14.92687     Lat.    37.25022     Lon.   -67.44571
- Timing for main: time 2000-01-24_12:09:00 on domain   1:    4.83500 elapsed seconds.
-     . . .

## Example: Compute a Diagnostic (complete)

-    REAL           :: sum_ws, max_ws, glat, glon, wind_vel
-    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
-    CHARACTER*256  :: outstring
-     . . .
**! Compute local maximum and sum of 10m wind-speed**
-    sum_ws = 0.
-    max_ws = 0.
-    DO j = jps, jpe
-      DO i = ips, ipe
-        wind_vel = sqrt( u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))
-        IF ( wind_vel .GT. max_ws ) THEN
-           max_ws = wind_vel
-           idex = i
-           jdex = j
-        ENDIF
-        sum_ws = sum_ws + wind_vel
-      ENDDO
-    ENDDO
**! Compute global sum**
-    sum_ws = wrf_dm_sum_real ( sum_ws )
**! Compute global maximum and associated i,j point**
-    CALL wrf_dm_maxval_real ( max_ws, idex, jdex )
- ! **Deterimine if i,j point of maximum is on this process**
**! and if so, set the lat and lon of that point, otherwise**
**!** set to an absolute minimum
-    IF ( ips .LE. idex .AND. idex .LE. ipe .AND.  &
-         jps .LE. jdex .AND. jdex .LE. jpe ) THEN
-       glat = xlat(idex,jdex)
-       glon = xlat(idex,jdex)
-    ELSE
-       glat = -99999.
-       glon = -99999.
-    ENDIF
**! Compute global maximum to find glat and glon**
-    glat = wrf_dm_max_real ( glat )
-    glon = wrf_dm_max_real ( glon )
**! Print out the result on the monitor process**
-    IF ( wrf_dm_on_monitor() ) THEN
-       WRITE(outstring,*)'Avg. ',sum_ws/((ide-ids*1)*(jde-jds+1))
-       CALL wrf_message ( TRIM(outstring) )
-       WRITE(outstring,*)'Max. ',max_ws,' Lat. ',glat,' Lon. ',glon
-       CALL wrf_message ( TRIM(outstring) )
-    ENDIF

## Summary