NOAA-OWP / noah-owp-modular

Modularized version of the NOAH-MP land surface model.
Other
8 stars 19 forks source link

homogenous gridded driver for noah-owp-modular #71

Closed GreyREvenson closed 1 year ago

GreyREvenson commented 1 year ago

PURPOSE

This PR is meant accomplish the first of three steps/milestones aimed at creating a gridded version of NOAA-NWS OWP's Noah-mp model (i.e., 'Noah-owp'). Step 1 is:

  1. Develop Noah-owp capacity to execute over a homogenous gridded domain

Steps 2 and 3 will be accomplished via subsequent PRs. Steps 2 and 3 are:

  1. Develop Noah-owp capacity to execute over a heterogenous gridded domain
  2. Develop Noah-owp capacity to execute using real-word spatial data

ADDITIONS

This PR adds a series of new 'grid' derived data types, each defined in a distinct module, i.e.:

The PR-added 'gridded' derived data types correspond the original derived data types within the Noah-owp column model, i.e.:

The PR-added 'gridded' derived data types contain the same variable names that exist within their corresponding original derived data types (e.g., both water_type and watergrid_type have a variable named qinsur); however, the 'gridded' version of the variables (e.g., qinsur in watergrid_type) have x and y dimensional indices if the variable varies over those dimensions. For example, qinsur is defined as real,allocatable,dimension(:,:) :: qinsur in watergrid_type but is defined as real :: qinsur in water_type.

The PR also adds noahowpgrid_type to RunModule.f90. noahowpgrid_type is composed of the new 'gridded' derived data types and the namelist_type, which is mostly unchanged by this PR. Thus, noahowpgrid_type is defined as:

type :: noahowpgrid_type

type(namelist_type)       :: namelist
type(levelsgrid_type)     :: levelsgrid
type(domaingrid_type)     :: domaingrid
type(optionsgrid_type)    :: optionsgrid
type(parametersgrid_type) :: parametersgrid
type(watergrid_type)      :: watergrid
type(forcinggrid_type)    :: forcinggrid
type(energygrid_type)     :: energygrid

end type

This PR leaves the noahowp_type (also in Runmodule.f90) unchanged except for the exclusion of the namelist_type. Thus, this PR defines noahowp_type as:

type :: noahowp_type

type(levels_type)     :: levels
type(domain_type)     :: domain
type(options_type)    :: options
type(parameters_type) :: parameters
type(water_type)      :: water
type(forcing_type)    :: forcing
type(energy_type)     :: energy

end type noahowp_type

This PR also adds a series of subroutines to copy/move or 'transfer' model variable values between the 'gridded' and original versions of the derived data types, i.e.:

This PR also adds a new subroutine called solve_noahowp_grid to RunModule.f90. The primary purpose of the solve_noahowp_grid subroutine is to iterate over the x and y dimensions of a gridded domain, transfer model variable values that correspond to the iterated grid cell from the noahowpgrid_type to an instance of noahowp_type via calls to transfer subroutines, then execute the Noah-owp column model for the iterated grid cell via a call to the solve_noahowp subroutine, and finally transfer model variable values from noahowp_type back to noahowpgrid_type.

CHANGES

This PR modifies the initialize_from_file subroutine in RunModule.f90 to initialize an instance of noahowpgrid_type. Previous to this PR, the initialize_from_file subroutine instead initialized an instance of noahowp_type. The PR-modified initialize_from_file subroutine first reads namelist.input via a call to namelist%ReadNamlist (i.e., call namelist%ReadNamelist(config_filename)). The initialize_from_file subroutine then makes a series of calls to type-bound initialization subroutines and subroutines that transfer namelist.input values to noahowpgrid_type and its member derived data types, i.e.:

call levelsgrid%Init(namelist)
call levelsgrid%InitTransfer(namelist)

call domaingrid%Init(namelist)
call domaingrid%InitTransfer(namelist)

call optionsgrid%Init(namelist)
call optionsgrid%InitTransfer(namelist)

call parametersgrid%Init(namelist)
call parametersgrid%paramRead(namelist,domaingrid)

call forcinggrid%Init(namelist)
call forcinggrid%InitTransfer(namelist)

call energygrid%Init(namelist)
call energygrid%InitTransfer(namelist)

call watergrid%Init(namelist)
call watergrid%InitTransfer(namelist)

After initializing and transfering namelist.input varaible values to noahowpgrid_type, the PR-modified initialize_from_file subroutine tracks closely with the pre-PR initialize_from_file subroutine, i.e.: it initializes important model variables to zero, sets time variables, opens the input forcing file, and finally initializes the output.nc file.

This PR removes initialization and transfer subroutines that were previously used to initialize and then transfer namelist.input variable values to noahowp_type. These subroutines were recycled/extended to initialize and then transfer namelist.input values to noahowpgrid_type instead and are now type-bound subroutines within each of the 'gridded' derived data types, i.e.:

This PR adds x and y dimensions to output.nc via modifications to OutputModule.f90.

This PR updates /test/analysis/noahowp_driver_test.f90 to work on the PR-modified code base as a means of testing the PR additions and changes (see TESTING section).

TESTING

The PR-modified Noah-owp model was successfully compiled on MacOS using GNU Fortran (Homebrew GCC 13.1.0) 13.1.0. Then, two tests were completed:

  1. A python script (compare_outputnc.py) was created to read-in the PR-modified model's output netcdf file (which gives output for every grid cell) as well as the output netcdf file from a clean build of the Noah-owp column model. The script checks that output for every grid cell in the PR-modified model matches the output given by the clean build for every time step and every variable in output.nc. As this PR assumes a homogenous grid, the output from every grid cell should be the same and match that of the column model, assuming the same namelist.input is used to initialize. The output from the py script is output_compare_outputnc.txt

  2. noahowp_driver_test.f90 was updated to work with the PR-modified code base. The output from the test program is test_output.txt

NOTES

Concern regarding performance: The general architectural approach of this PR is to have an object/type hold all model variables across the x, y, and z dimensions. For each time step, the grid is iterated in x and y dimensions, the iterated grid cell's values are transferred to a column model object/type, the column model is executed, then the model variables are transferred back to the grid object/type before moving to the next grid cell. This is the approach employed by NCAR's NoahmpDriverMainMod.f90 (in HRLDAS). This architectural approach requires a lot of copying/moving of variable values. Parallel processing should improve performance.

All model variables are included in the transfer subroutines to be on safe side. However, we could improve performance by commenting out / removing variables that are being unnecessarily transferred in these subroutines?

./EtFluxModule.f90, line 708: energy%LATHEA was changed to energy%LATHEAV because energy%LATHEA was set to huge(1.) and causing arithmetic error. Is this change acceptable/scientifically valid?

Checklist

SnowHydrology commented 1 year ago

@GreyEvenson-NOAA, please provide an illustrative title for the pull request and include a description (see https://github.com/NOAA-OWP/noah-owp-modular/pull/57 as an example, or this from CFE is even better https://github.com/NOAA-OWP/cfe/pull/62).

SnowHydrology commented 1 year ago

Next steps:

We also need to update the README.md file on the main level. This includes an updated link to the old Noah-MP codebase and new acknowledgments of NCAR's work on the grid routines (consistent with their license).

GreyREvenson commented 1 year ago

I believe we're ready for you to have another look, @SnowHydrology.

In my revisions to the parametersgrid_type member variables, I tried to error on the side of including x and y dim indices as I believe you mentioned in a call that member variables that do not include x and y dim indices would be the exception rather than the rule. However, please let me know if I misunderstood.

Following all revisions to the PR, I re-ran /test/analysis/compare_outputnc.py and its output was: output_compare_outputnc.txt. I also re-ran /test/noahowp_test.exe and its output was: noahowp_test_output.txt

I clicked on the button to 'Resolve conversation' for your comments that were easy fixes with little ambiguity. However, please let me know if I should leave those conversations open or have other feedback regarding my use of github.