diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 00000000..871f045a --- /dev/null +++ b/NEWS.md @@ -0,0 +1,3 @@ +## v4.0 + +* Public release after refactored code \ No newline at end of file diff --git a/README.md b/README.md index 1f6ac1ac..9fffe140 100644 --- a/README.md +++ b/README.md @@ -3,21 +3,11 @@ # rsofun -Provides a modelling framework that implements the P-model for leaf-level acclimation of photosynthesis for site-scale simulations, with SPLASH used for simulating the soil water balance (see also [Stocker et al., 2019 GMDD](https://www.geosci-model-dev-discuss.net/gmd-2019-200/)). The package provides the following functionalities: +A modelling framework for site-scale simulations of ecosystem processes, implemented as an R package (low-level routines in Fortran 90). Implements the following models: -- Calibrating model parameters -- Running the model and getting outputs directly back into R ("tidy" data) -- Evaluating outputs (benchmarking) - -Model forcing and calibration data is collected using the [ingestr](https://stineb.github.io/ingestr/) package. See [here](https://rpubs.com/stineb/rsofun) for an example. - -Parallelisation for a large number of site-level simulations is provided using the *multidplyr* R package. Calibration uses the *GenSA* R package. - -The P-model is implemented in different repositories for different purposes: - -- **rsofun** (this package): Is for site-scale simulations (large ensemble of sites can be run in parallelised mode), forced by time series of meteorological data. Acclimation of photosynthesis is assumed at a user-defined time scale. I.e., the P-model optimality criterion (Wang et al., 2017; Prentice et al., 2014) is solved daily with "damped" daily variations in the forcing data (similar to a running mean). All model code is implemented in Fortran. -- [**rpmodel**](https://stineb.github.io/rpmodel/): Implements the same equation as rsofun (all native R), but solves the optimality criterion for each time step independently (instantaneous acclimation). This can be used for hypothesis generation, exploration, and illustrations. Transient simulations of acclimation and GPP should be done using rsofun. -- [**SOFUN**](https://stineb.github.io/sofun/): This is for P-model simulations on a (global) spatial grid and is purely in Fortran. Forcing data is read from NetCDF files and outputs are written to NetCDF files. +- P-model for leaf-level acclimation of photosynthesis from [Stocker et al. (2019)](https://www.geosci-model-dev-discuss.net/gmd-2019-200/). +- SPLASH for bioclimatic variables, including the surface radiation budget and the soil water balance from [Davis et al. (2017)](https://doi.org/10.5194/gmd-10-689-2017). +- LM3-PPA for comprehensive simulations of ecosystem carbon and water cycling, tree growth, and tree cohort-explicit forest dynamics following the Perfect Plasticity Approximation, from [Weng et al., (2015)](https://doi.org/10.5194/bg-12-2655-2015). ## Installation diff --git a/docs/articles/benchmark_gpp_FLUXNET2015_ensemble_EXAMPLE_files/header-attrs-2.7/header-attrs.js b/docs/articles/benchmark_gpp_FLUXNET2015_ensemble_EXAMPLE_files/header-attrs-2.7/header-attrs.js new file mode 100644 index 00000000..dd57d92e --- /dev/null +++ b/docs/articles/benchmark_gpp_FLUXNET2015_ensemble_EXAMPLE_files/header-attrs-2.7/header-attrs.js @@ -0,0 +1,12 @@ +// Pandoc 2.9 adds attributes on both header and div. We remove the former (to +// be compatible with the behavior of Pandoc < 2.8). +document.addEventListener('DOMContentLoaded', function(e) { + var hs = document.querySelectorAll("div.section[class*='level'] > :first-child"); + var i, h, a; + for (i = 0; i < hs.length; i++) { + h = hs[i]; + if (!/^h[1-6]$/i.test(h.tagName)) continue; // it should be a header h1-h6 + a = h.attributes; + while (a.length > 0) h.removeAttribute(a[0].name); + } +}); diff --git a/docs/articles/example_lm3ppa_files/header-attrs-2.7/header-attrs.js b/docs/articles/example_lm3ppa_files/header-attrs-2.7/header-attrs.js new file mode 100644 index 00000000..dd57d92e --- /dev/null +++ b/docs/articles/example_lm3ppa_files/header-attrs-2.7/header-attrs.js @@ -0,0 +1,12 @@ +// Pandoc 2.9 adds attributes on both header and div. We remove the former (to +// be compatible with the behavior of Pandoc < 2.8). +document.addEventListener('DOMContentLoaded', function(e) { + var hs = document.querySelectorAll("div.section[class*='level'] > :first-child"); + var i, h, a; + for (i = 0; i < hs.length; i++) { + h = hs[i]; + if (!/^h[1-6]$/i.test(h.tagName)) continue; // it should be a header h1-h6 + a = h.attributes; + while (a.length > 0) h.removeAttribute(a[0].name); + } +}); diff --git a/docs/articles/example_pmodel_files/header-attrs-2.7/header-attrs.js b/docs/articles/example_pmodel_files/header-attrs-2.7/header-attrs.js new file mode 100644 index 00000000..dd57d92e --- /dev/null +++ b/docs/articles/example_pmodel_files/header-attrs-2.7/header-attrs.js @@ -0,0 +1,12 @@ +// Pandoc 2.9 adds attributes on both header and div. We remove the former (to +// be compatible with the behavior of Pandoc < 2.8). +document.addEventListener('DOMContentLoaded', function(e) { + var hs = document.querySelectorAll("div.section[class*='level'] > :first-child"); + var i, h, a; + for (i = 0; i < hs.length; i++) { + h = hs[i]; + if (!/^h[1-6]$/i.test(h.tagName)) continue; // it should be a header h1-h6 + a = h.attributes; + while (a.length > 0) h.removeAttribute(a[0].name); + } +}); diff --git a/docs/index.html b/docs/index.html index e45431cd..7a0803d6 100644 --- a/docs/index.html +++ b/docs/index.html @@ -98,24 +98,6 @@

devtools::install_github("stineb/rsofun") library(rsofun) -<<<<<<< HEAD -
-

-Environment

-

This is important

-

In order to successfully compile the package (Fortran source), you need to have gfortran installed and manually adjust the compiler flag specification. To do so, open the Makeconf file (you’ll find it by entering in R: file.path(R.home("etc"), "Makeconf")). In there, add the gfortran flag -ffree-line-length-0. The respective line then looks like this:

-
FCFLAGS = -Wall -g -O2 $(LTO) -ffree-line-length-0
-
-
-

-Dependencies

-

The rsofun package requires a large number of other R-packages (dependencies). Required dependencies are essential for rsofun functions and are:

- -
-======= ->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674

Example

@@ -125,41 +107,6 @@

Usage and contribution

This package is designed to be extendible to ingesting other data types (sources). The developer (Beni Stocker) would appreciate if you made sure that your developments can be fed back to this repository. To do so, please use git. See here for a brief introduction to git.

-<<<<<<< HEAD -
-

-Application only

-

I recommend the following steps if you would just like to use this package (no development):

-
    -
  • Install and load the library as described under ‘Installation’ above.
  • -
-
-
-

-For developers

-

I recommend the following steps if you would like to use and further develop the package (even if this is just some extension for your own application - Keep in mind: Others may benefit from your efforts too!):

-
    -
  1. Make sure you have a Github account.
  2. -
  3. Log on to Github, and go to https://github.com/stineb/rsofun and click on ‘Fork’ in the upper right corner. This makes a copy of the repository which then belongs to you, meaning that you can modify, commit, and push changes back to your forked repository as you please.
  4. -
  5. Clone your fork to your local computer by entering in your terminal (here, it’s cloned to a subdirectory ingestr placed in your home directory):
  6. -
-
cd home
-git clone https://github.com/<your_github_username>/rsofun.git
-
    -
  1. In RStudio, create a new project in your local directory ~/rsofun/. This opens the repository in RStudio and you have access to the code where all the functions of this package are implemented (see subdirectory ./R/).
  2. -
  3. In RStudio, after having edited code, select the ‘Build’ tab and click on ‘Install and Restart’ to build the package again. For quick edits and checks, you may simply source the edited files instead of re-building the whole package. If you like to add new functions, create a new source file that contains your function in subdirectory ./R/, write a nice roxygen header (see other source files as an example), then click on ‘Build’ -> ‘More’ -> ‘Document’, and then again on ‘Install and Restart’.
  4. -
  5. You can upload (commit and push) your edits and additions to your forked repository by
  6. -
-
git add -u  # adds all edits to your next commit
-git add <newfile>  # adds new file to the git repository
-git commit -m "a brief description of what you did"
-git push  # pushes the commit to your fork 
-
    -
  1. If you’re happy with your new edits and additions to the package, you may want to have it fed back from your fork to the original repository. To do so, please create a new pull request in GitHub: Click on ‘New pull request’ on the repository page and follow the inuitive steps. Thanks!
  2. -
-
-======= ->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index ac321f37..a330ec95 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,15 +1,6 @@ pandoc: 2.11.2 pkgdown: 1.6.1 pkgdown_sha: ~ -<<<<<<< HEAD -articles: - benchmark_gpp_FLUXNET2015_ensemble_EXAMPLE: benchmark_gpp_FLUXNET2015_ensemble_EXAMPLE.html - example_lm3ppa: example_lm3ppa.html - example_pmodel: example_pmodel.html - prepare_inputs_rsofun: prepare_inputs_rsofun.html -last_built: 2021-07-30T16:24Z -======= articles: {} -last_built: 2021-08-04T05:51Z ->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674 +last_built: 2021-09-27T11:51Z diff --git a/docs/reference/calc_net_assim.html b/docs/reference/calc_net_assim.html new file mode 100644 index 00000000..2eb05772 --- /dev/null +++ b/docs/reference/calc_net_assim.html @@ -0,0 +1,178 @@ + + + + + + + + +Calculate net assimilation — calc_net_assim • rsofun + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+ + +
+

Net assimilation calculation

+
+ +
calc_net_assim(par, args, iabs, kphio, a_unitcost)
+ +

Arguments

+ + + + + + + + + + + + + + + + + + + + + + +
par

parameters (vcmax, gs)

args

arguments (kmm, gammastart, ns_star, ca, vpd, beta)

iabs

iabs

kphio

kphio

a_unitcost

cost unit

+ +

Value

+ +

Net assimilation

+ +
+ +
+ + + +
+ + + + + + + + diff --git a/docs/reference/calib_sofun.html b/docs/reference/calib_sofun.html index 7238f047..39ba96ce 100644 --- a/docs/reference/calib_sofun.html +++ b/docs/reference/calib_sofun.html @@ -146,10 +146,6 @@

Value

the calibration settings and optimised parameter values.

Examples

-<<<<<<< HEAD -
if (FALSE) calib_sofun( df_drivers = filter(df_drivers, sitename %in% settings$sitenames), ddf_obs = ddf_obs_calib, settings = settings) - -=======
if (FALSE) { calib_sofun( df_drivers = filter(df_drivers, @@ -157,7 +153,6 @@

Examp ddf_obs = ddf_obs_calib, settings = settings) } ->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674

-<<<<<<< HEAD -
runread_lm3ppa_f(df_drivers, makecheck = TRUE, parallel = FALSE, ncores = 2)
-=======
runread_lm3ppa_f(drivers, makecheck = TRUE, parallel = FALSE, ncores = 2)
->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674

Arguments

-<<<<<<< HEAD - - -======= ->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674 @@ -153,10 +143,6 @@

Arg

- - - -
df_drivers

A nested data frame with one row for each site and columns named according to the -arguments of function runread_pmodel_f_bysite

drivers

A nested data frame with one row for each site and columns named according to the arguments of function `runread_pmodel_f_bysite()`

makecheck

An integer specifying the number of cores used for parallel computing. Defaults to 2.

params_modl

A named list of model parameters

Value

@@ -165,20 +151,11 @@

Value

in the nested column data

Examples

-<<<<<<< HEAD -
mod <- runread_pmodel_f( df_drivers, params_modl, makecheck = TRUE, parallel = FALSE, ncores = 2 ) -
#> Error: Problem with `mutate()` input `data`. -#> x Can't subset columns that don't exist. -#> x Column `rain` doesn't exist. -#> Input `data` is `purrr::pmap(...)`.
-
-=======
if (FALSE) { mod <- runread_pmodel_f( df_drivers, params_modl, makecheck = TRUE, parallel = FALSE, ncores = 2 ) }
->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674 -
runread_pmodel_f(drivers, par, makecheck = TRUE, parallel = FALSE, ncores = 2)
+
runread_pmodel_f(drivers, par, makecheck = TRUE, parallel = FALSE, ncores = 1)

Arguments

@@ -155,14 +155,6 @@

Value

in the nested column data

Examples

-<<<<<<< HEAD -
mod <- runread_pmodel_f( df_drivers, params_modl, makecheck = TRUE, parallel = FALSE, ncores = 2 ) -
#> Error: Problem with `mutate()` input `data`. -#> x Can't subset columns that don't exist. -#> x Column `rain` doesn't exist. -#> Input `data` is `purrr::pmap(...)`.
-
-=======
if (FALSE) { mod <- runread_pmodel_f( drivers, @@ -172,7 +164,6 @@

Examp ncores = 2 ) }

->>>>>>> 463af7b984d4541e03a5994c8b0069d37aadf674
+ + + + + + + + + +
object

input data generated using the run_*_f_bysite() function +or generally calls which return the rsofun_fit class

...

additional parameters to pass

+ +

Value

+ +

a table with fit statistics, a data frame with summary statistics

+ + + + + + + + + + + + + + + + diff --git a/src/biosphere_pmodel.mod.f90 b/src/biosphere_pmodel.mod.f90 index 192cedc1..f084399b 100644 --- a/src/biosphere_pmodel.mod.f90 +++ b/src/biosphere_pmodel.mod.f90 @@ -34,7 +34,6 @@ function biosphere_annual() result( out_biosphere ) ! contact: b.stocker@imperial.ac.uk !---------------------------------------------------------------- use md_interface_pmodel, only: myinterface, outtype_biosphere - use md_sofunutils, only: daily2monthly ! return variable type(outtype_biosphere) :: out_biosphere diff --git a/src/gpp_lm3ppa.mod.f90 b/src/gpp_lm3ppa.mod.f90 index 0fe1fb17..43bfe81b 100644 --- a/src/gpp_lm3ppa.mod.f90 +++ b/src/gpp_lm3ppa.mod.f90 @@ -619,8 +619,6 @@ subroutine getpar_modl_gpp() !//////////////////////////////////////////////////////////////// ! Subroutine reads module-specific parameters from input file. !---------------------------------------------------------------- - use md_sofunutils, only: getparreal - ! local variables integer :: pft diff --git a/src/gpp_pmodel.mod.f90 b/src/gpp_pmodel.mod.f90 index 83db443f..9292f064 100644 --- a/src/gpp_pmodel.mod.f90 +++ b/src/gpp_pmodel.mod.f90 @@ -510,8 +510,6 @@ subroutine getpar_modl_gpp() !//////////////////////////////////////////////////////////////// ! Subroutine reads module-specific parameters from input file. !---------------------------------------------------------------- - use md_sofunutils, only: getparreal - ! local variables integer :: pft diff --git a/src/plant_pmodel.mod.f90 b/src/plant_pmodel.mod.f90 index 53f9fa8f..cc6a83d9 100644 --- a/src/plant_pmodel.mod.f90 +++ b/src/plant_pmodel.mod.f90 @@ -180,8 +180,6 @@ subroutine getpar_modl_plant() ! Copyright (C) 2015, see LICENSE, Benjamin David Stocker ! contact: b.stocker@imperial.ac.uk !---------------------------------------------------------------- - use md_sofunutils, only: getparreal - ! local variables integer :: pft integer :: npft_site @@ -237,8 +235,6 @@ function getpftparams( pftname ) result( out_getpftparams ) !---------------------------------------------------------------- ! Read PFT parameters from respective file, given the PFT name !---------------------------------------------------------------- - use md_sofunutils, only: getparreal - ! arguments character(len=*), intent(in) :: pftname diff --git a/src/sofunutils.mod.f90 b/src/sofunutils.mod.f90 index 82e1f54d..06e1b123 100644 --- a/src/sofunutils.mod.f90 +++ b/src/sofunutils.mod.f90 @@ -87,846 +87,846 @@ function running( presval, inow, lenval, lenper, method, prevval ) result( runni end function running - function daily2monthly( dval, method ) result( mval ) - !///////////////////////////////////////////////////////////////////////// - ! Returns monthly values as a mean over daily values in each month. - !------------------------------------------------------------------------- - use md_params_core, only: ndaymonth, cumdaymonth, ndayyear, nmonth - - ! arguments - real, intent(in), dimension(ndayyear) :: dval ! vector containing 365 (366 in case lapyear is TRUE) daily values - character(len=*), intent(in) :: method ! true of monthly values represent total of daily values in resp. month - - ! function return variable - real, dimension(nmonth) :: mval - - ! local variables - integer :: moy - integer, dimension(nmonth) :: istart, iend - - istart = cumdaymonth - ndaymonth + 1 - iend = cumdaymonth - - ! loop over months and take sum/mean of days in that month - do moy=1,nmonth - if (method=="sum") then - mval(moy) = sum( dval( istart(moy) : iend(moy) )) - else if (method=="mean") then - mval(moy) = sum( dval( istart(moy) : iend(moy) )) / ndaymonth(moy) - else - stop 'DAILY2MONTHLY: select valid method (sum, mean)' - end if - end do - - end function daily2monthly - - - function monthly2daily_weather( mval_prec, mval_wet, prdaily_random ) result( dval_prec ) - !///////////////////////////////////////////////////////////////////////// - ! Weather (rain) generator. - ! Distributes monthly total precipitation to days, given number of - ! monthly wet days. Adopted from LPX. - !-------------------------------------------------------------------- - use md_params_core, only: nmonth, ndayyear, ndaymonth - - ! arguments - real, dimension(nmonth), intent(in) :: mval_prec ! monthly precipitation totals - real, dimension(nmonth) :: mval_wet ! monthly number of wet days - real, dimension(ndayyear,2), intent(in) :: prdaily_random ! random number seed - - ! local vars - integer :: doy, moy, dm, iloop - integer :: daysum !accumulated days at beginning of month - integer :: nwet !# of generated wet days + ! function daily2monthly( dval, method ) result( mval ) + ! !///////////////////////////////////////////////////////////////////////// + ! ! Returns monthly values as a mean over daily values in each month. + ! !------------------------------------------------------------------------- + ! use md_params_core, only: ndaymonth, cumdaymonth, ndayyear, nmonth + + ! ! arguments + ! real, intent(in), dimension(ndayyear) :: dval ! vector containing 365 (366 in case lapyear is TRUE) daily values + ! character(len=*), intent(in) :: method ! true of monthly values represent total of daily values in resp. month + + ! ! function return variable + ! real, dimension(nmonth) :: mval + + ! ! local variables + ! integer :: moy + ! integer, dimension(nmonth) :: istart, iend + + ! istart = cumdaymonth - ndaymonth + 1 + ! iend = cumdaymonth + + ! ! loop over months and take sum/mean of days in that month + ! do moy=1,nmonth + ! if (method=="sum") then + ! mval(moy) = sum( dval( istart(moy) : iend(moy) )) + ! else if (method=="mean") then + ! mval(moy) = sum( dval( istart(moy) : iend(moy) )) / ndaymonth(moy) + ! else + ! stop 'DAILY2MONTHLY: select valid method (sum, mean)' + ! end if + ! end do + + ! end function daily2monthly + + + ! function monthly2daily_weather( mval_prec, mval_wet, prdaily_random ) result( dval_prec ) + ! !///////////////////////////////////////////////////////////////////////// + ! ! Weather (rain) generator. + ! ! Distributes monthly total precipitation to days, given number of + ! ! monthly wet days. Adopted from LPX. + ! !-------------------------------------------------------------------- + ! use md_params_core, only: nmonth, ndayyear, ndaymonth + + ! ! arguments + ! real, dimension(nmonth), intent(in) :: mval_prec ! monthly precipitation totals + ! real, dimension(nmonth) :: mval_wet ! monthly number of wet days + ! real, dimension(ndayyear,2), intent(in) :: prdaily_random ! random number seed + + ! ! local vars + ! integer :: doy, moy, dm, iloop + ! integer :: daysum !accumulated days at beginning of month + ! integer :: nwet !# of generated wet days - logical :: dry - - real :: prob, vv, v1, rand - real, parameter :: c1 = 1.0 - real, parameter :: c2 = 1.2 - real, dimension(nmonth) :: prob_rain - real, dimension(nmonth) :: mprecave !average precipitation on wet days - real, dimension(nmonth) :: mprecip !acc. monthly generated precipitation - - ! function return variable - real, dimension(ndayyear) :: dval_prec - - doy = 0 - prob = 0.0 - do moy=1,nmonth - prob_rain(moy) = 0.0 - mprecave(moy) = 0.0 - mprecip(moy) = 0.0 - enddo - daysum = 0 - - ! write(0,*) 'A' - - ! Initialize 2nd random number generator - ! call srand( int( prdaily_random(1,1) * 100000 ) ) - - ! xxx this solves the problem (removing multiplication with 1e5. But number of - ! actual wet days is not perfectly consistent with prescribed number of wet days. - call srand( int( prdaily_random(1,1) ) ) + ! logical :: dry + + ! real :: prob, vv, v1, rand + ! real, parameter :: c1 = 1.0 + ! real, parameter :: c2 = 1.2 + ! real, dimension(nmonth) :: prob_rain + ! real, dimension(nmonth) :: mprecave !average precipitation on wet days + ! real, dimension(nmonth) :: mprecip !acc. monthly generated precipitation + + ! ! function return variable + ! real, dimension(ndayyear) :: dval_prec + + ! doy = 0 + ! prob = 0.0 + ! do moy=1,nmonth + ! prob_rain(moy) = 0.0 + ! mprecave(moy) = 0.0 + ! mprecip(moy) = 0.0 + ! enddo + ! daysum = 0 + + ! ! write(0,*) 'A' + + ! ! Initialize 2nd random number generator + ! ! call srand( int( prdaily_random(1,1) * 100000 ) ) + + ! ! xxx this solves the problem (removing multiplication with 1e5. But number of + ! ! actual wet days is not perfectly consistent with prescribed number of wet days. + ! call srand( int( prdaily_random(1,1) ) ) - do moy=1,nmonth - if ( mval_wet(moy)<=1.0 ) mval_wet(moy) = 1.0 - prob_rain(moy) = mval_wet(moy) / real( ndaymonth(moy) ) - mprecave(moy) = mval_prec(moy) / mval_wet(moy) - dry = .TRUE. - iloop = 0 - - ! write(0,*) 'B' - - do while( dry ) - ! write(0,*) 'aa' - iloop = iloop + 1 - nwet = 0 - do dm=1,ndaymonth(moy) - ! write(0,*) 'a' - doy = doy + 1 + ! do moy=1,nmonth + ! if ( mval_wet(moy)<=1.0 ) mval_wet(moy) = 1.0 + ! prob_rain(moy) = mval_wet(moy) / real( ndaymonth(moy) ) + ! mprecave(moy) = mval_prec(moy) / mval_wet(moy) + ! dry = .TRUE. + ! iloop = 0 + + ! ! write(0,*) 'B' + + ! do while( dry ) + ! ! write(0,*) 'aa' + ! iloop = iloop + 1 + ! nwet = 0 + ! do dm=1,ndaymonth(moy) + ! ! write(0,*) 'a' + ! doy = doy + 1 - ! Transitional probabilities (Geng et al. 1986) - if (doy>1) then - if (dval_prec(doy-1) < 0.1) then - prob = 0.75 * prob_rain(moy) - else - prob = 0.25 + (0.75 * prob_rain(moy)) - endif - endif - ! write(0,*) 'b' + ! ! Transitional probabilities (Geng et al. 1986) + ! if (doy>1) then + ! if (dval_prec(doy-1) < 0.1) then + ! prob = 0.75 * prob_rain(moy) + ! else + ! prob = 0.25 + (0.75 * prob_rain(moy)) + ! endif + ! endif + ! ! write(0,*) 'b' - ! Determine we randomly and use Krysanova / Cramer estimates of - ! parameter values (c1,c2) for an exponential distribution - if (iloop==1) then - ! write(0,*) 'getting prdaily_random' - vv = real(prdaily_random(doy,1)) - ! write(0,*) 'vv', vv - else - ! write(0,*) 'calling rand' - ! xxx problem: rand() generates a random number that leads to floating point exception - vv = rand() - ! write(0,*) 'vv' - ! write(0,*) vv - endif - ! write(0,*) 'c' - ! write(0,*) 'prob', prob - if (vv>prob) then - ! write(0,*) 'd1' - dval_prec(doy) = 0.0 - else - ! write(0,*) 'c1' - nwet = nwet + 1 - ! write(0,*) 'c2' - v1 = real( prdaily_random(doy,2) ) - ! write(0,*) 'c3' - dval_prec(doy) = ((-log(v1)) ** c2) * mprecave(moy) * c1 - ! write(0,*) 'c4' - if (dval_prec(doy) < 0.1) dval_prec(doy) = 0.0 - ! write(0,*) 'c5' - endif - ! write(0,*) 'd' - mprecip(moy) = mprecip(moy) + dval_prec(doy) - enddo + ! ! Determine we randomly and use Krysanova / Cramer estimates of + ! ! parameter values (c1,c2) for an exponential distribution + ! if (iloop==1) then + ! ! write(0,*) 'getting prdaily_random' + ! vv = real(prdaily_random(doy,1)) + ! ! write(0,*) 'vv', vv + ! else + ! ! write(0,*) 'calling rand' + ! ! xxx problem: rand() generates a random number that leads to floating point exception + ! vv = rand() + ! ! write(0,*) 'vv' + ! ! write(0,*) vv + ! endif + ! ! write(0,*) 'c' + ! ! write(0,*) 'prob', prob + ! if (vv>prob) then + ! ! write(0,*) 'd1' + ! dval_prec(doy) = 0.0 + ! else + ! ! write(0,*) 'c1' + ! nwet = nwet + 1 + ! ! write(0,*) 'c2' + ! v1 = real( prdaily_random(doy,2) ) + ! ! write(0,*) 'c3' + ! dval_prec(doy) = ((-log(v1)) ** c2) * mprecave(moy) * c1 + ! ! write(0,*) 'c4' + ! if (dval_prec(doy) < 0.1) dval_prec(doy) = 0.0 + ! ! write(0,*) 'c5' + ! endif + ! ! write(0,*) 'd' + ! mprecip(moy) = mprecip(moy) + dval_prec(doy) + ! enddo - ! If it never rained this month and mprec(moy)>0 and mval_wet(moy)>0, do - ! again - dry = (nwet==0 .and. iloop<50 .and. mval_prec(moy)>0.1) - if (iloop>50) then - write(0,*) 'Daily.F, prdaily: Warning stopped after 50 tries in cell' - endif - - ! Reset counter to start of month - if (dry) then - doy = doy-ndaymonth(moy) - endif - - enddo !while + ! ! If it never rained this month and mprec(moy)>0 and mval_wet(moy)>0, do + ! ! again + ! dry = (nwet==0 .and. iloop<50 .and. mval_prec(moy)>0.1) + ! if (iloop>50) then + ! write(0,*) 'Daily.F, prdaily: Warning stopped after 50 tries in cell' + ! endif + + ! ! Reset counter to start of month + ! if (dry) then + ! doy = doy-ndaymonth(moy) + ! endif + + ! enddo !while - ! write(0,*) 'C' - - - ! normalise generated precipitation by monthly CRU values - if ( moy > 1 ) daysum = daysum + ndaymonth(moy-1) - if ( mprecip(moy) < 1.0 ) mprecip(moy) = 1.0 - do dm=1,ndaymonth(moy) - doy = daysum + dm - dval_prec(doy) = dval_prec(doy) * (mval_prec(moy) / mprecip(moy)) - if ( dval_prec(doy) < 0.1 ) dval_prec(doy) = 0.0 - ! dval_prec(doy) = mval_prec(moy) / ndaymonth(moy) !no generator - enddo + ! ! write(0,*) 'C' + + + ! ! normalise generated precipitation by monthly CRU values + ! if ( moy > 1 ) daysum = daysum + ndaymonth(moy-1) + ! if ( mprecip(moy) < 1.0 ) mprecip(moy) = 1.0 + ! do dm=1,ndaymonth(moy) + ! doy = daysum + dm + ! dval_prec(doy) = dval_prec(doy) * (mval_prec(moy) / mprecip(moy)) + ! if ( dval_prec(doy) < 0.1 ) dval_prec(doy) = 0.0 + ! ! dval_prec(doy) = mval_prec(moy) / ndaymonth(moy) !no generator + ! enddo - ! Alternative: equal distribution of rain for fixed number of wet days - ! prob = prob_rain(moy) + prob - ! if (prob.ge.1.0) then - ! dval_prec(doy) = mprec(moy) - ! prob = prob-1.0 - ! else - ! dval_prec(doy) = 0.0 - ! prob = prob - ! endif + ! ! Alternative: equal distribution of rain for fixed number of wet days + ! ! prob = prob_rain(moy) + prob + ! ! if (prob.ge.1.0) then + ! ! dval_prec(doy) = mprec(moy) + ! ! prob = prob-1.0 + ! ! else + ! ! dval_prec(doy) = 0.0 + ! ! prob = prob + ! ! endif - enddo !month + ! enddo !month - end function monthly2daily_weather + ! end function monthly2daily_weather - function monthly2daily( mval, method, monthistotal, mval_pvy, mval_nxy ) result( dval ) - !///////////////////////////////////////////////////////////////////////// - ! Returns daily values based on monthly values, using a defined method. - !------------------------------------------------------------------------- - use md_params_core, only: middaymonth, ndayyear, ndaymonth, nmonth + ! function monthly2daily( mval, method, monthistotal, mval_pvy, mval_nxy ) result( dval ) + ! !///////////////////////////////////////////////////////////////////////// + ! ! Returns daily values based on monthly values, using a defined method. + ! !------------------------------------------------------------------------- + ! use md_params_core, only: middaymonth, ndayyear, ndaymonth, nmonth - ! arguments - real, dimension(nmonth), intent(in) :: mval ! vector containing 12 monthly values - character(len=*), intent(in) :: method - logical, intent(in), optional :: monthistotal ! true of monthly values represent total of daily values in resp. month - real, dimension(nmonth), intent(in), optional :: mval_pvy ! vector containing 12 monthly values of the previous year - real, dimension(nmonth), intent(in), optional :: mval_nxy ! vector containing 12 monthly values of the next year - - ! function return variable - real, dimension(ndayyear) :: dval + ! ! arguments + ! real, dimension(nmonth), intent(in) :: mval ! vector containing 12 monthly values + ! character(len=*), intent(in) :: method + ! logical, intent(in), optional :: monthistotal ! true of monthly values represent total of daily values in resp. month + ! real, dimension(nmonth), intent(in), optional :: mval_pvy ! vector containing 12 monthly values of the previous year + ! real, dimension(nmonth), intent(in), optional :: mval_nxy ! vector containing 12 monthly values of the next year + + ! ! function return variable + ! real, dimension(ndayyear) :: dval - ! local variables - integer :: moy, doy, today, dm - real :: dd, todaysval - - real, dimension(0:(nmonth+1)) :: mval_ext - !integer, dimension(0:(nmonth+1)) :: middaymonth_ext - real :: startt, endt, starttemp, endtemp, dt, d2t, d3t, dtold, & - dtnew, lastmonthtemp, nextmonthtemp, deltatemp, polya, polyb, polyc + ! ! local variables + ! integer :: moy, doy, today, dm + ! real :: dd, todaysval + + ! real, dimension(0:(nmonth+1)) :: mval_ext + ! !integer, dimension(0:(nmonth+1)) :: middaymonth_ext + ! real :: startt, endt, starttemp, endtemp, dt, d2t, d3t, dtold, & + ! dtnew, lastmonthtemp, nextmonthtemp, deltatemp, polya, polyb, polyc - if (method == "interpol") then - !-------------------------------------------------------------------- - ! LINEAR INTERPOLATION - ! of monthly to quasi-daily values. - ! If optional argument 'mval_pvy' is provided, take December-value - ! of previous year to interpolate to first 15 days of January, - ! otherwise, use the same year's December value to get first 15 days. - ! corresponds to subroutine 'daily' in LPX - !-------------------------------------------------------------------- - - ! define extended vector with monthly values for previous Dec and next Jan added - mval_ext(1:nmonth) = mval(1:nmonth) - - !middaymonth_ext(1:nmonth) = middaymonth(1:nmonth) - !middaymonth_ext(0) = middaymonth(nmonth) - !middaymonth_ext(nmonth+1) = 381 - - if (present(mval_pvy)) then - mval_ext(0) = mval_pvy(nmonth) ! Dec value of previous year - else - mval_ext(0) = mval(nmonth) ! take Dec value of this year ==> leads to jump! - end if - - if (present(mval_nxy)) then - mval_ext(nmonth+1) = mval_nxy(1) ! Jan value of next year - else - mval_ext(nmonth+1) = mval(1) ! take Jan value of this year ==> leads to jump! - end if - - do moy = 1,nmonth - dd = (mval_ext(moy+1)-mval_ext(moy)) / real(middaymonth(moy+1) - middaymonth(moy)) - todaysval = mval_ext(moy) - do doy = middaymonth(moy),middaymonth(moy+1)-1 - if (doy<=ndayyear) then - today = doy - else - today = doy-ndayyear - endif - dval(today) = todaysval - todaysval = todaysval + dd - enddo - enddo - - if (monthistotal) then - doy = 0 - do moy=1,nmonth - do dm=1,ndaymonth(moy) - doy = doy+1 - dval(doy) = dval(doy) / real(ndaymonth(moy)) - enddo - enddo - endif - - !doy=1 - !do moy=1,nmonth - ! do dm=1,ndaymonth(moy) - ! doy=doy+1 - ! if (doy>middaymonth(moy)) then - ! ! interpolate to next month - ! dval(doy) = mval_ext(moy) + (doy-middaymonth_ext(moy))/ndaymonth_ext(moy) * (mval_ext(moy+1)-mval_ext(moy)) - ! else if (doy leads to jump! + ! end if + + ! if (present(mval_nxy)) then + ! mval_ext(nmonth+1) = mval_nxy(1) ! Jan value of next year + ! else + ! mval_ext(nmonth+1) = mval(1) ! take Jan value of this year ==> leads to jump! + ! end if + + ! do moy = 1,nmonth + ! dd = (mval_ext(moy+1)-mval_ext(moy)) / real(middaymonth(moy+1) - middaymonth(moy)) + ! todaysval = mval_ext(moy) + ! do doy = middaymonth(moy),middaymonth(moy+1)-1 + ! if (doy<=ndayyear) then + ! today = doy + ! else + ! today = doy-ndayyear + ! endif + ! dval(today) = todaysval + ! todaysval = todaysval + dd + ! enddo + ! enddo + + ! if (monthistotal) then + ! doy = 0 + ! do moy=1,nmonth + ! do dm=1,ndaymonth(moy) + ! doy = doy+1 + ! dval(doy) = dval(doy) / real(ndaymonth(moy)) + ! enddo + ! enddo + ! endif + + ! !doy=1 + ! !do moy=1,nmonth + ! ! do dm=1,ndaymonth(moy) + ! ! doy=doy+1 + ! ! if (doy>middaymonth(moy)) then + ! ! ! interpolate to next month + ! ! dval(doy) = mval_ext(moy) + (doy-middaymonth_ext(moy))/ndaymonth_ext(moy) * (mval_ext(moy+1)-mval_ext(moy)) + ! ! else if (doy1.0d-8) + ! read(20, 100, err=999) (tmp(l), l=1,3) + ! enddo - end function read1year_monthly + ! else + ! ! find corresponding year in first column and read 3 values on this line + ! read(20, 100, err=999) (tmp(l), l=1,3) + ! do while (abs(realyear-tmp(1))>1.0d-8) + ! read(20, 100, err=999) (tmp(l), l=1,3) + ! enddo + ! endif - function getvalreal( filename, realyear, day ) result( valreal ) - !//////////////////////////////////////////////////////////////// - ! Function reads one (annual) value corresponding to the given - ! year from a time series ascii file. - !---------------------------------------------------------------- - use md_params_core, only: ndayyear + ! valreal = tmp(2) - ! arguments - character(len=*), intent(in) :: filename - integer, intent(in) :: realyear - integer, intent(in), optional :: day ! day in year (1:365) - ! integer, intent(in), optional :: dm ! day in month (1:31) - ! integer, intent(in), optional :: mo ! month in year (1:12) - - ! function return value - real :: valreal + ! 100 format (30d16.8) + ! close(20) - ! local variables - integer :: l - real :: tmp(3) ! 3 so that an additional value for this year could be read - real :: realyear_decimal + ! return - if (present(day)) then - ! convert day number into decimal number - realyear_decimal = real(realyear) + real(day)/real(ndayyear) - endif + ! 888 write(0,*) 'GETVALREAL: error opening file '//trim(filename)//'. Abort. ' + ! stop + ! 999 write(0,*) 'GETVALREAL: error reading file '//trim(filename)//'. Abort. ' + ! stop - open(20, file='./input/'//filename, status='old', form='formatted', err=888) - - if (present(day)) then - ! find corresponding day in first column and read 3 values on this line - read(20, 100, err=999) (tmp(l), l=1,3) - do while (abs(realyear_decimal-tmp(1))>1.0d-8) - read(20, 100, err=999) (tmp(l), l=1,3) - enddo - - else - ! find corresponding year in first column and read 3 values on this line - read(20, 100, err=999) (tmp(l), l=1,3) - do while (abs(realyear-tmp(1))>1.0d-8) - read(20, 100, err=999) (tmp(l), l=1,3) - enddo - - endif - - valreal = tmp(2) - - 100 format (30d16.8) - close(20) - - return - - 888 write(0,*) 'GETVALREAL: error opening file '//trim(filename)//'. Abort. ' - stop - 999 write(0,*) 'GETVALREAL: error reading file '//trim(filename)//'. Abort. ' - stop - - end function getvalreal + ! end function getvalreal - function getvalreal_STANDARD( filename, realyear, mo, dm, day, realyear_decimal ) result( valreal ) - !//////////////////////////////////////////////////////////////// - ! Reads one (annual) value corresponding to the given year - ! from a time series ascii file. File has to be located in - ! ./input/ and has to contain only rows formatted like - ! '2002 1 1 0.496632 0.054053', which represents - ! 'YYYY MM DM GPP GPP err.'. DM is the day within the month. - ! If 'realyear' is dummy (-9999), then it's interpreted as to - ! represent a mean climatology for the course of a year. - ! XXX THIS IS NOT IN USE IN SOFUN ANYMORE XXX - !---------------------------------------------------------------- - ! arguments - character(len=*), intent(in) :: filename - integer, intent(in), optional :: realyear ! year AD as integer - integer, intent(in), optional :: mo ! month in year (1:12) - integer, intent(in), optional :: dm ! day in month (1:31 or 1:31 or 1:28) - integer, intent(in), optional :: day ! day in year (1:365) - real, intent(in), optional :: realyear_decimal ! year AD as decimal number corresponding to day in the year - - ! function return value - real :: valreal - - ! local variables - integer :: inyear - integer :: inmo - integer :: indm - integer :: inday - real :: inyear_decimal - real :: inval1 - real :: inval2 - - !print*,'looking for realyear, mo, dm',realyear,mo,dm - - ! open file - open(20, file='./input/'//filename, status='old', form='formatted', err=888) - - if (present(realyear)) then - ! DATA FOR EACH YEAR - if (present(mo)) then - ! DATA FOR EACH MONTH - if (present(dm)) then - ! DATA FOR EACH DAY IN THE MONTH - ! read the 2 values for this day in this year - read(20, 100, err=999) inyear, inmo, indm, inval1, inval2 - do while ( (realyear-inyear).ne.0 .or. (mo-inmo).ne.0 .or. (dm-indm).ne.0 ) - read(20, 100, err=999) inyear, inmo, indm, inval1, inval2 - enddo - else - ! read the 2 values for this month in this year - read(20, 200, err=999) inyear, inmo, inval1, inval2 - do while ( (realyear-inyear).ne.0 .or. (mo-inmo).ne.0 ) - read(20, 200, err=999) inyear, inmo, inval1, inval2 - enddo - end if - else if (present(day)) then - ! DATA FOR EACH DAY IN THE YEAR - ! read the 2 values for this day in this year - read(20, 700, err=999) inyear, inday, inval1, inval2 - do while ( (realyear-inyear).ne.0 .or. (day-inday).ne.0 ) - read(20, 700, err=999) inyear, inday, inval1, inval2 - enddo - else - ! read the 2 values for this year - read(20, 300, err=999) inyear, inval1, inval2 - do while ( (realyear-inyear).ne.0 ) - read(20, 300, err=999) inyear, inval1, inval2 - enddo - end if - else if (present(realyear_decimal)) then - ! DATA PROVIDED FOR EACH DAY AS A DECIMAL OF REALYEAR - ! find corresponding day in first column and read 3 values on this line - read(20, 900, err=999) inyear_decimal, inval1, inval2 - do while (abs(realyear_decimal-inyear_decimal).gt.1.0d-8) - read(20, 900, err=999) inyear_decimal, inval1, inval2 - enddo - else - ! DATA AS AVERAGE OVER MULTIPLE YEARS (recycle climatology) - ! FOR EACH MONTH, AND DAY-IN-THE-MONTH - if (present(mo)) then - if (present(dm)) then - ! read the 2 values for this day - read(20, 400, err=999) inmo, indm, inval1, inval2 - !print*,'inmo, indm, inval1, inval2', inmo, indm, inval1, inval2 - do while ( (mo-inmo).ne.0 .or. (dm-indm).ne.0 ) - read(20, 400, err=999) inmo, indm, inval1, inval2 - !print*,'inmo, indm, inval1, inval2', inmo, indm, inval1, inval2 - enddo - else - ! read the 2 values for this month - read(20, 500, err=999) inmo, inval1, inval2 - do while ( (mo-inmo).ne.0 ) - read(20, 500, err=999) inmo, inval1, inval2 - enddo - - end if - else if (present(day)) then - ! DATA FOR EACH DAY IN THE YEAR - ! read the 2 values for this day - read(20, 800, err=999) inday, inval1, inval2 - do while ( (day-inday).ne.0 ) - read(20, 800, err=999) inday, inval1, inval2 - enddo - else - ! read the 2 values in this input file - read(20, 600, err=999) inval1, inval2 - end if - endif - - !print*,'found realyear, mo, dm ',inyear,inmo,indm,inval1 - - valreal = inval1 - - 100 format (I4,I3,I3,F10.7,F10.7) - 200 format (I4,I3,F10.7,F10.7) - 300 format (I4,F10.7,F10.7) - 400 format (I3,I3,F10.7,F10.7) - 500 format (I3,F10.7,F10.7) - 600 format (F10.7,F10.7) - 700 format (I4,I4,F10.7,F10.7) - 800 format (I4,F10.7,F10.7) - 900 format (30d16.8,F10.7,F10.7) - - close(20) - - return - - 888 write(0,*) 'GETVALREAL_STANDARD: error opening file '//trim(filename)//'. Abort. ' - stop - 999 write(0,*) 'GETVALREAL_STANDARD: error reading file '//trim(filename)//'. Abort. ' - stop - - end function getvalreal_STANDARD - - - function getparreal( filename, paraname, try ) result( paravalue ) - !//////////////////////////////////////////////////////////////// - ! Low-level function for reading real parameter value from text file. - !---------------------------------------------------------------- - use md_params_core, only: dummy - - ! arguments - character(len=*), intent(in) :: filename, paraname - logical, intent(in), optional :: try - - ! function return value - real :: paravalue - - ! local variables - integer :: filehandle - character(len=40) :: readname, readvalue - - filehandle = 111 - open(filehandle,status='old',err=19,file=filename) - 9 read(filehandle,12,end=10)readname,readvalue - if (trim(readname)==paraname) then - read(readvalue,*) paravalue - goto 11 - else - goto 9 - endif - 10 continue - - if (present(try)) then - if (try) then - paravalue = dummy - end if - else - write(0,*) 'getparreal: '//paraname//' of type real not found' - stop - end if + ! function getvalreal_STANDARD( filename, realyear, mo, dm, day, realyear_decimal ) result( valreal ) + ! !//////////////////////////////////////////////////////////////// + ! ! Reads one (annual) value corresponding to the given year + ! ! from a time series ascii file. File has to be located in + ! ! ./input/ and has to contain only rows formatted like + ! ! '2002 1 1 0.496632 0.054053', which represents + ! ! 'YYYY MM DM GPP GPP err.'. DM is the day within the month. + ! ! If 'realyear' is dummy (-9999), then it's interpreted as to + ! ! represent a mean climatology for the course of a year. + ! ! XXX THIS IS NOT IN USE IN SOFUN ANYMORE XXX + ! !---------------------------------------------------------------- + ! ! arguments + ! character(len=*), intent(in) :: filename + ! integer, intent(in), optional :: realyear ! year AD as integer + ! integer, intent(in), optional :: mo ! month in year (1:12) + ! integer, intent(in), optional :: dm ! day in month (1:31 or 1:31 or 1:28) + ! integer, intent(in), optional :: day ! day in year (1:365) + ! real, intent(in), optional :: realyear_decimal ! year AD as decimal number corresponding to day in the year + + ! ! function return value + ! real :: valreal + + ! ! local variables + ! integer :: inyear + ! integer :: inmo + ! integer :: indm + ! integer :: inday + ! real :: inyear_decimal + ! real :: inval1 + ! real :: inval2 + + ! !print*,'looking for realyear, mo, dm',realyear,mo,dm + + ! ! open file + ! open(20, file='./input/'//filename, status='old', form='formatted', err=888) + + ! if (present(realyear)) then + ! ! DATA FOR EACH YEAR + ! if (present(mo)) then + ! ! DATA FOR EACH MONTH + ! if (present(dm)) then + ! ! DATA FOR EACH DAY IN THE MONTH + ! ! read the 2 values for this day in this year + ! read(20, 100, err=999) inyear, inmo, indm, inval1, inval2 + ! do while ( (realyear-inyear).ne.0 .or. (mo-inmo).ne.0 .or. (dm-indm).ne.0 ) + ! read(20, 100, err=999) inyear, inmo, indm, inval1, inval2 + ! enddo + ! else + ! ! read the 2 values for this month in this year + ! read(20, 200, err=999) inyear, inmo, inval1, inval2 + ! do while ( (realyear-inyear).ne.0 .or. (mo-inmo).ne.0 ) + ! read(20, 200, err=999) inyear, inmo, inval1, inval2 + ! enddo + ! end if + ! else if (present(day)) then + ! ! DATA FOR EACH DAY IN THE YEAR + ! ! read the 2 values for this day in this year + ! read(20, 700, err=999) inyear, inday, inval1, inval2 + ! do while ( (realyear-inyear).ne.0 .or. (day-inday).ne.0 ) + ! read(20, 700, err=999) inyear, inday, inval1, inval2 + ! enddo + ! else + ! ! read the 2 values for this year + ! read(20, 300, err=999) inyear, inval1, inval2 + ! do while ( (realyear-inyear).ne.0 ) + ! read(20, 300, err=999) inyear, inval1, inval2 + ! enddo + ! end if + ! else if (present(realyear_decimal)) then + ! ! DATA PROVIDED FOR EACH DAY AS A DECIMAL OF REALYEAR + ! ! find corresponding day in first column and read 3 values on this line + ! read(20, 900, err=999) inyear_decimal, inval1, inval2 + ! do while (abs(realyear_decimal-inyear_decimal).gt.1.0d-8) + ! read(20, 900, err=999) inyear_decimal, inval1, inval2 + ! enddo + ! else + ! ! DATA AS AVERAGE OVER MULTIPLE YEARS (recycle climatology) + ! ! FOR EACH MONTH, AND DAY-IN-THE-MONTH + ! if (present(mo)) then + ! if (present(dm)) then + ! ! read the 2 values for this day + ! read(20, 400, err=999) inmo, indm, inval1, inval2 + ! !print*,'inmo, indm, inval1, inval2', inmo, indm, inval1, inval2 + ! do while ( (mo-inmo).ne.0 .or. (dm-indm).ne.0 ) + ! read(20, 400, err=999) inmo, indm, inval1, inval2 + ! !print*,'inmo, indm, inval1, inval2', inmo, indm, inval1, inval2 + ! enddo + ! else + ! ! read the 2 values for this month + ! read(20, 500, err=999) inmo, inval1, inval2 + ! do while ( (mo-inmo).ne.0 ) + ! read(20, 500, err=999) inmo, inval1, inval2 + ! enddo + + ! end if + ! else if (present(day)) then + ! ! DATA FOR EACH DAY IN THE YEAR + ! ! read the 2 values for this day + ! read(20, 800, err=999) inday, inval1, inval2 + ! do while ( (day-inday).ne.0 ) + ! read(20, 800, err=999) inday, inval1, inval2 + ! enddo + ! else + ! ! read the 2 values in this input file + ! read(20, 600, err=999) inval1, inval2 + ! end if + ! endif + + ! !print*,'found realyear, mo, dm ',inyear,inmo,indm,inval1 + + ! valreal = inval1 + + ! 100 format (I4,I3,I3,F10.7,F10.7) + ! 200 format (I4,I3,F10.7,F10.7) + ! 300 format (I4,F10.7,F10.7) + ! 400 format (I3,I3,F10.7,F10.7) + ! 500 format (I3,F10.7,F10.7) + ! 600 format (F10.7,F10.7) + ! 700 format (I4,I4,F10.7,F10.7) + ! 800 format (I4,F10.7,F10.7) + ! 900 format (30d16.8,F10.7,F10.7) + + ! close(20) + + ! return + + ! 888 write(0,*) 'GETVALREAL_STANDARD: error opening file '//trim(filename)//'. Abort. ' + ! stop + ! 999 write(0,*) 'GETVALREAL_STANDARD: error reading file '//trim(filename)//'. Abort. ' + ! stop + + ! end function getvalreal_STANDARD + + + ! function getparreal( filename, paraname, try ) result( paravalue ) + ! !//////////////////////////////////////////////////////////////// + ! ! Low-level function for reading real parameter value from text file. + ! !---------------------------------------------------------------- + ! use md_params_core, only: dummy + + ! ! arguments + ! character(len=*), intent(in) :: filename, paraname + ! logical, intent(in), optional :: try + + ! ! function return value + ! real :: paravalue + + ! ! local variables + ! integer :: filehandle + ! character(len=40) :: readname, readvalue + + ! filehandle = 111 + ! open(filehandle,status='old',err=19,file=filename) + ! 9 read(filehandle,12,end=10)readname,readvalue + ! if (trim(readname)==paraname) then + ! read(readvalue,*) paravalue + ! goto 11 + ! else + ! goto 9 + ! endif + ! 10 continue + + ! if (present(try)) then + ! if (try) then + ! paravalue = dummy + ! end if + ! else + ! write(0,*) 'getparreal: '//paraname//' of type real not found' + ! stop + ! end if - 11 continue - 12 format(2a40) + ! 11 continue + ! 12 format(2a40) - close(filehandle) + ! close(filehandle) - ! print*,'reading real value for ', paraname, ': ', paravalue + ! ! print*,'reading real value for ', paraname, ': ', paravalue - return + ! return - 19 write(0,*) 'getparreal: '//filename//' not found!' - stop + ! 19 write(0,*) 'getparreal: '//filename//' not found!' + ! stop - end function getparreal + ! end function getparreal - function getparint( filename, paraname ) result( paravalue ) - !//////////////////////////////////////////////////////////////// - ! Low-level function for reading integer parameter value from text file - !---------------------------------------------------------------- - ! arguments - character(len=*), intent(in) :: filename, paraname + ! function getparint( filename, paraname ) result( paravalue ) + ! !//////////////////////////////////////////////////////////////// + ! ! Low-level function for reading integer parameter value from text file + ! !---------------------------------------------------------------- + ! ! arguments + ! character(len=*), intent(in) :: filename, paraname - ! function return variable - integer :: paravalue - - ! local variables - integer :: filehandle - character(len=40) :: readname, readvalue - - filehandle = 111 - open(filehandle,status='old',err=19,file=filename) - 9 read(filehandle,12,end=10)readname,readvalue - if (trim(readname)==paraname) then - read(readvalue,*) paravalue - goto 11 - else - goto 9 - endif - 10 continue - write(0,*) 'getparint: in file '//filename//':' - write(0,*) 'getparint: '//paraname//' of type integer not found' - stop + ! ! function return variable + ! integer :: paravalue + + ! ! local variables + ! integer :: filehandle + ! character(len=40) :: readname, readvalue + + ! filehandle = 111 + ! open(filehandle,status='old',err=19,file=filename) + ! 9 read(filehandle,12,end=10)readname,readvalue + ! if (trim(readname)==paraname) then + ! read(readvalue,*) paravalue + ! goto 11 + ! else + ! goto 9 + ! endif + ! 10 continue + ! write(0,*) 'getparint: in file '//filename//':' + ! write(0,*) 'getparint: '//paraname//' of type integer not found' + ! stop + + ! 11 continue + ! 12 format(2a40) + + ! close(filehandle) + + ! ! print*,'reading integer for ', paraname, ': ', paravalue + + ! return + + ! 19 write(0,*) 'getparint: file '//filename//' not found:' + ! write(0,*) filename + ! stop + + ! end function getparint + + + ! function getparlogical( filename, paraname ) result( paravalue ) + ! !//////////////////////////////////////////////////////////////// + ! ! Low-level function for reading boolean parameter value from text file + ! !---------------------------------------------------------------- + ! ! arguments + ! character(len=*), intent(in) :: filename, paraname + ! ! function return variable + ! logical :: paravalue + + ! ! local variables + ! integer :: filehandle + ! character(len=40) :: readname, readvalue + + ! filehandle = 111 + ! open(filehandle,status='old',err=19,file=filename) + ! 9 read(filehandle,12,end=10)readname,readvalue + ! if (trim(readname)==paraname) then + ! read(readvalue,*) paravalue + ! goto 11 + ! else + ! goto 9 + ! endif + ! 10 continue + ! write(0,*) 'GETPARLOGICAL: '//paraname//' of type logical not found' + ! stop + + ! 11 continue + ! 12 format(2a40) + + ! close(filehandle) + + ! ! print*,'reading logical for ', paraname, ': ', paravalue + + ! return + + ! 19 write(0,*) 'GETPARLOGICAL: '//filename//' not found!' + ! stop + + ! end function getparlogical + + + ! subroutine getparstring( filename, paraname, paravalue ) + ! !//////////////////////////////////////////////////////////////// + ! ! Low-level function for reading parameter value (string) from text file + ! !---------------------------------------------------------------- + ! character(len=*), intent(in) :: filename, paraname + ! character(len=*) :: paravalue + + ! integer :: filehandle,i + ! character(len=40) :: readname + ! character(len=1024) :: readvalue + + ! filehandle = 111 + ! open(filehandle,status='old',err=19,file=filename) + ! 9 read(filehandle,12,end=10)readname,readvalue + ! if (trim(readname)==paraname) then + ! ! Strip leading and trailing whitespace from readvalue + ! i=1 + ! do while (readvalue(i:i)==' ' .and. i.lt.len(readvalue)) + ! i=i+1 + ! enddo + ! paravalue = trim(readvalue(i:)) + ! goto 11 + ! else + ! goto 9 + ! endif + ! 10 continue + ! write(0,*) 'GETSTRING: '//paraname//' of type string not found' + ! stop + + ! 11 continue + ! 12 format(a40,a) + + ! close(filehandle) + + ! print*, 'reading string for ', paraname, ': ', paravalue - 11 continue - 12 format(2a40) + ! return - close(filehandle) - - ! print*,'reading integer for ', paraname, ': ', paravalue - - return - - 19 write(0,*) 'getparint: file '//filename//' not found:' - write(0,*) filename - stop - - end function getparint - - - function getparlogical( filename, paraname ) result( paravalue ) - !//////////////////////////////////////////////////////////////// - ! Low-level function for reading boolean parameter value from text file - !---------------------------------------------------------------- - ! arguments - character(len=*), intent(in) :: filename, paraname - ! function return variable - logical :: paravalue - - ! local variables - integer :: filehandle - character(len=40) :: readname, readvalue - - filehandle = 111 - open(filehandle,status='old',err=19,file=filename) - 9 read(filehandle,12,end=10)readname,readvalue - if (trim(readname)==paraname) then - read(readvalue,*) paravalue - goto 11 - else - goto 9 - endif - 10 continue - write(0,*) 'GETPARLOGICAL: '//paraname//' of type logical not found' - stop - - 11 continue - 12 format(2a40) - - close(filehandle) - - ! print*,'reading logical for ', paraname, ': ', paravalue - - return - - 19 write(0,*) 'GETPARLOGICAL: '//filename//' not found!' - stop - - end function getparlogical - - - subroutine getparstring( filename, paraname, paravalue ) - !//////////////////////////////////////////////////////////////// - ! Low-level function for reading parameter value (string) from text file - !---------------------------------------------------------------- - character(len=*), intent(in) :: filename, paraname - character(len=*) :: paravalue - - integer :: filehandle,i - character(len=40) :: readname - character(len=1024) :: readvalue - - filehandle = 111 - open(filehandle,status='old',err=19,file=filename) - 9 read(filehandle,12,end=10)readname,readvalue - if (trim(readname)==paraname) then - ! Strip leading and trailing whitespace from readvalue - i=1 - do while (readvalue(i:i)==' ' .and. i.lt.len(readvalue)) - i=i+1 - enddo - paravalue = trim(readvalue(i:)) - goto 11 - else - goto 9 - endif - 10 continue - write(0,*) 'GETSTRING: '//paraname//' of type string not found' - stop - - 11 continue - 12 format(a40,a) - - close(filehandle) - - print*, 'reading string for ', paraname, ': ', paravalue - - return - - 19 write(0,*) 'GETSTRING: '//filename//' not found!' - stop - - end subroutine getparstring - - - function getparchar( filename, paraname ) result( paravalue ) - !//////////////////////////////////////////////////////////////// - ! Low-level function for reading parameter value from text file - !---------------------------------------------------------------- - character(len=*), intent(in) :: filename, paraname - character(len=1) :: paravalue - - integer :: filehandle - character(len=40) :: readname, readvalue - - filehandle = 111 - open(filehandle,status='old',err=19,file=filename) - 9 read(filehandle,12,end=10)readname,readvalue - if (trim(readname)==paraname) then - read(readvalue,*) paravalue - goto 11 - else - goto 9 - endif - 10 continue - write(0,*) 'getparchar: '//paraname//' of type char not found' - stop + ! 19 write(0,*) 'GETSTRING: '//filename//' not found!' + ! stop - 11 continue - 12 format(2a40) + ! end subroutine getparstring - close(filehandle) - return + ! function getparchar( filename, paraname ) result( paravalue ) + ! !//////////////////////////////////////////////////////////////// + ! ! Low-level function for reading parameter value from text file + ! !---------------------------------------------------------------- + ! character(len=*), intent(in) :: filename, paraname + ! character(len=1) :: paravalue + + ! integer :: filehandle + ! character(len=40) :: readname, readvalue + + ! filehandle = 111 + ! open(filehandle,status='old',err=19,file=filename) + ! 9 read(filehandle,12,end=10)readname,readvalue + ! if (trim(readname)==paraname) then + ! read(readvalue,*) paravalue + ! goto 11 + ! else + ! goto 9 + ! endif + ! 10 continue + ! write(0,*) 'getparchar: '//paraname//' of type char not found' + ! stop + + ! 11 continue + ! 12 format(2a40) + + ! close(filehandle) - 19 write(0,*) 'getparchar: '//filename//' not found!' - stop + ! return + + ! 19 write(0,*) 'getparchar: '//filename//' not found!' + ! stop - end function getparchar + ! end function getparchar function area( lat, dx, dy ) result( out_area ) diff --git a/src/soiltemp_sitch.mod.f90 b/src/soiltemp_sitch.mod.f90 index bd5ad3ba..13960782 100644 --- a/src/soiltemp_sitch.mod.f90 +++ b/src/soiltemp_sitch.mod.f90 @@ -44,7 +44,8 @@ subroutine soiltemp( soil, dtemp, ngridcells, jpngr, moy, doy ) ! Calculates soil temperature based on. !------------------------------------------------------------------------- use md_params_core, only: ndayyear, nlu, ndaymonth, pi - use md_sofunutils, only: running, daily2monthly + use md_sofunutils, only: running + ! use md_sofunutils, only: daily2monthly use md_tile_pmodel, only: soil_type use md_interface_pmodel, only: myinterface diff --git a/src/waterbal_splash.mod.f90 b/src/waterbal_splash.mod.f90 index c8e05eff..ac26c1a1 100644 --- a/src/waterbal_splash.mod.f90 +++ b/src/waterbal_splash.mod.f90 @@ -27,7 +27,7 @@ module md_waterbal use md_forcing_pmodel, only: climate_type use md_grid, only: gridtype use md_interface_pmodel, only: myinterface - use md_sofunutils, only: daily2monthly, radians, dgsin, dgcos, degrees + use md_sofunutils, only: radians, dgsin, dgcos, degrees implicit none