Day 2 Adding a New Scheme to the CCPP

Day 2 Adding a New Scheme to the CCPP carson Mon, 03/04/2019 - 15:14

Adding a new scheme to the CCPP

In this exercise, we will implement a new scheme in CCPP that calculates the Modified Red Flag Threat Index (RFTIM) following https://www.weather.gov/media/ama/rftimod_primer.pdf

We will also be adding the RFTIM to the diagnostic output, run a short forecast and display the results.

The exercise is designed to touch all the relevant bits and pieces for writing a new scheme, adding interstitial variables and a diagnostic variable, and outputting the diagnostic variable to disk.

Prepare for the GitHub workflow

Prepare for the GitHub workflow carson Mon, 03/04/2019 - 14:51

Set up your development forks

In order to submit a pull request (PR) for the code changes at a later stage, it is necessary to create a GitHub user fork of the authoritative repositories for which changes are made. For this exercise, changes are made to the umbrella repository NEMSfv3gfs and the two submodules ccpp-physics and FV3. Creating a fork is a one-off task.

  1. Navigate to http://github.com/NCAR/NEMSfv3gfs, locate and click on the the "fork" button towards the top right of the page.
  2. Navigate to http://github.com/NCAR/FV3 locate and click on the the "fork" button 
  3. Navigate to http://github.com/NCAR/ccpp-physics locate and click on the the "fork" button 

If you already forked this repository in the past, you will see a message "You've already forked NEMSfv3gfs" and a link that you can click on to get to your fork (http://github.com/YOUR_GITHUB_USERNAME/NEMSfv3gfs). Otherwise, follow the instructions to create your user fork. Repeat for ccpp-physics and FV3.

Set up the directory structure for this hands-on session and checkout the source code.  

   mkdir ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training
   cd ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training
   git clone --branch=emc_training_march_2019 --recursive https://github.com/NCAR/NEMSfv3gfs develop

Note. Potential errors while checking out the submodule NCEPLIBS-pyprodutil (e.g., "gerrit: Name or service not known") can be ignored.

At this point, you should have a directory tree named "develop" which contains sub-directories for FV3, FMS, NEMS, ccpp/physics, ccpp/framework.  Each of these contains the tag "emc_training_march_2019", as a starting point for these exercises.

Configure remotes and create branches

Next, configure the remote destinations (remotes) so that your development is pushed to your fork and not (inadvertently) to the GMTB repository at NCAR.  Also, create branches in the three repositories where changes will be made:

   cd develop
   # Submodule ccpp-physics
   cd ccpp/physics
   git remote rename origin upstream
   git remote add origin https://github.com/YOUR_GITHUB_USERNAME/ccpp-physics
   git remote -v show
   git remote update
   git checkout -b my_rftim_implementation
   cd ../../
   # Submodule FV3
   cd FV3
   git remote rename origin upstream
   git remote add origin https://github.com/YOUR_GITHUB_USERNAME/FV3
   git remote -v show
   git remote update
   git checkout -b my_rftim_implementation
   cd ../
   # Umbrella repository NEMSfv3gfs
   git remote rename origin upstream
   git remote add origin https://github.com/YOUR_GITHUB_USERNAME/NEMSfv3gfs
   git remote -v show
   git remote update
   git checkout -b my_rftim_implementation
   

Preparing to Calculate the Modified Red Flag Threat Index

Preparing to Calculate the Modified Red Flag Threat Index carson Mon, 03/04/2019 - 14:53

The calculation of the RFTIM is based on the 2-m relative humidity in percent (rh2m) and on the 6-m wind speed in miles per hour (ws6m). It uses a geographically-dependent climatology of thresholds for the two quantities to return an integer value between 0 and 5 for each of the two variables. In our simplified version, we use the climatological values derived for KAMA (Rick Husband, Amarillo International Airport, Amarillo, TX), independent of the geographical location.  This calculation could be made more sophisticated by using a spatially-varying climatology and by considering additional variables, such as snow cover.

Add the new variables 

This calculation is done best in the existing scheme sfc_diag (ccpp/physics/physics/sfc_diag.f), because all the necessary quantities to derive relative humidity and wind speed are computed in the scheme. Before making changes to sfc_diag.f, note that the required input variables, rh2m and ws6m, are not provided by the model. The first task therefore is to add these two variables to FV3/gfsphysics/GFS_layer/GFS_typedefs.F90.

First let’s work with rh2m and ws6m Both rh2m and ws6m are pure interstitial variables, i.e. they do not need to be persistent in memory from one time step to another. It is thus recommended to add them to the GFS_interstitial data type.

Variable rh2m:
   standard_name : relative_humidity_at_2m
   units         : percent
   rank          : 1
   type          : real
   kind          : kind_phys

Variable ws6m:
   standard_name : wind_speed_at_6m_in_miles_per_hour
   units         : mi h-1
   rank          : 1
   type          : real
   kind          : kind_phys

FYI: When variables have non-standard units, such as miles per hour, the variables standard name are decorated with information about the units, This is the case with variable ws6m above.

Each variable that is added to FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 requires a minimum of 4 additions to the file:

  • a metadata table entry
  • a variable declaration inside one of the derived type definitions
  • an allocation statement within a *_create method for the derived type
  • an initialization statement.

For expediency, you can duplicate the entries pertaining to the variable zt1d for rh2m and ws6m, making the necessary adjustments to use the metadata above, that is, use the proper standard_name and units for each variable added. You will note that zt1d appears a 5th time, which is a print statement.

FYI: For some variables added to the GFS_interstitial_type or GFS_diag_type types, a few more additions are necessary for things like resetting their values in a *_reset method and adding a line to a *_print method.

Hint 1: Refer to presentation day2_5_Host-Side_Coding, section “Adding a variable to the host+CCPP”, “Case 1: the new variable is a member of a derived data type that already exists on the host model side”.

Hint 2: Using the standard_name and units exactly as listed above is critical for this development to work, These are the two aspects of metadata used in communication between the CCPP and the host model. Regarding the longname, you can use any description of your preference.

Hint 3: Use optional=F for all variable declarations in GFS_typedefs.F90

Next let’s work with rftim.  Variable rftim also needs to be added to  FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 in 4 locations:  metadata table, variable declaration, variable allocation, and variable initialization.

Variable rftim:
   standard_name : modified_red_flag_threat_index
   units         : index
   rank          : 1
   type          : real
   kind          : kind_phys

Since variable rftim will be output by the model, it needs to be added to the GFS_diag_type data type. You need to modify the metadata table where the various GFS_diag_type variables are listed and add a line for rftim. For example, you could add a line at the end of the list, right after the line for atmosphere_momentum_diffusivity_for_mynnpbl. Note: if adding to the end of the metadata table list, be sure to leave a blank line starting with !! at the bottom of the table -- this is needed for parsing by the CCPP framework.

Next, you need to declare variable rftim. Add the following to the list of variables in section “Output - only in physics”. You can place this after the declaration for “MP quantities for 3D

   !--- Modified Read Flag Thread Index
    real (kind=kind_phys), pointer :: rftim(:)       => null() !<

You need to allocate the variable rftim in subroutine diag_create and then initialize it. You can enter the text below right before the line call Diag%rad_zero  (Model)

   ! Modified Red Flag Threat Index calculation
   allocate (Diag%rftim (IM))
   Diag%rftim = clear_val

 

Calculate the New Variables

 

To calculate these variables in sfc_diag, you can use the following equations:

   !> - For Modified Red Flag Threat Index calculation,
   !! relative humidity 2m above ground in percent
   rh2m(i) = 100.0*q2m(i)/qss

   !> - Wind speed 6m above ground in miles per hour,
   !! conversion factor from 10m to 6m is 0.886
   !! (Lindley et al., 2011), and m/s to mi/h is 2.23694
   ws6m(i) = sqrt(u10m(i)**2 + v10m(i)**2) * 0.886 * 2.23694

Edit file ccpp/physics/physics/sfc_diag.f. You will need to do the following modifications:

  • Add the new variables (rh2m and ws6m) to the metadata table.
  • Add the new variables as arguments to the subroutine.
  • Declare the variables as real(kind=kind_phys), dimension(im), intent(out)
  • Compute the variables using the equations provided above
  • Print out the variable (optional)

To make these changes, you can follow what is done to variable u10m, since it also goes through the 5 steps above.

Hint 1: Don’t forget to add the new variables to the metadata table.  You can copy/paste the entries from GFS_typdefs.F90 and remember to change the intent to OUT.

Hint 2: Don’t forget to pass your variables in and out of the subroutine. Arguments in the subroutine call must be in the same order as in the table of metadata.  Also don’t forgot to declare them appropriately in the declaration section of the subroutine using intent(out) and dimension (im).

Calculate the Index and Add it to the Output

Calculate the Index and Add it to the Output carson Mon, 03/04/2019 - 14:57

Add diagnostic variables for output

In order to add RFTIM to the output files, we need to add a new diagnostic variable to NEMSfv3gfs that stores the RFTIM. Because the NEMSfv3gfs diagnostic output doesn't work with integers, we need to declare rftim as a real array constituent of GFS_diag_type:

Variable rftim:
   standard_name : modified_red_flag_threat_index
   units         : index
   rank          : 1
   type          : real
   kind          : kind_phys

Let's register this variable in the GFS diagnostics module. Open FV3/gfsphysics/GFS_layer/GFS_diagnostics.F90 and add the following blob to the end of routine GFS_externaldiag_populate (be sure to enter this after the last #endif preprocessor directive):

   idx = idx + 1
   ExtDiag(idx)%axes = 2
   ExtDiag(idx)%name = 'rftim'
   ExtDiag(idx)%desc = 'modified red flag threat index'
   ExtDiag(idx)%unit = 'n/a'
   ExtDiag(idx)%mod_name = 'gfs_sfc'
   ExtDiag(idx)%intpl_method = 'bilinear'
   allocate (ExtDiag(idx)%data(nblks))
   do nb = 1,nblks
     ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%rftim(:)
   enddo

Calculate the index

Now we can begin the task of calculating the RFTIM (variable name rftim in the following). A skeleton CCPP scheme called rtfim, as well as a module containing the simplified climatology to map rh2m and ws6m to red flag threat indices, are provided for you to copy into your development codebase.  Both files need to be copied into ccpp/physics/physics:

      /scratch4/BMC/gmtb/emc_training_march_2019/handson_day2/templates/module_rftim_clima.F90
      /scratch4/BMC/gmtb/emc_training_march_2019/handson_day2/templates/rftim.F90

Please take a look at module_rftim_clima.F90, but note that you do not need to make any modification in this file. You will see two functions in this file. The functions rfti_rh_clima and rfti_ws_clima can be called for arrays or scalars, and return integers. The interfaces are for relative humidity and wind speed are:

   elemental function rfti_rh_clima(lat, lon, rh)  result(rfti_rh)
      real(kind_phys), intent(in) :: lat
      real(kind_phys), intent(in) :: lon
      real(kind_phys), intent(in) :: rh      
      integer :: rfti_rh
   end function rfti_rh_clima

   elemental function rfti_ws_clima(lat, lon, ws) result(rfti_ws)
      real(kind_phys), intent(in) :: lat
      real(kind_phys), intent(in) :: lon
      real(kind_phys), intent(in) :: ws
      integer :: rfti_ws
   end function rfti_ws_clima

Next you need to edit rftim.F90 to compute RFTIM.

Note that in file rftim.F90 there are 3 subroutines wrapped in a module: rftim_init, rftim_finalize, and rftim_run. This is typical of CCPP schemes.  No changes are needed in rftim_init and rftim_finalize, but you need to edit rftim_run and do the following:

  • Update the metadata table in the module

  • Declare local variable i

  • Calculate RFTIM

The argument list for subroutine rftim_run is:

  subroutine rftim_run(im, lat, lon, islmsk, rh2, ws6, rftim, errmsg, errflg)

The metadata table needs to be populated with the same variables and in the same order as the subroutine call. The local variable names, standard names, intent, and optional columns are already pre-populated. Your job is to add the remaining metadata. Note that units are critical because they are used during the build system to compare the variables requested by the physics against those supplied by the host model. Therefore, the units needs to be exactly what the host model provides.

Variable im:
   standard_name : horizontal_loop_extent
   units         : count
   rank          : 0
   type          : integer
   kind          : 

Variable lat:
   standard_name : latitude
   units         : radians
   rank          : 1
   type          : real
   kind          : kind_phys

Variable lon:
   standard_name : longitude
   units         : radians
   rank          : 1
   type          : real
   kind          : kind_phys

Variable islmsk:
   standard_name : sea_land_ice_mask
   units         : flag
   rank          : 1
   type          : integer
   kind          : 

Variable rh2m:
   standard_name : relative_humidity_at_2m
   units         : percent
   rank          : 1
   type          : real
   kind          : kind_phys

Variable ws6m:
   standard_name : wind_speed_at_6m_in_miles_per_hour
   units         : mi h-1
   rank          : 1
   type          : real
   kind          : kind_phys

Variable rftim:
   standard_name : modified_red_flag_threat_index
   units         : index
   rank          : 1
   type          : real
   kind          : kind_phys

Variable errmsg:
   standard_name : ccpp_error_message
   units         : none
   rank          : 0
   type          : character
   kind          : len=*

Variable errflg:
   standard_name : ccpp_error_flag
   units         : flag
   rank          : 0
   type          : integer
   kind          : 

Hint: For unit “mi h-1” make sure you have the white space in the middle. If not, there will be a mismatch between the variable provided by the host model and what scheme rtfim requests.

Now that you are done with the metadata table, add to the subroutine

integer :: i

do i=1,im

        rftim(i) = real(rfti_rh_clima(lat(i), lon(i), rh2(i)) + rfti_ws_clima(lat(i), lon(i), ws6(i)), kind_phys)

enddo

Scientific Documentation

Scientific Documentation carson Wed, 03/06/2019 - 12:58

Prepare the Code

In this part of the training you will see how to add documentation to the new scheme, rtfim. The documentation uses Doxygen markup, so it can later be converted for display as an html or PDF document.  The doxygen markup have been added to the template. At the end of this section, students are expected to

  • finish the metadata table 
  • insert the calculation of RFTIM 
  • add the citation in library.bib file.

Subroutines rftim_init() and rftim_finalize() are just empty placeholders, we do not need to document them. Since rftim_run() is the entry point of the RFTIM scheme, Doxygen documentation should include:

1. In the first line of the Fortran file, add a brief one-sentence overview of the file purpose following “\file”, e.g.,

!>\file rftim.F90
!! This file contains the modified Red Flag Threat Index calculation.

The Doxygen code block begins with “!>”, and subsequent lines begin with “!!”.

2. Define a parent module using “\defgroup”, e.g., this module implements a particular scheme functionality. Also add a  one-line description of the scheme with reference using “\cite” if applicable.  Add the following lines after “end subroutine rftim_finalize”

!>\defgroup rftim_cal The Modified Red Flag Threat Index calculation
!!\brief This module calculates RFTIM using the 2-m relative humidity in percent
!! and 6m wind speed in miles per hour Lindley et al. (2011) \cite lindley_et_al_2011.

Doxygen can handle in-line paper citations and link to an automatically created bibliography page. A ccpp/physics/physics/docs/library.bib file in BibTex format is provided in the repository. Insert  the following item in the file (order is not important):

          @article {lindley_et_al_2011,
                Author  = {T.T. Lindley and J.D. Vitale and W.S. Burgett and M.-J. Beierle},
                Title      = {Proximity meteorological observations for wind-driven grassland wildfire starts on the southern High Plains},
                Journal = {Electronic J. Severe Storms Meteor.},
                Year     = {2011},
                Volume= {6},
                Number  = {1},
                Pages   = {1-27}}

3. In rftim_run(), which is the entry point to the scheme, further documentation will include:

     An argument table with:

  • Local variable name “local_name”
  • CF-standard compliant “standard_name”: this is used as a key to communicate variables between the CCPP and the host model
  • Short description “long_name”
  • Units (format follows “unit exponent”, i.e., m2 s-2 for m2/s2)
  • Rank (0 for scalar, 1 for 1-D array, 2 for 2-D array, etc)

 

Part of the argument table is intentionally left blank. Student are expected to fill in the remaining metadata.

     The scheme general algorithm:

  • Start with “\section sectionID XYZ General Algorithm”
  • List in-line calculation steps by using “-#” markers

 

4. For those child subroutines that are not an entry point to the scheme, we can associate a subroutine or function with its parent module by using “\ingroup” to enable the auto-generation of a  calling hierarchy within the scheme.

!>\ingroup rftim_cal
!!This function calculates the climatological relative humidity to integer.
!!...(detailed description could go here)...
!!\param lon  real, longitude in radiance
!!\param rh   real, relative humidity in percentage
!!\param rfti_rh  integer, the climatological relative humidity
    elemental function rfti_rh_clima(lat,lon,rh) result(rfti_rh)

 

Following a brief description of this function or subroutine, use “\param” to define each parameter with local name, provide a short description and list the unit (optional but highly recommended).

 

Run Doxygen

 

  • Get set up to use Doxygen v1.8.10 by adding  the following line into .cshrc file under your home directory:

alias doxygen /scratch4/BMC/gmtb/doxygen-1.8.10/bin/doxygen

  • Source your .cshrc file.
  • Edit file ccpplatex_dox and add the new files to the INPUT variable: 

                        ../calpreciptype.f90   \
### EMC_CCPP_TRAINING rftim
                        ../rftim.F90       \
                        .. module_rftim_clima.F90    \

Here RFTIM scheme is inserted after calpreciptype.f90

Under ccpp/physics/physics/docs,  run Doxygen by typing:

bash:

doxygen ccpplatex_dox 2>&1 | tee ccppdoc.log

tcsh:

doxygen ccpplatex_dox | & tee ccppdoc.log

  • The generated html  will be located under ccpp-physics/physics/doc/html/.
  • In order to visualize the output, later you can transfer the html/ directory to your local machine. It can be viewed by pointing a browser to the index.html file.

 

Add the scheme to CCPP

Add the scheme to CCPP carson Mon, 03/04/2019 - 15:00

Add the scheme to CCPP

Once the CCPP scheme rftim is complete, it is time to add the scheme rftim.F90 and its prerequisite module_rftim_clima.F90 to the CCPP prebuild configuration in ccpp/config/ccpp_prebuild_config.py.

Hint: Follow the procedure outlined in talk day2_3_prebuild.ppt section “Modifying CCPP prebuild config”. Note that scheme rftim will be added to the “slow physics” group.  rftim.F90 goes in SCHEME_FILES and module_rftim_clima.F90 goes in SCHEME_FILES_DEPENDENCIES.

 

The next step is to create a new Suite Definition File and edit it so that it runs the scheme rftim just after the sfc_diag scheme:

cd ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/develop

cp ccpp/suites/suite_FV3_GFS_2017_updated_gfdlmp.xml ccpp/suites/suite_FV3_GFS_2017_rftim.xml

Edit ccpp/suites/suite_FV3_GFS_2017_rftim.xml and

  • change the suite name to 'FV3_GFS_2017_rftim'
  • add the scheme "rftim" after sfc_diag.

Copy the development directory to a test directory and change into that directory:

   rsync -av ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/develop/ ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/compile/

   cd ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/compile/tests

Compile the code (this example is for a dynamic build since the STATIC variable is not specified):

  • for bash shell:

  ./compile.sh $PWD/../FV3 theia.intel “CCPP=Y HYBRID=N” '' NO NO 2>&1 | tee compile.log

  • for t/csh shell:

  ./compile.sh $PWD/../FV3 theia.intel “CCPP=Y HYBRID=N” '' NO NO  |& tee compile.log

Hint 1: After HYBRID=NO, there is a double quote plus two single quotes.

Hint 2: In file compile.log, look for the following statement: INFO: CCPP prebuild step completed successfully. If not found, some problem occurred in the CCPP prebuild stage. If found, prebuild was successful and you can look at other messages in compile.log pertaining to prebuild.

If the entire prebuild and build process was successful, this creates an executable 'fv3.exe' in the same directory.

  • The compile.sh script takes up to 6 arguments. 

PATHTR BUILD_TARGET [ MAKE_OPT [ BUILD_NR ] [ clean_before ] [ clean_after ]

  • so, in the example command-line given above,
    • PATHTR = $PWD/../FV3
    • BUILD_TARGET = theia.intel
    • MAKE_OPT = "CCPP=Y HYBRID=N"
    • BUILD_NR (name of the fv3.exe executable) = null (so no additional string, just fv3.exe)
    • clean_before = NO
    • clean_after = NO

If successful, this creates an executable 'fv3.exe' in the same directory.

 

Run the model

Run the model carson Mon, 03/04/2019 - 15:00

Run the model

Now copy the run directory template to your directory:

   rsync -av /scratch4/BMC/gmtb/emc_training_march_2019/handson_day2/fv3_rftim/ ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/run/

   cd ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/run/

Edit job_card to:

  • Adjust the account name in line "# PBS -A ...".
  • Adjust CODEDIR accordingly - this will also adjust the LD_LIBRARY_PATH later in the job_card script.

Take a quick look at the additional steps required for CCPP:

  • Link the Suite Definition File you previously created (develop/ccpp/suites/suite_FV3_GFS_2017_rftim.xml) to the run directory

ln -s ${YOUR_WORK_DIR}/develop/ccpp/suites/suite_FV3_GFS_2017_rftim.xml ccpp_suite.xml

  • Edit file input.nml and set it to use your Suite Definition File

ccpp_suite = 'ccpp_suite.xml'

Note that in this example the CCPP library will be loaded at runtime since the compilation command was for the dynamic CCPP build.

Before we run the case, we need to add the new diagnostic variable to the list in file diag_table. Edit diag_table and add the following line just below the 'qrain' line:

  "gfs_sfc",     "rftim", "rftim",     "fv3_history2d", "all", .false.,  "none", 2

Submit the job script

 

   qsub job_card && touch out && tail -f out

When finished, plot the new diagnostic variable rftim:

   module load intel/15.1.133 ncl/6.5.0
   ncl plot_rftim.ncl

Contribute your changes to GitHub

Contribute your changes to GitHub carson Mon, 03/04/2019 - 15:02

Assuming you are happy with your development, commit the changes and push them to your GitHub fork.

cd ${YOUR_WORK_DIR}/NEMSfv3gfs_CCPP_training/develop

   # ccpp-physics
   cd ccpp/physics
   git add .
   git commit -m "Add scheme rftim and dependencies"
   git push origin my_rftim_implementation
   cd ../..
   # FV3
   cd FV3
   git add .
   git commit -m "Add variables and configure diagnostic output for RFTIM calculation"
   git push origin my_rftim_implementation
   cd ..
   # NEMSfv3gfs
   git add .
   git commit -m "Add scheme rftim and dependencies to CCPP prebuild config, add SDF"
   git push origin my_rftim_implementation

At this point, your changes have been committed to your forks of the repositories.

To propose these changes for inclusion into the GMTB authoritative repository on NCAR's GitHub organization, you need to create three pull requests (one for each of the repositories). For the example of ccpp-physics:

  1. Navigate to http://github.com/NCAR/ccpp-physics and click on 'new pull request'
  2. Select 'compare across forks', set base repository to 'NCAR/ccpp-physics', base to 'master', head repository to 'YOUR_GITHUB_USERNAME/ccpp-physics', compare to 'my_rftim_implementation'
  3. Add a descriptive title and short description in the text boxes
  4. This concludes the exercise
  5. If this were a real development for inclusion in master, you would now click on 'Create pull request' . PLEASE DO NOT DO THAT FOR THIS EXERCISE!