From 8b0a9ea3601c9d47f1f39aba70313af54b85eb0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20Schr=C3=B6n?= Date: Mon, 30 Nov 2015 23:05:18 +0000 Subject: [PATCH] B:trunk: merged with B:neutron_calib --- mhm.nml | 101 +++++---- pre-proc/enlargegrid.pl | 1 + src/mHM/mhm_driver.f90 | 48 ++-- src/mHM/mo_global_variables.f90 | 142 ++++++------ src/mHM/mo_init_states.f90 | 140 ++++++------ src/mHM/mo_mhm.f90 | 109 +++++---- src/mHM/mo_mhm_constants.f90 | 28 +-- src/mHM/mo_mhm_eval.f90 | 193 +++++++++------- src/mHM/mo_objective_function.f90 | 346 ++++++++++++++++++++--------- src/mHM/mo_read_config.f90 | 37 +-- src/mHM/mo_read_optional_data.f90 | 147 ++++++++++-- src/mHM/mo_write_ascii.f90 | 89 ++++---- src/mHM/mo_write_fluxes_states.f90 | 130 +++++------ 13 files changed, 912 insertions(+), 599 deletions(-) diff --git a/mhm.nml b/mhm.nml index fcab8fe0..e7d4722c 100644 --- a/mhm.nml +++ b/mhm.nml @@ -22,11 +22,11 @@ ! Matthias Zink, Nov 2014 - added multiple options for process 5 - PET ! Matthias Zink, Dec 2014 - adopted inflow gauges to ignore headwater cells ! Matthias Zink, Mar 2015 - added optional soil mositure read in for calibration - + !****************************************************************************************** ! !****************************************************************************************** -! MAIN +! MAIN !****************************************************************************************** !> Main namelist !> Most of the variables (if not all) given in this namelist are common @@ -37,13 +37,13 @@ !----------------------------------------------------------------------------- timestep = 1 !----------------------------------------------------------------------------- -!> input data & model run cordinate system +!> input data & model run cordinate system !> 0 -> regular X & Y coordinate system (e.g., GK-4 or Lambert equal area system) !> 1 -> regular lat & lon coordinate system !----------------------------------------------------------------------------- iFlag_cordinate_sys = 0 !----------------------------------------------------------------------------- -!> Number of basins to be modeled. \n +!> Number of basins to be modeled. \n !> Number given here should correspond to one given in "gaugeinfo.txt" file.\n !> All gauging stations within those basins will be taken for the optimization.\n !> IF routing process is ON then give nBasins = 1, for this case, mHM will internally @@ -52,19 +52,19 @@ iFlag_cordinate_sys = 0 nBasins = 2 !----------------------------------------------------------------------------- !> resolution of Level-1 hydrological simulations in mHM [m or degree] per basin -!> NOTE: if iFlag_cordinate_sys = 0, then resolution_Hydrology is in [m] \n -!> if iFlag_cordinate_sys = 1, then resolution_Hydrology is in [degree-decimal] \n +!> NOTE: if iFlag_cordinate_sys = 0, then resolution_Hydrology is in [m] \n +!> if iFlag_cordinate_sys = 1, then resolution_Hydrology is in [degree-decimal] \n !----------------------------------------------------------------------------- -resolution_Hydrology(1) = 12000 -resolution_Hydrology(2) = 24000 +resolution_Hydrology(1) = 12000 +resolution_Hydrology(2) = 24000 !----------------------------------------------------------------------------- !> resolution of Level-11 discharge routing [m or degree] per basin \n !> this level-11 discharge routing resolution must be >= and multiple of the !> level-1 hydrological simulations resolution \n -!> NOTE: if iFlag_cordinate_sys = 0, then resolution_Routing is in [m] \n -!> if iFlag_cordinate_sys = 1, then resolution_Routing is in [degree-decimal] \n +!> NOTE: if iFlag_cordinate_sys = 0, then resolution_Routing is in [m] \n +!> if iFlag_cordinate_sys = 1, then resolution_Routing is in [degree-decimal] \n !----------------------------------------------------------------------------- -resolution_Routing(1) = 12000 +resolution_Routing(1) = 12000 resolution_Routing(2) = 24000 !---------------------------------------------------------------------------- !> specify same index for basins to share L0_data to save memory \n @@ -75,7 +75,7 @@ L0Basin(1) = 1 L0Basin(2) = 2 !----------------------------------------------------------------------------- !> flag for optimization: .TRUE.: optimization -!> or .FALSE.: no optimazition +!> or .FALSE.: no optimazition !----------------------------------------------------------------------------- optimize = .FALSE. !> Optimization shall be restarted from ./mo_.restart file, which @@ -103,6 +103,7 @@ opti_method = 1 !> (14) SO: Q: sum[((1.0-KGE_i)/ nGauges)**6]**(1/6) > combination of KGE of every gauging station based on a power-6 norm \n !> (15) SO: Q + basin_avg_TWS: [1.0-KGE(Q)]*RMSE(basin_avg_TWS) - objective function using Q and basin average (standard score) TWS\n !> (16) (reserved) please use the next number when implementing a new one +!> (17) SO: N: 1.0 - KGE of spatio-temporal neutron data, catchment-average !> further functions can be implemented in mo_objective_function opti_function = 3 !----------------------------------------------------------------------------- @@ -126,12 +127,12 @@ mrm_coupling_mode = 2 / !****************************************************************************************** -! DIRECTORIES +! DIRECTORIES !****************************************************************************************** -!> Namelist with all directories for common file as well as separate file for every basin.\n +!> Namelist with all directories for common file as well as separate file for every basin.\n !> Number in brackets indicates basin number.\n !> This number HAS TO correspond with the number of basin given below in the "mainconfig" -!> namelist as well as the indices given in "evaluation_gauges" namelist.\n +!> namelist as well as the indices given in "evaluation_gauges" namelist.\n ! directories used only by mRM &directories_mRM @@ -152,7 +153,7 @@ dir_Total_Runoff(2) = 'test_basin_2/output/' &directories_general ! !----------------------------------------------------- -! Directory for common files (to all basins) +! Directory for common files (to all basins) !----------------------------------------------------- !> config run out file common to all modeled basins should be written to directory dirConfigOut = "test_basin/" @@ -206,7 +207,7 @@ inputFormat_meteo_forcings = "nc" dir_Precipitation(1) = "test_basin/input/meteo/pre/" dir_Temperature(1) = "test_basin/input/meteo/tavg/" !> paths depending on PET process (processCase(5)) -!> 0 - PET is input +!> 0 - PET is input !> 1 - Hargreaves-Sammani method !> 2 - Priestley-Taylor mehtod !> 3 - Penman-Monteith method @@ -228,7 +229,7 @@ dir_gridded_LAI(1) = "test_basin/input/lai/" dir_Precipitation(2) = "test_basin_2/input/meteo/pre/" dir_Temperature(2) = "test_basin_2/input/meteo/tavg/" !> paths depending on PET process (processCase(5)) -!> 0 - PET is input +!> 0 - PET is input !> 1 - Hargreaves-Sammani method !> 2 - Priestley-Taylor mehtod !> 3 - Penman-Monteith method @@ -270,22 +271,26 @@ timeStep_sm_input = -2 !> file name including path with timeseries of GRACE-based data file_tws(1) = "test_basin/input/optional_data/tws_basin_1.txt" file_tws(2) = "test_basin/input/optional_data/tws_basin_2.txt" +! +!> directory to neutron data +! expected file name: neutrons.nc, expected variable name: neutrons +dir_neutrons(1) = "test_basin/input/optional_data/" / &time_periods !----------------------------------------------------------------------------- !> specification of number of warming days [d] and the simulation period.\n -!> All dynamic data sets(e.g., meteo. forcings, landcover scenes) should start +!> All dynamic data sets(e.g., meteo. forcings, landcover scenes) should start !> from warming days and ends at the last day of the evaluation period. \n ! !> 1---------2-------------------3 !> !> 1-> Starting of the effective modeling period (including the warming days) !> 2-> Starting of the given simulation period -!> 3-> Ending of the given simulation period (= end of the effective modeling period) +!> 3-> Ending of the given simulation period (= end of the effective modeling period) ! !> IF you want to run the model from 2002/01/01 (Starting of the given simulation -!> period=2) to 2003/12/31 (End of the given simulation period=3) with 365 warming +!> period=2) to 2003/12/31 (End of the given simulation period=3) with 365 warming !> day, which is 2001/01/01 = 1), THEN all dynamic datasets should be given for !> the effective modeling period of 2001/01/01 to 2003/12/31. !----------------------------------------------------------------------------- @@ -309,7 +314,7 @@ eval_Per(2)%mEnd = 12 !> last day of wanted simulation period eval_Per(1)%dEnd = 31 eval_Per(2)%dEnd = 31 -!> switch to control read input frequency of the gridded meteo input, +!> switch to control read input frequency of the gridded meteo input, !> i.e. precipitation, potential evapotransiration, and temperature\n !> >0: after each days\n !> 0: only at beginning of the run\n @@ -323,12 +328,12 @@ time_step_model_inputs(2) = 0 / !****************************************************************************************** -! mHM SOIL LAYERING +! mHM SOIL LAYERING !****************************************************************************************** !> Namelist controlling the layering of the soil !> Variables given in this namelist are common to all basins to be modeled. &soilLayer -!> [mm] soil depth down to which organic matter is possible +!> [mm] soil depth down to which organic matter is possible tillageDepth = 200 !> No. of soil horizons to be modeled nSoilHorizons_mHM = 2 @@ -338,16 +343,16 @@ soil_Depth(1) = 200 / !****************************************************************************************** -! LAND COVER +! LAND COVER !****************************************************************************************** &LCover !> Variables given in this namelist are common to all basins to be modeled. -!> Please make sure that the land cover periods are covering the simulation period. +!> Please make sure that the land cover periods are covering the simulation period. !>fraction of area within city assumed to be fully sealed [0.0-1.0] fracSealed_cityArea = 0.6 !> number of land cover scenes to be used\n !> The land cover scene periods are shared by all catchments. -!> The names should be equal for all basins. The land cover scnes have to be ordered +!> The names should be equal for all basins. The land cover scnes have to be ordered !> chronologically. nLcover_scene = 3 ! indicate period with brackets behind variable @@ -380,16 +385,16 @@ LCoverfName(3) = 'lc_1990.asc' &LAI_data_information ! !----------------------------------------------------------------------------------- -!> Flag timeStep_LAI_input identifies how LAI is read in mHM. +!> Flag timeStep_LAI_input identifies how LAI is read in mHM. !> This flag is unique and valid for all basins.\n ! !> timeStep_LAI_input !> 0: read LAI from longterm monthly mean lookup table (related to land cover file). !> The filename (LAI_classdefinition.txt) for the LUT is hard coded in mo_file.f90 !> Information regarding long-term monthly mean LAI for land cover classes -!> appearing in all modeled basins should be included in this LUT file. +!> appearing in all modeled basins should be included in this LUT file. !> This is an unique file applicable to all basins to be modeled. \n -!> The respective plant functional type is in LAI_class.asc, which must be +!> The respective plant functional type is in LAI_class.asc, which must be !> located in each basin's morph directory. !> <0: Read gridded LAI files.\n !> -1: gridded LAI are daily values\n @@ -408,13 +413,13 @@ inputFormat_gridded_LAI = "nc" ! PROCESSES !****************************************************************************************** !> This matrix manages which processes and process descriptions are used for simulation. -!> The number of processes and its corresponding numbering are fixed. The process description can be +!> The number of processes and its corresponding numbering are fixed. The process description can be !> chosen from the options listed above the name of the particular process case. This number has to be !> given for processCase(*). ! &processSelection !> interception -!> 1 - maximum Interception +!> 1 - maximum Interception processCase(1) = 1 !> snow !> 1 - degree-day approach @@ -432,7 +437,7 @@ processCase(4) = 1 !> 3 - Penman-Monteith processCase(5) = 0 !> interflow -!> 1 - storage reservoir with one outflow threshold and nonlinear response +!> 1 - storage reservoir with one outflow threshold and nonlinear response processCase(6) = 1 !> percolation !> 1 - GW assumed as linear reservoir @@ -447,7 +452,7 @@ processCase(9) = 1 !> ground albedo of cosmic-ray neutrons !> THIS IS WORK IN PROGRESS, DO NOT USE FOR RESEARCH !> 0 - deactivated -!> 1 - inverse N0 based on Desilets et al. 2010 +!> 1 - inverse N0 based on Desilets et al. 2010 !> 2 - COSMIC forward operator by Shuttlworth et al. 2013 processCase(10) = 0 / @@ -457,7 +462,7 @@ processCase(10) = 0 !****************************************************************************************** !> namelist controlling the gauging station information !> The ID has to correspond to the ID's given in the 'gaugelocation.asc' and -!> to the filename containing the time series +!> to the filename containing the time series &evaluation_gauges !> Gauges for model evaluation ! @@ -471,26 +476,26 @@ nGaugesTotal = 2 !> basin 1 !> number of gauges for subbasin (1) NoGauges_basin(1) = 1 -!> id of gauge(1) for subbasin(1) --> (1,1) +!> id of gauge(1) for subbasin(1) --> (1,1) Gauge_id(1,1) = 398 -!> name of file with timeseries of gauge(1) for subbasin(1) --> (1,1) +!> name of file with timeseries of gauge(1) for subbasin(1) --> (1,1) gauge_filename(1,1) = "00398.txt" ! !> basin 2 !> number of gauges for subbasin (2) NoGauges_basin(2) = 1 -!> id of gauge(1) for subbasin(2) --> (2,1) +!> id of gauge(1) for subbasin(2) --> (2,1) Gauge_id(2,1) = 45 -!> name of file with timeseries of gauge(1) for subbasin(2) --> (2,1) +!> name of file with timeseries of gauge(1) for subbasin(2) --> (2,1) Gauge_filename(2,1) = "45.txt" / &inflow_gauges !> Gauges / gridpoints used for inflow to the model domain -!> e.g. in the case of upstream/headwater areas which are +!> e.g. in the case of upstream/headwater areas which are !> not included in the model domain ! -!> Total number of inflow gauges (sum of all gauges in all subbains) +!> Total number of inflow gauges (sum of all gauges in all subbains) nInflowGaugesTotal = 0 !> structure of gauge_id(i,j) & gauge_filename(i,j): !> 1st dimension is the number of the subbasin i @@ -500,11 +505,11 @@ nInflowGaugesTotal = 0 !> basin 1 !> number of gauges for subbasin (1) NoInflowGauges_basin(1) = 0 -!> id of inflow gauge(1) for subbasin(1) --> (1,1) +!> id of inflow gauge(1) for subbasin(1) --> (1,1) InflowGauge_id(1,1) = -9 -!> name of file with timeseries of inflow gauge(1) for subbasin(1) --> (1,1) +!> name of file with timeseries of inflow gauge(1) for subbasin(1) --> (1,1) InflowGauge_filename(1,1) = "" -!> consider flows from upstream/headwater cells of inflow gauge(1) for subbasin(1) --> (1,1) +!> consider flows from upstream/headwater cells of inflow gauge(1) for subbasin(1) --> (1,1) InflowGauge_Headwater(1,1) = .FALSE. / @@ -518,13 +523,13 @@ evap_coeff = 1.30, 1.20, 0.72, 0.75, 1.00, 1.00, 1.00, 1.00, 1.00, / !****************************************************************************************** -! ANNUAL CYCLE METEOROLOGICAL FORCINGS +! ANNUAL CYCLE METEOROLOGICAL FORCINGS !****************************************************************************************** &nightDayRatio !> night ratio for precipitation !> only night values required because day values are the opposite fnight_prec = 0.46, 0.50, 0.52, 0.51, 0.48, 0.50, 0.49, 0.48, 0.52, 0.56, 0.50, 0.47 -!> night ratio for PET +!> night ratio for PET fnight_pet = 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10 !> night correction factor for temperature fnight_temp = -0.76, -1.30, -1.88, -2.38, -2.72, -2.75, -2.74, -3.04, -2.44, -1.60, -0.94, -0.53 @@ -568,8 +573,8 @@ sce_nps = -9 ! ------------------------------------- !> MCMC specific: \n ! ------------------------------------- -!> .true.: use MCMC for optimisation and estimation of parameter uncertainty -!> .false.: use MCMC for estimation of parameter uncertainty +!> .true.: use MCMC for optimisation and estimation of parameter uncertainty +!> .false.: use MCMC for estimation of parameter uncertainty mcmc_opti = .false. !> Parameters of error model if mcmc_opti=.false. !> e.g. for opti_function=8: two parameters a and b: err = a + b*Q diff --git a/pre-proc/enlargegrid.pl b/pre-proc/enlargegrid.pl index fd729a13..a1bb03c8 100755 --- a/pre-proc/enlargegrid.pl +++ b/pre-proc/enlargegrid.pl @@ -20,6 +20,7 @@ { local $^I = ".bak"; # backup files local @ARGV = glob("*.asc"); + print "$_\n" foreach @ARGV; while (<>) { if (/^ncols/) { $datastart=0;$xdiff=-1;$ydiff=-1; print STDOUT "$ARGV\n"; } if ($datastart) { s/[\r\n]+//; print $_.(" $NODATA" x $xdiff)."\n"; } diff --git a/src/mHM/mhm_driver.f90 b/src/mHM/mhm_driver.f90 index 16b08508..5766e268 100644 --- a/src/mHM/mhm_driver.f90 +++ b/src/mHM/mhm_driver.f90 @@ -140,7 +140,7 @@ ! Stephan Thober, Aug 2015 - removed routing related variables ! Stephan Thober, Oct 2015 - reorganized optimization (now compatible with mRM) ! Oldrich Rakovec, Rohini Kumar, Oct 2015 - added reading of basin averaged TWS and objective function 15 -! for simultaneous calibration based on runoff and TWS +! for simultaneous calibration based on runoff and TWS ! ! -------------------------------------------------------------------------- @@ -179,7 +179,9 @@ PROGRAM mhm_driver USE mo_meteo_forcings, ONLY : prepare_meteo_forcings_data USE mo_mhm_eval, ONLY : mhm_eval USE mo_prepare_gridded_LAI, ONLY : prepare_gridded_daily_LAI_data ! prepare daily LAI gridded fields - USE mo_read_optional_data, ONLY : read_soil_moisture, read_basin_avg_TWS ! optional soil moisture reader, basin_avg_TWS reader + USE mo_read_optional_data, ONLY : read_soil_moisture, & ! optional soil moisture reader, basin_avg_TWS reader + read_basin_avg_TWS, & + read_neutrons USE mo_read_config, ONLY : read_config ! Read main configuration files USE mo_read_wrapper, ONLY : read_data ! Read all input data USE mo_read_latlon, ONLY : read_latlon @@ -192,13 +194,13 @@ PROGRAM mhm_driver write_configfile, & ! Writing Configuration file write_optifile, & ! Writing optimized parameter set and objective write_optinamelist ! Writing optimized parameter set to a namelist - USE mo_objective_function, ONLY : objective ! objective functions and likelihoods + USE mo_objective_function, ONLY : objective ! objective functions and likelihoods USE mo_optimization, ONLY : optimization #ifdef mrm2mhm USE mo_mrm_objective_function_runoff, only: single_objective_runoff USE mo_mrm_init, ONLY : mrm_init USE mo_mrm_write, only : mrm_write -#endif +#endif !$ USE omp_lib, ONLY : OMP_GET_NUM_THREADS ! OpenMP routines IMPLICIT NONE @@ -212,7 +214,7 @@ PROGRAM mhm_driver real(dp) :: funcbest ! best objective function achivied during optimization logical, dimension(:), allocatable :: maskpara ! true = parameter will be optimized = parameter(i,4) = 1 ! ! false = parameter will not be optimized = parameter(i,4) = 0 - + ! -------------------------------------------------------------------------- ! START ! -------------------------------------------------------------------------- @@ -224,7 +226,7 @@ PROGRAM mhm_driver call message(' ', trim(version_date)) call message() call message('Originally by L. Samaniego & R. Kumar') - + call message(separator) call message() @@ -264,16 +266,16 @@ PROGRAM mhm_driver call message(' Temperature directory: ', trim(dirTemperature(ii) )) select case (processMatrix(5,1)) case(0) - call message(' PET directory: ', trim(dirReferenceET(ii) )) + call message(' PET directory: ', trim(dirReferenceET(ii) )) case(1) - call message(' Min. temperature directory: ', trim(dirMinTemperature(ii) )) - call message(' Max. temperature directory: ', trim(dirMaxTemperature(ii) )) + call message(' Min. temperature directory: ', trim(dirMinTemperature(ii) )) + call message(' Max. temperature directory: ', trim(dirMaxTemperature(ii) )) case(2) call message(' Net radiation directory: ', trim(dirNetRadiation(ii) )) case(3) call message(' Net radiation directory: ', trim(dirNetRadiation(ii) )) - call message(' Abs. vap. press. directory: ', trim(dirabsVapPressure(ii) )) - call message(' Windspeed directory: ', trim(dirwindspeed(ii) )) + call message(' Abs. vap. press. directory: ', trim(dirabsVapPressure(ii) )) + call message(' Windspeed directory: ', trim(dirwindspeed(ii) )) end select call message(' Output directory: ', trim(dirOut(ii) )) if (timeStep_LAI_input < 0) then @@ -337,6 +339,12 @@ PROGRAM mhm_driver call read_soil_moisture(ii) endif + ! read optional spatio-temporal neutrons data + if ( (opti_function .EQ. 17) .AND. optimize ) then + call read_neutrons(ii) + call message(' neutrons data read') + endif + end do ! read optional basin average TWS data at once, therefore outside of the basin loop to ensure same time for all basins @@ -346,7 +354,6 @@ PROGRAM mhm_driver call message(' basin_avg TWS data read') endif - ! The following block is for testing of the restart <<<<<<<<<<<<<<<<<<<<<<<<<< ! if ( write_restart ) then ! itimer = itimer + 1 @@ -364,7 +371,7 @@ PROGRAM mhm_driver ! READ and INITIALISE mRM ROUTING ! -------------------------------------------------------------------------- if (processMatrix(8, 1) .eq. 1) call mrm_init(basin%L0_mask, L0_elev, L0_LCover) -#endif +#endif !this call may be moved to another position as it writes the master config out file for all basins call write_configfile() @@ -375,7 +382,7 @@ PROGRAM mhm_driver iTimer = iTimer + 1 call message() if ( optimize ) then -#ifdef mrm2mhm +#ifdef mrm2mhm ! call optimization against only runoff (no other variables) if ((opti_function .eq. 1) .or. & (opti_function .eq. 2) .or. & @@ -389,13 +396,14 @@ PROGRAM mhm_driver (opti_function .eq. 14)) & call optimization(single_objective_runoff, dirConfigOut, funcBest, maskpara) #endif - + ! call optimization for other variables if ((opti_function .eq. 10) .or. & (opti_function .eq. 11) .or. & (opti_function .eq. 12) .or. & (opti_function .eq. 13) .or. & - (opti_function .eq. 15)) & + (opti_function .eq. 15) .or. & + (opti_function .eq. 17)) & call optimization(objective, dirConfigOut, funcBest, maskpara) ! write a file with final objective function and the best parameter set @@ -403,8 +411,8 @@ PROGRAM mhm_driver ! write a file with final best parameter set in a namlist format call write_optinamelist(processMatrix, global_parameters, maskpara, global_parameters_name(:)) deallocate(maskpara) - else - + else + ! -------------------------------------------------------------------------- ! call mHM ! get runoff timeseries if possible (i.e. when processMatrix(8,1) > 0) @@ -416,7 +424,7 @@ PROGRAM mhm_driver call timer_stop(itimer) call message(' in ', trim(num2str(timer_get(itimer),'(F12.3)')), ' seconds.') - end if + end if ! -------------------------------------------------------------------------- ! WRITE RESTART files @@ -431,7 +439,7 @@ PROGRAM mhm_driver call message(' in ', trim(num2str(timer_get(itimer),'(F9.3)')), ' seconds.') end if -#ifdef mrm2mhm +#ifdef mrm2mhm ! -------------------------------------------------------------------------- ! WRITE RUNOFF (INCLUDING RESTART FILES, has to be called after mHM restart ! files are written) diff --git a/src/mHM/mo_global_variables.f90 b/src/mHM/mo_global_variables.f90 index 75e234c2..f2f216fc 100644 --- a/src/mHM/mo_global_variables.f90 +++ b/src/mHM/mo_global_variables.f90 @@ -2,16 +2,16 @@ !> \brief Global variables ONLY used in reading, writing and startup. -!> \details +!> \details !> \authors Luis Samaniego !> \date Dec 2012 MODULE mo_global_variables - ! This module provides + ! This module provides - ! + ! ! Written Luis Samaniego, Dec 2005 ! Modified Luis Samaniego, Feb 2013 - new variable names, new modules, units ! Rohini Kumar, Jul 2013 - fraction of perfectly sealed area within city added @@ -19,7 +19,7 @@ MODULE mo_global_variables ! Rohini Kumar, Aug 2013 - name changed from "L0_LAI" to "L0_LCover_LAI" ! Rohini Kumar, Aug 2013 - added dirSoil_LUT and dirGeology_LUT ! Luis Samaniego, Nov 2013 - documentation of dimensions - ! Matthias Zink, Nov 2013 - added "InflowGauge" and inflow gauge variabels in basin + ! Matthias Zink, Nov 2013 - added "InflowGauge" and inflow gauge variabels in basin ! Rohini Kumar, May 2014 - added options for the model run cordinate system ! Stephan Thober, Jun 2014 - added timeStep_model_inputs and readPer ! Stephan Thober, Jun 2014 - added perform_mpr, updated restart flags @@ -28,7 +28,7 @@ MODULE mo_global_variables ! Matthias Zink, Mar 2015 - added optional soil mositure readin: dirSoil_moisture, L1_sm ! Stephan Thober, Aug 2015 - moved routing related variables to mRM ! Oldrich Rakovec,Oct 2015 - added definition of basin averaged TWS data - + USE mo_kind, ONLY: i4, i8, dp use mo_common_variables, ONLY: period USE mo_mhm_constants, ONLY: nOutFlxState, YearMonths, maxNoBasins, maxNLCovers @@ -54,11 +54,11 @@ MODULE mo_global_variables ! ! process 6 :: interflow ! ! process 7 :: percolation ! ! process 8 :: routing - ! ! process 9 :: baseflow - ! ! process 10:: neutrons + ! ! process 9 :: baseflow + ! ! process 10:: neutrons integer(i4), dimension(nProcesses, 3), public :: processMatrix ! Info about which process runs in which option and ! ! number of parameters necessary for this option - ! ! col1: process_switch + ! ! col1: process_switch ! ! col2: no. of parameters ! ! col3: cum. no. of parameters @@ -70,15 +70,17 @@ MODULE mo_global_variables real(dp), dimension(:), allocatable, public :: resolutionHydrology ! [m or degree] resolution of hydrology - Level 1 real(dp), dimension(:), allocatable, public :: resolutionRouting ! [m or degree] resolution of routing - Level 11 integer(i4), dimension(:), allocatable, public :: L0_Basin - logical, public :: read_restart ! flag - logical, public :: write_restart ! flag + logical, public :: read_restart ! flag + logical, public :: write_restart ! flag logical, public :: perform_mpr ! switch for performing ! multiscale parameter regionalization character(256),public :: inputFormat_meteo_forcings ! format of meteo input data(bin or nc) ! LAI information character(256), public :: inputFormat_gridded_LAI ! format of gridded LAI data(bin or nc) integer(i4), public :: timeStep_LAI_input ! time step of gridded LAI input - integer(i4), public :: timeStep_sm_input ! time step of optional data: soil moisture sm + ! Optional data + integer(i4), public :: timeStep_sm_input ! time step of optional data: soil moisture sm + integer(i4), public :: timeStep_neutrons_input ! time step of optional data: soil moisture sm integer(i4), public :: iFlag_cordinate_sys ! options model for the run cordinate system ! ------------------------------------------------------------------ @@ -104,8 +106,9 @@ MODULE mo_global_variables character(256), dimension(:), allocatable, public :: dirSoil_moisture ! File of monthly soil moisture character(256), dimension(:), allocatable, public :: fileTWS ! File of tws data - - ! directory common to all basins + character(256), dimension(:), allocatable, public :: dirNeutrons ! File of spatio-temporal neutron data + + ! directory common to all basins character(256), public :: dirConfigOut ! Directory where config run output is written to character(256), public :: dirCommonFiles ! directory where common input files should be located ! ! for all modeled basins @@ -114,13 +117,13 @@ MODULE mo_global_variables ! ------------------------------------------------------------------ real(dp), dimension(:), allocatable, target, public :: yCoor ! GK4 (DHDN3-zone 4) easting real(dp), dimension(:), allocatable, target, public :: xCoor ! GK4 (DHDN3-zone 4) northing - real(dp), dimension(:,:), allocatable, target, public :: lons ! WGS84 lons + real(dp), dimension(:,:), allocatable, target, public :: lons ! WGS84 lons real(dp), dimension(:,:), allocatable, target, public :: lats ! WGS84 lats ! ------------------------------------------------------------------ - ! CONSTANT + ! CONSTANT ! ------------------------------------------------------------------ - real(dp), public :: c2TSTu ! Unit transformation = timeStep/24 + real(dp), public :: c2TSTu ! Unit transformation = timeStep/24 integer(i4), public :: nTstepDay ! Number of time intervals per day ! ! (was previously NAGG) integer(i4), public, parameter :: routingStates = 2 ! [-] Routing states (2=current, 1=past) @@ -142,30 +145,30 @@ MODULE mo_global_variables ! input data integer(i4), dimension(:), allocatable :: id ! Soil Id integer(i4), dimension(:), allocatable :: nHorizons ! Number of horizons - integer(i4), dimension(:), allocatable :: is_present ! Wether this soil type is present in + integer(i4), dimension(:), allocatable :: is_present ! Wether this soil type is present in ! ! this basin or not - real(dp), dimension(:,:), allocatable :: UD ! [mm] Upper Bound of depth - real(dp), dimension(:,:), allocatable :: LD ! [mm] Lower Bound of depth - real(dp), dimension(:,:), allocatable :: clay ! [%] Clay content - real(dp), dimension(:,:), allocatable :: sand ! [%] Sand content + real(dp), dimension(:,:), allocatable :: UD ! [mm] Upper Bound of depth + real(dp), dimension(:,:), allocatable :: LD ! [mm] Lower Bound of depth + real(dp), dimension(:,:), allocatable :: clay ! [%] Clay content + real(dp), dimension(:,:), allocatable :: sand ! [%] Sand content real(dp), dimension(:,:), allocatable :: DbM ! [g/cm2] Mineral Bulk density - real(dp), dimension(:,:), allocatable :: depth ! [mm] Depth of the soil Horizon - real(dp), dimension(:), allocatable :: RZdepth ! [mm] Total soil depth - real(dp), dimension(:,:,:), allocatable :: Wd ! [1] Weights of mHM Horizons according to - ! ! horizons provided in soil database - integer(i4), dimension(:), allocatable :: nTillHorizons ! [1] Number of tillage horizons + real(dp), dimension(:,:), allocatable :: depth ! [mm] Depth of the soil Horizon + real(dp), dimension(:), allocatable :: RZdepth ! [mm] Total soil depth + real(dp), dimension(:,:,:), allocatable :: Wd ! [1] Weights of mHM Horizons according to + ! ! horizons provided in soil database + integer(i4), dimension(:), allocatable :: nTillHorizons ! [1] Number of tillage horizons - ! derived soil hydraulic properties + ! derived soil hydraulic properties real(dp), dimension(:,:,:), allocatable :: thetaS_Till ! [1] Saturated water content of soil horizons ! ! tillage depth - f(OM, management) - real(dp), dimension(:,:), allocatable :: thetaS ! [1] Saturated water content of soil horizons + real(dp), dimension(:,:), allocatable :: thetaS ! [1] Saturated water content of soil horizons ! ! after tillage depth real(dp), dimension(:,:,:), allocatable :: Db ! [g/cm2] Bulk density, LUC dependent ! ! = f( OM, management) - real(dp), dimension(:,:,:), allocatable :: thetaFC_Till ! [1] Field capacity of tillage layers; + real(dp), dimension(:,:,:), allocatable :: thetaFC_Till ! [1] Field capacity of tillage layers; ! ! LUC dependent - f(OM, management) real(dp), dimension(:,:), allocatable :: thetaFC ! [1] Field capacity of deeper layers - real(dp), dimension(:,:,:), allocatable :: thetaPW_Till ! [1] Permament wilting point of tillage layers; + real(dp), dimension(:,:,:), allocatable :: thetaPW_Till ! [1] Permament wilting point of tillage layers; ! ! LUC dependent - f(OM, management) real(dp), dimension(:,:), allocatable :: thetaPW ! [1] Permanent wilting point of deeper layers real(dp), dimension(:,:,:), allocatable :: Ks ! [cm/d] Saturated hydaulic conductivity @@ -198,7 +201,7 @@ MODULE mo_global_variables ! ----------------------------------------------------------------- ! Land cover information real(dp), public :: fracSealed_cityArea ! fraction of area within city assumed to be - ! perfectly sealed [0-1] + ! perfectly sealed [0-1] integer(i4), public :: nLCoverScene ! Number of land cover scene character(256), public, dimension(:), allocatable :: LCfilename ! file names for the different land cover scenes integer(i4), public, dimension(:,:), allocatable :: LCyearId ! Mapping of landcover scenes (1, 2, ...) for each basin @@ -208,7 +211,7 @@ MODULE mo_global_variables integer(i4), public :: nLAIclass ! Number of LAI classes integer(i4), public, dimension(:), allocatable :: LAIUnitList ! List of ids of each LAI class in LAILUT real(dp), public, dimension(:,:), allocatable :: LAILUT ! [m2/m2] Leaf area index for LAIUnit - ! ! dim1=land cover class, dim2=month of year + ! ! dim1=land cover class, dim2=month of year ! ------------------------------------------------------------------- ! GRID description @@ -233,7 +236,7 @@ MODULE mo_global_variables real(dp), dimension(:), allocatable, public :: L1_latitude ! 1d latitude array for active grid cells real(dp), dimension(:), allocatable, public :: L1_rect_longitude ! 1d longitude array for whole basin rectangle real(dp), dimension(:), allocatable, public :: L1_rect_latitude ! 1d latitude array for whole basin rectangle - + ! ------------------------------------------------------------------- ! PERIOD description ! ------------------------------------------------------------------- @@ -249,7 +252,7 @@ MODULE mo_global_variables integer(i4), public :: nBasins ! Number of basins for multi-basin optimization type basinInfo ! dim1 = basin id - ! for remapping + ! for remapping integer(i4), dimension(:), allocatable :: L0_iStart ! Starting cell index of a given basin at L0 integer(i4), dimension(:), allocatable :: L0_iEnd ! Ending cell index of a given basin at L0 integer(i4), dimension(:), allocatable :: L0_iStartMask ! Starting cell index of mask a given basin at L0 @@ -273,51 +276,51 @@ MODULE mo_global_variables logical, dimension(:), allocatable, target :: L0_mask ! target variable for mRM ! ------------------------------------------------------------------- - ! L0 DOMAIN description -> + ! L0 DOMAIN description -> ! ------------------------------------------------------------------- ! dim1 = number grid cells ! input data - morphological variables - real(dp), public, dimension(:), allocatable, target :: L0_elev ! [m] Elevation (sinks removed) + real(dp), public, dimension(:), allocatable, target :: L0_elev ! [m] Elevation (sinks removed) ! ! target variable for coupling to mRM - real(dp), public, dimension(:), allocatable :: L0_slope ! [%] Slope + real(dp), public, dimension(:), allocatable :: L0_slope ! [%] Slope real(dp), public, dimension(:), allocatable :: L0_asp ! [degree] Aspect degree integer(i4), public, dimension(:), allocatable :: L0_soilId ! Soil id integer(i4), public, dimension(:), allocatable :: L0_geoUnit ! Geologic formation (unit) - + ! input data - land cover integer(i4), public, dimension(:), allocatable :: L0_LCover_LAI ! Special landcover id (upto 10 classes) only for the LAI index - integer(i4), public, dimension(:,:), allocatable, target :: L0_LCover ! Normal landcover id (upto 3 classes) + integer(i4), public, dimension(:,:), allocatable, target :: L0_LCover ! Normal landcover id (upto 3 classes) ! ! dim1=number grid cells, dim2=Number of land cover scenes ! ! target variable for coupling to mRM ! mHM derived variables ! dim1 = number grid cells L0 - integer(i4), public :: L0_nCells ! Number of valid cells + integer(i4), public :: L0_nCells ! Number of valid cells integer(i4), public, dimension(:,:), allocatable :: L0_cellCoor ! Cell coordinates (row,col) for each grid cell, dim2=2 integer(i4), public, dimension(:), allocatable :: L0_Id ! Level-0 id real(dp), public, dimension(:), allocatable :: L0_slope_emp ! Empirical quantiles of slope ! real(dp), public, dimension(:,:), allocatable :: L0_gridded_LAI ! gridded LAI data used when timeStep_LAI_input<0 ! ! dim1=number of grid cells, dim2=number of LAI time steps - real(dp), public, dimension(:), allocatable :: L0_areaCell ! [m2] Area of a cell at level-0 + real(dp), public, dimension(:), allocatable :: L0_areaCell ! [m2] Area of a cell at level-0 ! ------------------------------------------------------------------- ! L1 DOMAIN description ! ------------------------------------------------------------------- ! dim1 = number grid cells L1 ! dim2 = 2 - integer(i4), public :: L1_nCells ! No. of cells at this level + integer(i4), public :: L1_nCells ! No. of cells at this level integer(i4), public, dimension(:,:), allocatable :: L1_cellCoor ! Cell coordinates (row,col) ! ! -> L1 modelling integer(i4), public, dimension(:), allocatable :: L1_Id ! Level-1 id real(dp), public, dimension(:), allocatable :: L1_areaCell ! [km2] Effective area of cell at this level - integer(i4), public, dimension(:), allocatable :: L1_upBound_L0 ! Row start at finer level-0 scale - integer(i4), public, dimension(:), allocatable :: L1_downBound_L0 ! Row end at finer level-0 scale - integer(i4), public, dimension(:), allocatable :: L1_leftBound_L0 ! Col start at finer level-0 scale - integer(i4), public, dimension(:), allocatable :: L1_rightBound_L0 ! Col end at finer level-0 scale + integer(i4), public, dimension(:), allocatable :: L1_upBound_L0 ! Row start at finer level-0 scale + integer(i4), public, dimension(:), allocatable :: L1_downBound_L0 ! Row end at finer level-0 scale + integer(i4), public, dimension(:), allocatable :: L1_leftBound_L0 ! Col start at finer level-0 scale + integer(i4), public, dimension(:), allocatable :: L1_rightBound_L0 ! Col end at finer level-0 scale integer(i4), public, dimension(:), allocatable :: L1_nTCells_L0 ! Total number of valid L0 cells - ! ! in a given L1 cell + ! ! in a given L1 cell ! Forcings ! dim1 = number grid cells L1 ! dim2 = number of meteorological time steps @@ -335,37 +338,40 @@ MODULE mo_global_variables ! dim2 = number of meteorological time steps ! soil moisture real(dp), public, dimension(:,:), allocatable :: L1_sm ! [-] soil moisture input for optimization - logical, public, dimension(:,:), allocatable :: L1_sm_mask ! [-] mask for valid data in L1_sm + logical, public, dimension(:,:), allocatable :: L1_sm_mask ! [-] mask for valid data in L1_sm integer(i4) :: nTimeSteps_L1_sm ! [-] number of time steps in L1_sm_mask - integer(i4) :: nSoilHorizons_sm_input ! No. of mhm soil horizons equivalent to - ! ! soil moisture input + integer(i4) :: nSoilHorizons_sm_input ! No. of mhm soil horizons equivalent to sm input + ! neutrons + real(dp), public, dimension(:,:), allocatable :: L1_neutronsdata ! [cph] ground albedo neutrons input + logical, public, dimension(:,:), allocatable :: L1_neutronsdata_mask ! [cph] mask for valid data in L1_neutrons + integer(i4) :: nTimeSteps_L1_neutrons ! [-] number of time steps in L1_neutrons_mask ! Land cover ! dim1 = number grid cells L1 real(dp), public, dimension(:), allocatable :: L1_fSealed ! [1] Fraction of sealed area real(dp), public, dimension(:), allocatable :: L1_fForest ! [1] Fraction of forest cover real(dp), public, dimension(:), allocatable :: L1_fPerm ! [1] Fraction of permeable cover - + ! State variables ! dim1 = number grid cells L1 ! dim2 = number model soil horizons - real(dp), public, dimension(:), allocatable :: L1_inter ! [mm] Interception + real(dp), public, dimension(:), allocatable :: L1_inter ! [mm] Interception real(dp), public, dimension(:), allocatable :: L1_snowPack ! [mm] Snowpack real(dp), public, dimension(:), allocatable :: L1_sealSTW ! [mm] Retention storage of impervious areas real(dp), public, dimension(:,:), allocatable :: L1_soilMoist ! [mm] Soil moisture of each horizon real(dp), public, dimension(:), allocatable :: L1_unsatSTW ! [mm] upper soil storage real(dp), public, dimension(:), allocatable :: L1_satSTW ! [mm] groundwater storage - real(dp), public, dimension(:), allocatable :: L1_neutrons ! [cph] ground albedo neutrons + real(dp), public, dimension(:), allocatable :: L1_neutrons ! [mm] Ground Albedo Neutrons ! Fluxes ! dim1 = number grid cells L1 ! dim2 = number model soil horizons - real(dp), public, dimension(:), allocatable :: L1_pet_calc ! [mm TS-1] estimated/corrected potential evapotranspiration + real(dp), public, dimension(:), allocatable :: L1_pet_calc ! [mm TS-1] estimated/corrected potential evapotranspiration real(dp), public, dimension(:,:), allocatable :: L1_aETSoil ! [mm TS-1] Actual ET from soil layers real(dp), public, dimension(:), allocatable :: L1_aETCanopy ! [mm TS-1] Real evaporation intensity from canopy real(dp), public, dimension(:), allocatable :: L1_aETSealed ! [mm TS-1] Real evap. from free water surfaces real(dp), public, dimension(:), allocatable :: L1_baseflow ! [mm TS-1] Baseflow - real(dp), public, dimension(:,:), allocatable :: L1_infilSoil ! [mm TS-1] Infiltration intensity each soil horizon + real(dp), public, dimension(:,:), allocatable :: L1_infilSoil ! [mm TS-1] Infiltration intensity each soil horizon real(dp), public, dimension(:), allocatable :: L1_fastRunoff ! [mm TS-1] Fast runoff component real(dp), public, dimension(:), allocatable :: L1_melt ! [mm TS-1] Melting snow depth. real(dp), public, dimension(:), allocatable :: L1_percol ! [mm TS-1] Percolation. @@ -374,7 +380,7 @@ MODULE mo_global_variables real(dp), public, dimension(:), allocatable :: L1_runoffSeal ! [mm TS-1] Direct runoff from impervious areas real(dp), public, dimension(:), allocatable :: L1_slowRunoff ! [mm TS-1] Slow runoff component real(dp), public, dimension(:), allocatable :: L1_snow ! [mm TS-1] Snow precipitation depth - real(dp), public, dimension(:), allocatable :: L1_Throughfall ! [mm TS-1] Throughfall. + real(dp), public, dimension(:), allocatable :: L1_Throughfall ! [mm TS-1] Throughfall. real(dp), public, dimension(:), allocatable :: L1_total_runoff ! [m3 TS-1] Generated runoff ! Effective parameters @@ -382,8 +388,8 @@ MODULE mo_global_variables ! dim2 = number model soil horizons real(dp), public, dimension(:), allocatable :: L1_alpha ! [1] Exponent for the upper reservoir real(dp), public, dimension(:), allocatable :: L1_degDayInc ! [d-1 degC-1] Increase of the Degree-day factor - ! ! per mm of increase in precipitation - real(dp), public, dimension(:), allocatable :: L1_degDayMax ! [mm-1 degC-1] Maximum Degree-day factor + ! ! per mm of increase in precipitation + real(dp), public, dimension(:), allocatable :: L1_degDayMax ! [mm-1 degC-1] Maximum Degree-day factor real(dp), public, dimension(:), allocatable :: L1_degDayNoPre ! [mm-1 degC-1] Degree-day factor with no precipitation. real(dp), public, dimension(:), allocatable :: L1_degDay ! [mm d-1degC-1] Degree-day factor. real(dp), public, dimension(:), allocatable :: L1_karstLoss ! [1] Karstic percolation loss @@ -393,22 +399,22 @@ MODULE mo_global_variables real(dp), public, dimension(:,:), allocatable :: L1_aeroResist ! [s m-1] aerodynamical resitance real(dp), public, dimension(:,:), allocatable :: L1_surfResist ! [s m-1] bulk surface resitance ! ! dim1 = No cells for basin, dim2 = No of Months in year - real(dp), public, dimension(:,:), allocatable :: L1_fRoots ! [1] Fraction of roots in soil horizons - real(dp), public, dimension(:), allocatable :: L1_maxInter ! [mm] Maximum interception - real(dp), public, dimension(:), allocatable :: L1_kfastFlow ! [d-1] Fast interflow recession coefficient - real(dp), public, dimension(:), allocatable :: L1_kSlowFlow ! [d-1] Slow interflow recession coefficient - real(dp), public, dimension(:), allocatable :: L1_kBaseFlow ! [d-1] Baseflow recession coefficient + real(dp), public, dimension(:,:), allocatable :: L1_fRoots ! [1] Fraction of roots in soil horizons + real(dp), public, dimension(:), allocatable :: L1_maxInter ! [mm] Maximum interception + real(dp), public, dimension(:), allocatable :: L1_kfastFlow ! [d-1] Fast interflow recession coefficient + real(dp), public, dimension(:), allocatable :: L1_kSlowFlow ! [d-1] Slow interflow recession coefficient + real(dp), public, dimension(:), allocatable :: L1_kBaseFlow ! [d-1] Baseflow recession coefficient real(dp), public, dimension(:), allocatable :: L1_kPerco ! [d-1] percolation coefficient real(dp), public, dimension(:,:), allocatable :: L1_soilMoistFC ! [mm] Soil moisture below which actual ET ! ! is reduced linearly till PWP - real(dp), public, dimension(:,:), allocatable :: L1_soilMoistSat ! [mm] Saturation soil moisture for each horizon [mm] + real(dp), public, dimension(:,:), allocatable :: L1_soilMoistSat ! [mm] Saturation soil moisture for each horizon [mm] real(dp), public, dimension(:,:), allocatable :: L1_soilMoistExp ! [1] Exponential parameter to how non-linear ! ! is the soil water retention real(dp), public, dimension(:), allocatable :: L1_tempThresh ! [degC] Threshold temperature for snow/rain real(dp), public, dimension(:), allocatable :: L1_unsatThresh ! [mm] Threshhold water depth controlling fast interflow real(dp), public, dimension(:), allocatable :: L1_sealedThresh ! [mm] Threshhold water depth for surface runoff ! ! in sealed surfaces - real(dp), public, dimension(:,:), allocatable :: L1_wiltingPoint ! [mm] Permanent wilting point: below which neither + real(dp), public, dimension(:,:), allocatable :: L1_wiltingPoint ! [mm] Permanent wilting point: below which neither ! ! plant can take water nor water can drain in ! ------------------------------------------------------------------- @@ -425,12 +431,12 @@ MODULE mo_global_variables real(dp), public, dimension(int(YearMonths,i4)) :: fnight_temp ! [-] Night factor mean temp ! ------------------------------------------------------------------- - ! DEFINE OUTPUTS + ! DEFINE OUTPUTS ! ------------------------------------------------------------------- ! integer(i4) :: timeStep_model_outputs ! timestep for writing model outputs logical, dimension(nOutFlxState) :: outputFlxState ! Define model outputs see "mhm_outputs.nml" - ! dim1 = number of output variables to be written + ! dim1 = number of output variables to be written ! - + END MODULE mo_global_variables diff --git a/src/mHM/mo_init_states.f90 b/src/mHM/mo_init_states.f90 index 0089db7a..22601f70 100644 --- a/src/mHM/mo_init_states.f90 +++ b/src/mHM/mo_init_states.f90 @@ -2,7 +2,7 @@ !> \brief Initialization of all state variables of mHM. -!> \details This module initializes all state variables required to run mHM. +!> \details This module initializes all state variables required to run mHM. !> Two options are provided: !> - (1) default values !> - (2) from nc file @@ -38,14 +38,14 @@ MODULE mo_init_states !> \brief Allocation of space for mHM related L1 and L11 variables. !> \details Allocation of space for mHM related L1 and L11 variables (e.g., states, - !> fluxes, and parameters) for a given basin. Variables allocated here is + !> fluxes, and parameters) for a given basin. Variables allocated here is !> defined in them mo_global_variables.f90 file. After allocating any variable - !> in this routine, initalize them in the following variables_default_init + !> in this routine, initalize them in the following variables_default_init !> subroutine: - !> + !> ! ! CALLING SEQUENCE - ! call variables_alloc(iBasin) + ! call variables_alloc(iBasin) ! INTENT(IN) !> \param[in] "integer(i4) :: iBasin" - basin id @@ -78,18 +78,18 @@ MODULE mo_init_states ! Modified, R. Kumar, Sep 2013 - documentation added according to the template ! S. Thober, Aug 2015 - removed routing related variables - subroutine variables_alloc(iBasin) + subroutine variables_alloc(iBasin) use mo_global_variables, only: nSoilHorizons_mHM, & - L1_fSealed, L1_fForest, L1_fPerm, L1_inter, L1_snowPack, L1_sealSTW, & + L1_fSealed, L1_fForest, L1_fPerm, L1_inter, L1_snowPack, L1_sealSTW, & L1_soilMoist, L1_unsatSTW, L1_satSTW, & L1_pet_calc, L1_aETSoil, L1_aETCanopy, L1_aETSealed, & - L1_baseflow, L1_infilSoil, L1_fastRunoff, L1_melt, & - L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, L1_slowRunoff, & - L1_snow, L1_Throughfall, L1_total_runoff, L1_alpha, L1_degDayInc, & + L1_baseflow, L1_infilSoil, L1_fastRunoff, L1_melt, & + L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, L1_slowRunoff, & + L1_snow, L1_Throughfall, L1_total_runoff, L1_alpha, L1_degDayInc, & L1_degDayMax, L1_degDayNoPre, L1_degDay, L1_karstLoss, L1_fAsp, & L1_HarSamCoeff, L1_PrieTayAlpha, L1_aeroResist, L1_surfResist, & - L1_fRoots, L1_maxInter, L1_kfastFlow, L1_kSlowFlow, L1_kBaseFlow, & + L1_fRoots, L1_maxInter, L1_kfastFlow, L1_kSlowFlow, L1_kBaseFlow, & L1_kPerco, L1_soilMoistFC, L1_soilMoistSat, L1_soilMoistExp, & L1_tempThresh, L1_unsatThresh, L1_sealedThresh, L1_wiltingPoint, & L1_neutrons @@ -109,7 +109,7 @@ subroutine variables_alloc(iBasin) real(dp), dimension(:,:), allocatable :: dummy_Matrix_months ! level-1 information - call get_basin_info( iBasin, 1, nrows1, ncols1, ncells=ncells1 ) + call get_basin_info( iBasin, 1, nrows1, ncols1, ncells=ncells1 ) ! for appending and intialization allocate( dummy_Vector (nCells1) ) @@ -139,7 +139,7 @@ subroutine variables_alloc(iBasin) dummy_Vector(:) = 0.0_dp call append( L1_sealSTW, dummy_Vector ) - ! Soil moisture of each horizon + ! Soil moisture of each horizon dummy_Matrix(:,:) = 0.0_dp call append( L1_soilMoist, dummy_Matrix ) @@ -150,7 +150,7 @@ subroutine variables_alloc(iBasin) ! groundwater storage dummy_Vector(:) = 0.0_dp call append( L1_satSTW, dummy_Vector ) - + ! ground albedo neutrons dummy_Vector(:) = 0.0_dp call append( L1_neutrons, dummy_Vector ) @@ -162,7 +162,7 @@ subroutine variables_alloc(iBasin) ! calculated / corrected potential evapotranspiration dummy_Vector(:) = 0.0_dp call append( L1_pet_calc, dummy_Vector ) - + ! soil actual ET dummy_Matrix(:,:) = 0.0_dp call append( L1_aETSoil, dummy_Matrix ) @@ -207,19 +207,19 @@ subroutine variables_alloc(iBasin) dummy_Vector(:) = 0.0_dp call append( L1_runoffSeal, dummy_Vector ) - ! slow runoff + ! slow runoff dummy_Vector(:) = 0.0_dp call append( L1_slowRunoff, dummy_Vector ) - ! snow (solid water) + ! snow (solid water) dummy_Vector(:) = 0.0_dp call append( L1_snow, dummy_Vector ) - ! throughfall + ! throughfall dummy_Vector(:) = 0.0_dp call append( L1_Throughfall, dummy_Vector ) - ! throughfall + ! throughfall dummy_Vector(:) = 0.0_dp call append( L1_total_runoff, dummy_Vector ) @@ -235,7 +235,7 @@ subroutine variables_alloc(iBasin) dummy_Vector(:) = 0.0_dp call append( L1_degDayInc, dummy_Vector ) - ! maximum degree-day factor + ! maximum degree-day factor dummy_Vector(:) = 0.0_dp call append( L1_degDayMax, dummy_Vector ) @@ -271,27 +271,27 @@ subroutine variables_alloc(iBasin) dummy_Matrix_months = 0.0_dp call append( L1_surfResist, dummy_Matrix_months ) - ! Fraction of roots in soil horizons + ! Fraction of roots in soil horizons dummy_Matrix(:,:) = 0.0_dp call append( L1_fRoots, dummy_Matrix ) - ! Maximum interception + ! Maximum interception dummy_Vector(:) = 0.0_dp call append( L1_maxInter, dummy_Vector(:) ) - ! fast interflow recession coefficient + ! fast interflow recession coefficient dummy_Vector(:) = 0.0_dp call append( L1_kfastFlow, dummy_Vector ) - ! slow interflow recession coefficient + ! slow interflow recession coefficient dummy_Vector(:) = 0.0_dp call append( L1_kSlowFlow, dummy_Vector ) - ! baseflow recession coefficient + ! baseflow recession coefficient dummy_Vector(:) = 0.0_dp call append( L1_kBaseFlow, dummy_Vector ) - ! percolation coefficient + ! percolation coefficient dummy_Vector(:) = 0.0_dp call append( L1_kPerco, dummy_Vector ) @@ -307,7 +307,7 @@ subroutine variables_alloc(iBasin) dummy_Matrix(:,:) = 0.0_dp call append( L1_soilMoistExp, dummy_Matrix ) - ! Threshold temperature for snow/rain + ! Threshold temperature for snow/rain dummy_Vector(:) = 0.0_dp call append( L1_tempThresh, dummy_Vector ) @@ -324,27 +324,27 @@ subroutine variables_alloc(iBasin) call append( L1_wiltingPoint, dummy_Matrix ) ! free space - if ( allocated( dummy_Vector ) ) deallocate( dummy_Vector ) + if ( allocated( dummy_Vector ) ) deallocate( dummy_Vector ) if ( allocated( dummy_Matrix ) ) deallocate( dummy_Matrix ) if ( allocated( dummy_Matrix_months ) ) deallocate( dummy_Matrix_months ) end subroutine variables_alloc - + ! ------------------------------------------------------------------ ! NAME ! variables_default_init !> \brief Default initalization mHM related L1 variables - + !> \details Default initalization of mHM related L1 variables (e.g., states, !> fluxes, and parameters) as per given constant values given in mo_mhm_constants. - !> Variables initalized here is defined in the mo_global_variables.f90 file. - !> Only Variables that are defined in the variables_alloc subroutine are + !> Variables initalized here is defined in the mo_global_variables.f90 file. + !> Only Variables that are defined in the variables_alloc subroutine are !> intialized here. ! !> If a variable is added or removed here, then it also has to be added or removed - !> in the subroutine state_variables_set in the module mo_restart and in the + !> in the subroutine state_variables_set in the module mo_restart and in the !> subroutine set_state in the module mo_set_netcdf_restart. ! ! CALLING SEQUENCE @@ -384,15 +384,15 @@ subroutine variables_default_init() use mo_global_variables, only: & nSoilHorizons_mHM, HorizonDepth_mHM, & - L1_fSealed, L1_fForest, L1_fPerm, L1_inter, L1_snowPack, L1_sealSTW, & + L1_fSealed, L1_fForest, L1_fPerm, L1_inter, L1_snowPack, L1_sealSTW, & L1_soilMoist, L1_unsatSTW, L1_satSTW, & L1_pet_calc, L1_aETSoil, L1_aETCanopy, L1_aETSealed, & L1_baseflow, L1_infilSoil, L1_fastRunoff, L1_melt, & - L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, L1_slowRunoff, & - L1_snow, L1_Throughfall, L1_total_runoff, L1_alpha, L1_degDayInc, & - L1_degDayMax, L1_degDayNoPre, L1_degDay, L1_karstLoss, L1_fAsp, & + L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, L1_slowRunoff, & + L1_snow, L1_Throughfall, L1_total_runoff, L1_alpha, L1_degDayInc, & + L1_degDayMax, L1_degDayNoPre, L1_degDay, L1_karstLoss, L1_fAsp, & L1_HarSamCoeff, L1_PrieTayAlpha, L1_aeroResist, L1_surfResist, & - L1_fRoots, L1_maxInter, L1_kfastFlow, L1_kSlowFlow, L1_kBaseFlow, & + L1_fRoots, L1_maxInter, L1_kfastFlow, L1_kSlowFlow, L1_kBaseFlow, & L1_kPerco, L1_soilMoistFC, L1_soilMoistSat, L1_soilMoistExp, & L1_tempThresh, L1_unsatThresh, L1_sealedThresh, L1_wiltingPoint, & L1_neutrons @@ -410,7 +410,7 @@ subroutine variables_default_init() !------------------------------------------- ! LAND COVER variables !------------------------------------------- - L1_fSealed = P1_InitStateFluxes + L1_fSealed = P1_InitStateFluxes L1_fForest = P1_InitStateFluxes L1_fPerm = P1_InitStateFluxes @@ -426,7 +426,7 @@ subroutine variables_default_init() !Retention storage of impervious areas L1_sealSTW = P1_InitStateFluxes - ! Soil moisture of each horizon + ! Soil moisture of each horizon do i = 1, nSoilHorizons_mHM-1 if (i .eq. 1) then L1_soilMoist(:,i) = HorizonDepth_mHM(i)*C1_InitStateSM @@ -442,7 +442,7 @@ subroutine variables_default_init() ! groundwater storage L1_satSTW = P4_InitStateFluxes - + ! ground albedo neutrons, initially zero L1_neutrons = P1_InitStateFluxes @@ -452,7 +452,7 @@ subroutine variables_default_init() ! corrected / calculated potential ET L1_pet_calc = P1_InitStateFluxes - + ! soil actual ET L1_aETSoil = P1_InitStateFluxes @@ -486,16 +486,16 @@ subroutine variables_default_init() ! runoff from impervious area L1_runoffSeal = P1_InitStateFluxes - ! slow runoff + ! slow runoff L1_slowRunoff = P1_InitStateFluxes - ! snow (solid water) + ! snow (solid water) L1_snow = P1_InitStateFluxes - ! throughfall + ! throughfall L1_Throughfall = P1_InitStateFluxes - ! throughfall + ! throughfall L1_total_runoff = P1_InitStateFluxes !------------------------------------------- @@ -508,7 +508,7 @@ subroutine variables_default_init() ! increase of the Degree-day factor per mm of increase in precipitation L1_degDayInc = P1_InitStateFluxes - ! maximum degree-day factor + ! maximum degree-day factor L1_degDayMax = P1_InitStateFluxes ! degree-day factor with no precipitation @@ -529,29 +529,29 @@ subroutine variables_default_init() ! PET Priestley Taylor coefficient L1_PrieTayAlpha = P1_InitStateFluxes - ! PET aerodynamical resistance + ! PET aerodynamical resistance L1_aeroResist = P1_InitStateFluxes - ! PET bulk surface resistance + ! PET bulk surface resistance L1_surfResist = P1_InitStateFluxes - ! Fraction of roots in soil horizons + ! Fraction of roots in soil horizons L1_fRoots = P1_InitStateFluxes - ! Maximum interception + ! Maximum interception L1_maxInter = P1_InitStateFluxes - ! fast interflow recession coefficient + ! fast interflow recession coefficient L1_kfastFlow = P1_InitStateFluxes - ! slow interflow recession coefficient + ! slow interflow recession coefficient L1_kSlowFlow = P1_InitStateFluxes - ! baseflow recession coefficient + ! baseflow recession coefficient L1_kBaseFlow = P1_InitStateFluxes - ! percolation coefficient + ! percolation coefficient L1_kPerco = P1_InitStateFluxes ! Soil moisture below which actual ET is reduced linearly till PWP @@ -563,7 +563,7 @@ subroutine variables_default_init() ! Exponential parameter to how non-linear is the soil water retention L1_soilMoistExp = P1_InitStateFluxes - ! Threshold temperature for snow/rain + ! Threshold temperature for snow/rain L1_tempThresh = P1_InitStateFluxes ! Threshhold water depth controlling fast interflow @@ -582,13 +582,13 @@ end subroutine variables_default_init ! get_basin_info !> \brief Get basic basin information (e.g., nrows, ncols, indices, mask) - + !> \details Get basic basin information (e.g., nrows, ncols, indices, mask) for !> different levels (L0, L1, and L2). ! ! CALLING SEQUENCE ! call get_basin_info(iBasin, iLevel,nrows,ncols, ncells, iStart, iEnd, & - ! iStartMask, iEndMask, mask, xllcorner, yllcorner, cellsize) + ! iStartMask, iEndMask, mask, xllcorner, yllcorner, cellsize) ! INTENT(IN) !> \param[in] "integer(i4) :: iBasin" basin id @@ -601,7 +601,7 @@ end subroutine variables_default_init !> \param[out] "integer(i4) :: nRows" no. of rows !> \param[out] "integer(i4) :: nCols" no. of coloums !> \param[out] "integer(i4) :: ncells" no. of cells - !> \param[out] "integer(i4) :: iStart" start cell index of a given basin at a given level + !> \param[out] "integer(i4) :: iStart" start cell index of a given basin at a given level !> \param[out] "integer(i4) :: iEnd" end cell index of a given basin at a given level !> \param[out] "integer(i4) :: iStartMask" start cell index of mask a given basin at a given level !> \param[out] "integer(i4) :: iEndMask" end cell index of mask a given basin at a given level @@ -612,7 +612,7 @@ end subroutine variables_default_init ! INTENT(INOUT), OPTIONAL ! None - ! INTENT(OUT), OPTIONAL + ! INTENT(OUT), OPTIONAL !> \param[out] "logical, optional :: mask" Mask at a given level !> \param[out] "real(dp), optional :: xllcorner" x coordinate of the lowerleft corner at a given level !> \param[out] "real(dp), optional :: yllcorner" y coordinate of the lowerleft corner at a given level @@ -633,7 +633,7 @@ end subroutine variables_default_init ! Stephan Thober, Aug 2015 - moved L11 and L110 to mRM subroutine get_basin_info(iBasin, iLevel, nrows, ncols, ncells, iStart, iEnd, & - iStartMask, iEndMask, mask, xllcorner, yllcorner, cellsize) + iStartMask, iEndMask, mask, xllcorner, yllcorner, cellsize) use mo_message, only: message use mo_global_variables, only: basin, level0, level1, level2 @@ -724,19 +724,19 @@ end subroutine get_basin_info ! calculate_grid_properties ! PURPOSE - !> \brief Calculates basic grid properties at a required coarser level using + !> \brief Calculates basic grid properties at a required coarser level using !> information of a given finer level. - !> \brief Calculates basic grid properties at a required coarser level (e.g., L11) using - !> information of a given finer level (e.g., L0). Basic grid properties such as + !> \brief Calculates basic grid properties at a required coarser level (e.g., L11) using + !> information of a given finer level (e.g., L0). Basic grid properties such as !> nrows, ncols, xllcorner, yllcorner cellsize are estimated in this - !> routine. + !> routine. ! CALLING SEQUENCE ! call calculate_grid_properties( nrowsIn, ncolsIn, xllcornerIn, & ! yllcornerIn, cellsizeIn, nodata_valueIn, & ! aimingResolution, nrowsOut, ncolsOut, xllcornerOut, & - ! yllcornerOut, cellsizeOut, nodata_valueOut ) + ! yllcornerOut, cellsizeOut, nodata_valueOut ) ! INTENT(IN) !> \param[in] "integer(i4) :: nrowsIn" no. of rows at an input level !> \param[in] "integer(i4) :: ncolsIn" no. of cols at an input level @@ -786,7 +786,7 @@ end subroutine get_basin_info subroutine calculate_grid_properties( & nrowsIn, ncolsIn, xllcornerIn, yllcornerIn, cellsizeIn, nodata_valueIn, & aimingResolution, & - nrowsOut, ncolsOut, xllcornerOut, yllcornerOut, cellsizeOut, nodata_valueOut ) + nrowsOut, ncolsOut, xllcornerOut, yllcornerOut, cellsizeOut, nodata_valueOut ) use mo_message, only: message ! for print out use mo_string_utils, only: num2str @@ -805,7 +805,7 @@ subroutine calculate_grid_properties( & integer(i4), intent(out) :: ncolsOut real(dp), intent(out) :: xllcornerOut real(dp), intent(out) :: yllcornerOut - real(dp), intent(out) :: cellsizeOut + real(dp), intent(out) :: cellsizeOut real(dp), intent(out) :: nodata_valueOut ! local variables @@ -824,8 +824,8 @@ subroutine calculate_grid_properties( & cellsizeOut = cellsizeIn * cellFactor ncolsOut = ceiling( real(ncolsIn, dp) / cellFactor) nrowsOut = ceiling( real(nrowsIn, dp) / cellFactor) - xllcornerOut = xllcornerIn + real(ncolsIn,dp) * cellsizeIn - real(ncolsOut,dp) * cellsizeOut - yllcornerOut = yllcornerIn + real(nrowsIn,dp) * cellsizeIn - real(nrowsOut,dp) * cellsizeOut + xllcornerOut = xllcornerIn + real(ncolsIn,dp) * cellsizeIn - real(ncolsOut,dp) * cellsizeOut + yllcornerOut = yllcornerIn + real(nrowsIn,dp) * cellsizeIn - real(nrowsOut,dp) * cellsizeOut nodata_valueOut = nodata_valueIn end subroutine calculate_grid_properties diff --git a/src/mHM/mo_mhm.f90 b/src/mHM/mo_mhm.f90 index 0691fbc8..e7495367 100644 --- a/src/mHM/mo_mhm.f90 +++ b/src/mHM/mo_mhm.f90 @@ -9,26 +9,26 @@ !> The processes are executed in ascending order. At the moment only !> process 5 and 8 have options.\n !> The MPR technique is only called either if the land cover has been -!> changed or for very first time step.\n +!> changed or for very first time step.\n !> !> Currently the following processes are implemented: \n !> -!> Process | Name | Flag | Description +!> Process | Name | Flag | Description !> ---------- | ------------------------- | ----- | ------------------------------------------ -!> 1 | interception | 1 | Maximum interception +!> 1 | interception | 1 | Maximum interception !> 2 | snow and melting | 1 | Degree-day -!> 3 | soil moisture | 1 | Infiltration capacity, Brooks-Corey -!> 4 | direct runoff | 1 | Linear reservoir exceedance -!> 5 | PET | 0 | PET is read as input +!> 3 | soil moisture | 1 | Infiltration capacity, Brooks-Corey +!> 4 | direct runoff | 1 | Linear reservoir exceedance +!> 5 | PET | 0 | PET is read as input !> 5 | '' | 1 | Hargreaves-Samani !> 5 | '' | 2 | Priestley-Taylor !> 5 | '' | 3 | Penman-Monteith !> 6 | interflow | 1 | Nonlinear reservoir with saturation excess -!> 7 | percolation and base flow | 1 | GW linear reservoir +!> 7 | percolation and base flow | 1 | GW linear reservoir !> 8 | routing | 0 | no routing -!> 8 | '' | 1 | use mRM i.e. Muskingum +!> 8 | '' | 1 | use mRM i.e. Muskingum !> - + !> \author Luis Samaniego !> \date Dec 2012 @@ -36,9 +36,9 @@ MODULE mo_mHM use mo_kind, only: i4, dp use mo_mhm_constants, only: nodata_dp - use mo_message, only: message + use mo_message, only: message !$ USE omp_lib - + IMPLICIT NONE PUBLIC :: mHM ! initialization sequence @@ -94,18 +94,18 @@ MODULE mo_mHM ! Modified Luis Samaniego, Rohini Kumar, Dec 2012 - modularization ! Luis Samaniego, Feb 2013 - call routine - ! Rohini Kumar, Feb 2013 - MPR call and other pre-requisite + ! Rohini Kumar, Feb 2013 - MPR call and other pre-requisite ! variables for this call ! Rohini Kumar, May 2013 - Error checks ! Rohini Kumar, Jun 2013 - sealed area correction in total runoff ! - initalization of soil moist. at first timestep - ! Rohini Kumar, Aug 2013 - dynamic LAI option included, and changed within + ! Rohini Kumar, Aug 2013 - dynamic LAI option included, and changed within ! the code made accordingly (e.g., canopy intecpt.) ! - max. canopy interception is estimated outside of MPR ! call ! Matthias Zink, Feb 2014 - added PET calculation: Hargreaves-Samani (Process 5) ! Matthias Zink, Mar 2014 - added inflow from upstream areas - ! Matthias Zink, Apr 2014 - added PET calculation: Priestley-Taylor and Penamn-Monteith + ! Matthias Zink, Apr 2014 - added PET calculation: Priestley-Taylor and Penamn-Monteith ! and its parameterization (Process 5) ! Rohini Kumar, Apr 2014 - mHM run with a single L0 grid cell, also in the routing mode ! Stephan Thober, Jun 2014 - added flag for switching of MPR @@ -238,7 +238,7 @@ subroutine mHM( & unsat_thresh , & ! Threshold water depth in upper reservoir water_thresh_sealed , & ! Threshold water depth in impervious areas wilting_point ) ! Permanent wilting point for each horizon - + ! subroutines required to estimate variables prior to the MPR call use mo_upscaling_operators, only: L0_fractionalCover_in_Lx ! land cover fraction use mo_multi_param_reg, only: mpr,canopy_intercept_param ! reg. and scaling @@ -250,11 +250,11 @@ subroutine mHM( & use mo_soil_moisture, only: soil_moisture use mo_neutrons, only: DesiletsN0, COSMIC use mo_runoff, only: runoff_unsat_zone - use mo_runoff, only: runoff_sat_zone, L1_total_runoff + use mo_runoff, only: runoff_sat_zone, L1_total_runoff use mo_julian, only: dec2date, date2dec use mo_string_utils, only: num2str use mo_mhm_constants, only: HarSamConst ! parameters for Hargreaves-Samani Equation - + implicit none ! Intent @@ -263,8 +263,8 @@ subroutine mHM( & real(dp), intent(in) :: fSealedInCity ! fraction of perfectly sealed area within cities integer(i4), intent(in) :: timeStep_LAI_input ! time step of gridded LAI input integer(i4), intent(in) :: counter_year ! counter to tackle the change of year - integer(i4), intent(in) :: counter_month ! counter to tackle the change of month - integer(i4), intent(in) :: counter_day ! counter to tackle the change of day + integer(i4), intent(in) :: counter_month ! counter to tackle the change of month + integer(i4), intent(in) :: counter_day ! counter to tackle the change of day integer(i4), intent(in) :: tt real(dp), intent(in) :: time integer(i4), dimension(:,:), intent(in) :: processMatrix @@ -280,7 +280,7 @@ subroutine mHM( & integer(i4), intent(in) :: LCyearId integer(i4), dimension(:), intent(in) :: GeoUnitList integer(i4), dimension(:), intent(in) :: GeoUnitKar - integer(i4), dimension(:), intent(in) :: LAIUnitList + integer(i4), dimension(:), intent(in) :: LAIUnitList real(dp), dimension(:,:), intent(in) :: LAILUT ! Physiographic L0 @@ -398,8 +398,8 @@ subroutine mHM( & integer(i4) :: doy ! doy of the year [1-365 or 1-366] integer(i4) :: k ! cell index - real(dp) :: pet ! - real(dp) :: prec ! + real(dp) :: pet ! + real(dp) :: prec ! real(dp) :: temp ! ! temporary arrays so that inout of routines is contiguous array @@ -414,9 +414,9 @@ subroutine mHM( & !------------------------------------------------------------------- ! MPR CALL - ! -> call only when LC had changed - ! -> or for very first time step - ! + ! -> call only when LC had changed + ! -> or for very first time step + ! ! Variables required prior to the MPR call: ! 1) LC fractions ! --> time independent variable: to be initalized every time @@ -426,7 +426,7 @@ subroutine mHM( & ! with landcover change !------------------------------------------------------------------- if( (LCyearId .NE. yId) .or. (tt .EQ. 1) ) then - + ! abort if land cover change is there and mpr is switched off if ( (tt .ne. 1) .and. (.not. perform_mpr) ) then call message() @@ -434,9 +434,9 @@ subroutine mHM( & stop end if - ! update yId to keep track of LC change - yId = LCyearId - + ! update yId to keep track of LC change + yId = LCyearId + ! estimate land cover fractions for dominant landcover class ! --> time independent variable: to be initalized every time ! with landcover change @@ -459,7 +459,7 @@ subroutine mHM( & L0rightBound_inL1, & nTCells0_inL1 ) !--------------------------------------------------------- - ! Update fractions of sealed area fractions + ! Update fractions of sealed area fractions ! based on the sealing fraction[0-1] in cities !--------------------------------------------------------- fSealed1(:) = fSealedInCity*fSealed1(:) @@ -469,7 +469,7 @@ subroutine mHM( & fForest1(:) = fForest1(:) / ( fForest1(:) + fSealed1(:) + fPerm1(:) ) fSealed1(:) = fSealed1(:) / ( fForest1(:) + fSealed1(:) + fPerm1(:) ) fPerm1(:) = fPerm1(:) / ( fForest1(:) + fSealed1(:) + fPerm1(:) ) - + !------------------------------------------------------------------- ! NOW call MPR !------------------------------------------------------------------- @@ -490,8 +490,8 @@ subroutine mHM( & temp_thresh, unsat_thresh, water_thresh_sealed, wilting_point ) end if !------------------------------------------------------------------- - ! Update the inital states of soil water content for the first time - ! step and when perform_mpr = FALSE + ! Update the inital states of soil water content for the first time + ! step and when perform_mpr = FALSE ! based on the half of the derived values of Field capacity ! other states are kept at their inital values !------------------------------------------------------------------- @@ -500,7 +500,7 @@ subroutine mHM( & end if end if - + !------------------------------------------------------------------- ! CALL regionalization of parameters related to LAI ! IT is now outside of mHM since LAI is now dynamic variable @@ -509,18 +509,18 @@ subroutine mHM( & case(0) ! Estimate max. intecept. capacity based on long term monthly mean LAI values ! Max. interception is updated every month rather than every day - if( (tt .EQ. 1) .OR. (month .NE. counter_month) ) then + if( (tt .EQ. 1) .OR. (month .NE. counter_month) ) then call canopy_intercept_param( processMatrix, global_parameters(:), & - LAI0, nTCells0_inL1, L0upBound_inL1, & + LAI0, nTCells0_inL1, L0upBound_inL1, & L0downBound_inL1, L0leftBound_inL1, & L0rightBound_inL1, cellId0, mask0, & - nodata_dp, interc_max ) + nodata_dp, interc_max ) end if ! Estimate max. inteception based on daily LAI values case(-1) ! daily if ( (tt .EQ. 1) .OR. (day .NE. counter_day) ) then call canopy_intercept_param( processMatrix, global_parameters(:), & - LAI0, nTCells0_inL1, L0upBound_inL1, & + LAI0, nTCells0_inL1, L0upBound_inL1, & L0downBound_inL1, L0leftBound_inL1, & L0rightBound_inL1, cellId0, mask0, & nodata_dp, interc_max ) @@ -528,7 +528,7 @@ subroutine mHM( & case(-2) ! monthly if ( (tt .EQ. 1) .OR. (month .NE. counter_month) ) then call canopy_intercept_param( processMatrix, global_parameters(:), & - LAI0, nTCells0_inL1, L0upBound_inL1, & + LAI0, nTCells0_inL1, L0upBound_inL1, & L0downBound_inL1, L0leftBound_inL1, & L0rightBound_inL1, cellId0, mask0, & nodata_dp, interc_max ) @@ -536,7 +536,7 @@ subroutine mHM( & case(-3) ! yearly if ( (tt .EQ. 1) .OR. (year .NE. counter_year) ) then call canopy_intercept_param( processMatrix, global_parameters(:), & - LAI0, nTCells0_inL1, L0upBound_inL1, & + LAI0, nTCells0_inL1, L0upBound_inL1, & L0downBound_inL1, L0leftBound_inL1, & L0rightBound_inL1, cellId0, mask0, & nodata_dp, interc_max ) @@ -549,9 +549,9 @@ subroutine mHM( & ! flag for day or night depending on hours of the day !------------------------------------------------------------------- isday = ( hour .gt. 6 ) .AND. ( hour .le. 18 ) - + !------------------------------------------------------------------- - ! HYDROLOGICAL PROCESSES at L1-LEVEL + ! HYDROLOGICAL PROCESSES at L1-LEVEL !------------------------------------------------------------------- !$OMP parallel default(shared) & @@ -570,19 +570,19 @@ subroutine mHM( & ! if (tmax_in(k) .LE. tmin_in(k)) call message('WARNING: tmax smaller tmin at doy ', & num2str(doy), ' in year ', num2str(year),' at cell', num2str(k),'!') - + pet = fAsp(k) * pet_hargreaves(HarSamCoeff(k), HarSamConst, temp_in(k), tmax_in(k), & - tmin_in(k), latitude(k), doy) + tmin_in(k), latitude(k), doy) case(2) ! Priestley-Taylor ! Priestley Taylor is not defined for values netrad < 0.0_dp - pet = pet_priestly( PrieTayAlpha(k,month), max(netrad_in(k), 0.0_dp), temp_in(k)) + pet = pet_priestly( PrieTayAlpha(k,month), max(netrad_in(k), 0.0_dp), temp_in(k)) case(3) ! Penman-Monteith pet = pet_penman (max(netrad_in(k), 0.0_dp), temp_in(k), absvappres_in(k)/1000.0_dp, & aeroResist(k,month) / windspeed_in(k), surfResist(k,month)) end select - + ! temporal disaggreagtion of forcing variables call temporal_disagg_forcing( isday, ntimesteps_day, prec_in(k), & ! Intent IN pet, temp_in(k), fday_prec(month), fday_pet(month), & ! Intent IN @@ -615,7 +615,7 @@ subroutine mHM( & infiltration(k, nHorizons_mHM), unsat_thresh(k), & ! Intent IN satStorage(k), unsatStorage(k), & ! Intent INOUT slow_interflow(k), fast_interflow(k), perc(k) ) ! Intent OUT - + call runoff_sat_zone( k2(k), & ! Intent IN satStorage(k), & ! Intent INOUT baseflow(k) ) ! Intent OUT @@ -631,8 +631,8 @@ subroutine mHM( & !------------------------------------------------------------------- ! Nested model: Neutrons state variable, related to soil moisture !------------------------------------------------------------------- - - ! based on soilMoisture + + ! based on soilMoisture ! TODO they again loop over all cells. Maybe move this to line 680 in the loop used above? if ( processMatrix(10, 1) .eq. 1 ) & call DesiletsN0( soilMoisture(:,:), horizon_depth(:), & @@ -642,16 +642,7 @@ subroutine mHM( & call COSMIC( soilMoisture(:,:), horizon_depth(:), & global_parameters(processMatrix(10,3)-processMatrix(10,2)+2:processMatrix(10,3)), & neutrons(:)) - + end subroutine mHM END MODULE mo_mHM - - - - - - - - - diff --git a/src/mHM/mo_mhm_constants.f90 b/src/mHM/mo_mhm_constants.f90 index ca6d26f9..7692db70 100644 --- a/src/mHM/mo_mhm_constants.f90 +++ b/src/mHM/mo_mhm_constants.f90 @@ -14,8 +14,8 @@ MODULE mo_mhm_constants IMPLICIT NONE PRIVATE - - ! natural + + ! natural real(dp), public, parameter :: H2Odens = 1000.0_dp ! Density of water (kg/m3) ! computational @@ -28,7 +28,7 @@ MODULE mo_mhm_constants integer(i4), public, parameter :: nColPars = 5_i4 ! number of properties of the global variables integer(i4), public, parameter :: maxNoSoilHorizons = 10_i4 ! maximum number of allowed soil layers integer(i4), public, parameter :: maxNoBasins = 50_i4 ! maximum number of allowed basins - integer(i4), public, parameter :: maxNLcovers = 50_i4 ! maximum number of allowed LCover scenes + integer(i4), public, parameter :: maxNLcovers = 50_i4 ! maximum number of allowed LCover scenes integer(i4), public, parameter :: maxGeoUnit = 25_i4 ! maximum number of allowed geological classes ! default inital values for states and fluxes as well as parameter fields @@ -48,16 +48,16 @@ MODULE mo_mhm_constants integer(i4), public, parameter :: YearMonths_i4 = 12 ! months per year real(dp), public, parameter :: YearDays = 365.0_dp ! days in a year real(dp), public, parameter :: DaySecs = 86400.0_dp ! sec in a day - + ! soil paramterization (mo_mpr_soilmoist) - ! organic matter constant for calculation of mineral bulk density following RAWL + ! organic matter constant for calculation of mineral bulk density following RAWL real(dp), public, parameter :: BulkDens_OrgMatter = 0.224_dp ! [g/cm3] from W.R. RAWLS ! constants for determinination of the field capacity following Twarakavi real(dp), public, parameter :: field_cap_c1 = -0.60_dp ! field capacity constant 1 real(dp), public, parameter :: field_cap_c2 = 2.0_dp ! field capacity constant 2 ! constants for determinination of the van Genuchten parameter n and sand treshold real(dp), public, parameter :: vGenuchten_sandtresh = 66.5_dp ! van Genuchten snad treshold - real(dp), public, parameter :: vGenuchtenN_c1 = 1.392_dp ! constants for van Genuchten n + real(dp), public, parameter :: vGenuchtenN_c1 = 1.392_dp ! constants for van Genuchten n real(dp), public, parameter :: vGenuchtenN_c2 = 0.418_dp real(dp), public, parameter :: vGenuchtenN_c3 = -0.024_dp real(dp), public, parameter :: vGenuchtenN_c4 = 1.212_dp @@ -82,13 +82,13 @@ MODULE mo_mhm_constants real(dp), public, parameter :: PWP_matPot_ThetaR = 15000.0_dp ! [hPa] matrix potential of -1500 kPa, assumed as thetaR=0 - !> Stefan-Boltzmann constant [W m^-2 K^-4] - real(dp), public, parameter :: StBoltzmann = 5.67e-08_dp + !> Stefan-Boltzmann constant [W m^-2 K^-4] + real(dp), public, parameter :: StBoltzmann = 5.67e-08_dp !> Hargreaves-Samani ref. ET formula [deg C] - real(dp), public, parameter :: HarSamConst = 17.800_dp + real(dp), public, parameter :: HarSamConst = 17.800_dp !> assumed meteorol. measurement hight for estimation of aeroResist and surfResist - real(dp), public, parameter :: WindMeasHeight = 10.0_dp - !> von karman constant + real(dp), public, parameter :: WindMeasHeight = 10.0_dp + !> von karman constant real(dp), public, parameter :: karman = 0.41_dp !> LAI factor for bulk surface resistance formulation @@ -110,16 +110,16 @@ MODULE mo_mhm_constants real(dp), public, parameter :: tetens_c3 = 237.30_dp !> calculation of the slope of the saturation vapour pressure curve following Tetens real(dp), public, parameter :: satpressureslope1 = 4098.0_dp - + !> Neutrons and moisture: N0 formula, Desilets et al. 2010 real(dp), public, parameter :: Desilets_a0 = 0.0808_dp real(dp), public, parameter :: Desilets_a1 = 0.372_dp real(dp), public, parameter :: Desilets_a2 = 0.115_dp - + !> Neutrons and moisture: COSMIC, Shuttleworth et al. 2013 real(dp), public, parameter :: COSMIC_bd = 1.4020_dp ! Dry soil bulk density (g/m3) real(dp), public, parameter :: COSMIC_vwclat = 0.0753_dp ! Volumetric "lattice" water content (m3/m3) - real(dp), public, parameter :: COSMIC_N = 510.51737902_dp ! High energy neutron flux (-) + real(dp), public, parameter :: COSMIC_N = 348.33_dp ! High energy neutron flux (cph), original was 510.51737902_dp real(dp), public, parameter :: COSMIC_alpha = 0.2392421548_dp ! Ratio of Fast Neutron Creation Factor (Soil to Water) real(dp), public, parameter :: COSMIC_L1 = 161.98621864_dp ! High Energy Soil Attenuation Length (g/cm2) real(dp), public, parameter :: COSMIC_L2 = 129.14558985_dp ! High Energy Water Attenuation Length (g/cm2) diff --git a/src/mHM/mo_mhm_eval.f90 b/src/mHM/mo_mhm_eval.f90 index 776ff6fd..0aabb204 100644 --- a/src/mHM/mo_mhm_eval.f90 +++ b/src/mHM/mo_mhm_eval.f90 @@ -15,7 +15,7 @@ MODULE mo_mhm_eval PRIVATE - PUBLIC :: mhm_eval + PUBLIC :: mhm_eval ! ------------------------------------------------------------------ @@ -31,7 +31,7 @@ MODULE mo_mhm_eval !> \details Runs mhm with a specific parameter set and returns required variables, e.g. runoff. ! INTENT(IN) - !> \param[in] "real(dp), dimension(:) :: parameterset" + !> \param[in] "real(dp), dimension(:) :: parameterset" !> a set of global parameter (gamma) to run mHM, DIMENSION [no. of global_Parameters] ! INTENT(INOUT) @@ -47,11 +47,11 @@ MODULE mo_mhm_eval ! None ! INTENT(OUT), OPTIONAL - !> \param[out] "real(dp), dimension(:,:), optional :: runoff" + !> \param[out] "real(dp), dimension(:,:), optional :: runoff" !> returns runoff time series, DIMENSION [nTimeSteps, nGaugesTotal] - !> \param[out] "real(dp), dimension(:,:), optional :: sm_opti" + !> \param[out] "real(dp), dimension(:,:), optional :: sm_opti" !> returns soil moisture time series for all grid cells (of multiple basins concatenated), DIMENSION [nCells, nTimeSteps] - !> \param[out] "real(dp), dimension(:,:), optional :: basin_avg_tws" + !> \param[out] "real(dp), dimension(:,:), optional :: basin_avg_tws" !> returns basin averaged total water storage time series, DIMENSION [nTimeSteps, nBasins] ! RETURN @@ -70,7 +70,7 @@ MODULE mo_mhm_eval !> \author Juliane Mai, Rohini Kumar !> \date Feb 2013 ! Modified, R. Kumar, Jun 2013 - restart_flag_states_read is passed to mhm call - ! for the soil moisture initalisation + ! for the soil moisture initalisation ! R. Kumar, Jun 2013 - frac_sealed_city_area is added ! R. Kumar & S. Thober, Aug 2013 - code change to incorporate output timestep ! during writing of the netcdf file @@ -91,7 +91,7 @@ MODULE mo_mhm_eval ! Stephan Thober, Sep 2015 - updated mrm_routing call ! Oldrich Rakovec, Rohini Kumar, Oct 2015 - added optional output for basin averaged TWS - SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) + SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws, neutrons_opti) use mo_init_states, only : get_basin_info use mo_init_states, only : variables_default_init ! default initalization of variables @@ -108,46 +108,47 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) read_restart, perform_mpr, fracSealed_CityArea, & timeStep_model_inputs, & timeStep, nBasins, simPer, readPer, & ! [h] simulation time step, No. of basins - processMatrix, c2TSTu, HorizonDepth_mHM, & - nSoilHorizons_mHM, NTSTEPDAY, timeStep, & - LCyearId, LAIUnitList, LAILUT, & + processMatrix, c2TSTu, HorizonDepth_mHM, & + nSoilHorizons_mHM, NTSTEPDAY, timeStep, & + LCyearId, LAIUnitList, LAILUT, & GeoUnitList, GeoUnitKar, soilDB, & - L0_Id, L0_soilId, & + L0_Id, L0_soilId, & L0_LCover, L0_asp, L0_LCover_LAI, L0_geoUnit, & - soilDB, L1_nTCells_L0, & + soilDB, L1_nTCells_L0, & L0_slope_emp, & - L1_upBound_L0, L1_downBound_L0, L1_leftBound_L0, & + L1_upBound_L0, L1_downBound_L0, L1_leftBound_L0, & L1_rightBound_L0, L0_latitude, L1_latitude, & - evap_coeff, fday_prec, & - fnight_prec, fday_pet, fnight_pet, fday_temp, & + evap_coeff, fday_prec, & + fnight_prec, fday_pet, fnight_pet, fday_temp, & fnight_temp, L1_pet, L1_tmin, L1_tmax, L1_netrad, & L1_absvappress, L1_windspeed, & - L1_pre, L1_temp , L1_fForest, & - L1_fPerm, L1_fSealed, L1_inter, & - L1_snowPack, L1_sealSTW, L1_soilMoist, L1_unsatSTW, & + L1_pre, L1_temp , L1_fForest, & + L1_fPerm, L1_fSealed, L1_inter, & + L1_snowPack, L1_sealSTW, L1_soilMoist, L1_unsatSTW, & L1_satSTW, L1_pet_calc, & L1_aETSoil, L1_aETCanopy, L1_aETSealed, & - L1_baseflow, L1_infilSoil, L1_fastRunoff, L1_melt, & - L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, & - L1_slowRunoff, L1_snow, L1_Throughfall, & + L1_baseflow, L1_infilSoil, L1_fastRunoff, L1_melt, & + L1_percol, L1_preEffect, L1_rain, L1_runoffSeal, & + L1_slowRunoff, L1_snow, L1_Throughfall, & L1_total_runoff, L1_alpha, L1_degDayInc, & - L1_degDayMax, & - L1_degDayNoPre, L1_degDay, L1_fAsp, L1_HarSamCoeff, & + L1_degDayMax, & + L1_degDayNoPre, L1_degDay, L1_fAsp, L1_HarSamCoeff, & L1_PrieTayAlpha, L1_aeroResist, L1_surfResist, & - L1_fRoots, L1_maxInter, L1_karstLoss, L1_kfastFlow, & - L1_kSlowFlow, L1_kBaseFlow, L1_kPerco, & - L1_soilMoistFC, L1_soilMoistSat, L1_soilMoistExp, & - L1_tempThresh, L1_unsatThresh, L1_sealedThresh, & + L1_fRoots, L1_maxInter, L1_karstLoss, L1_kfastFlow, & + L1_kSlowFlow, L1_kBaseFlow, L1_kPerco, & + L1_soilMoistFC, L1_soilMoistSat, L1_soilMoistExp, & + L1_tempThresh, L1_unsatThresh, L1_sealedThresh, & L1_wiltingPoint, L1_neutrons, & basin_avg_TWS_sim, & warmingDays, & timeStep_LAI_input, & ! flag on how LAI data has to be read L0_gridded_LAI, dirRestartIn, & ! restart directory location timeStep_sm_input, & ! time step of soil moisture input (day, month, year) - nSoilHorizons_sm_input, & ! no. of mhm soil horizons equivalent to sm input - nTimeSteps_L1_sm ! total number of timesteps in soil moisture input + nSoilHorizons_sm_input, & ! no. of mhm soil horizons equivalent to sm input + nTimeSteps_L1_sm, & ! total number of timesteps in soil moisture input + nTimeSteps_L1_neutrons ! total number of timesteps in neutrons input use mo_common_variables, only: & - optimize + optimize #ifdef mrm2mhm use mo_utils, only: ge use mo_mrm_global_variables, only: & @@ -188,18 +189,19 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) use mo_mrm_write, only: mrm_write_output_fluxes use mo_mrm_init, only: variables_default_init_routing #endif - + implicit none real(dp), dimension(:), intent(in) :: parameterset - real(dp), dimension(:,:), allocatable, optional, intent(out) :: runoff ! dim1=time dim2=gauge - real(dp), dimension(:,:), allocatable, optional, intent(out) :: sm_opti ! dim1=ncells, dim2=time - real(dp), dimension(:,:), allocatable, optional, intent(out) :: basin_avg_tws! dim1=time dim2=nBasins + real(dp), dimension(:,:), allocatable, optional, intent(out) :: runoff ! dim1=time dim2=gauge + real(dp), dimension(:,:), allocatable, optional, intent(out) :: sm_opti ! dim1=ncells, dim2=time + real(dp), dimension(:,:), allocatable, optional, intent(out) :: basin_avg_tws ! dim1=time dim2=nBasins + real(dp), dimension(:,:), allocatable, optional, intent(out) :: neutrons_opti ! dim1=ncells, dim2=time ! ------------------------------------- ! local variables ! - ! FOR WRITING GRIDDED STATES AND FLUXES + ! FOR WRITING GRIDDED STATES AND FLUXES integer(i4) :: tIndex_out ! for writing netcdf file real(dp), dimension(size(L1_fSealed)) :: L1_fNotSealed type(OutputDataset) :: nc @@ -233,25 +235,25 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) logical :: writeout ! if true write out netcdf files integer(i4) :: writeout_counter ! write out time step ! -#ifdef mrm2mhm +#ifdef mrm2mhm ! for routing logical :: do_mpr integer(i4) :: s11, e11 ! start and end index at L11 integer(i4) :: s110, e110 ! start and end index of L11 at L0 logical, allocatable :: mask11(:,:) -#endif +#endif ! ! for basin average tws timeseries integer(i4) :: gg real(dp), dimension(:), allocatable :: TWS_field ! field of TWS real(dp) :: area_basin - + ! LAI options integer(i4) :: day_counter integer(i4) :: month_counter real(dp), dimension(:), allocatable :: LAI ! local variable for leaf area index - + !---------------------------------------------------------- ! Check optionals and initialize !---------------------------------------------------------- @@ -269,19 +271,27 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) allocate(sm_opti(size(L1_pre, dim=1), nTimeSteps_L1_sm)) sm_opti(:,:) = 0.0_dp ! has to be intialized with zero because later summation end if + ! neutrons optimization + !-------------------------- + if ( present(neutrons_opti) ) then + ! ! total No of cells, No of timesteps + ! ! of all basins , in neutrons input + allocate(neutrons_opti(size(L1_pre, dim=1), nTimeSteps_L1_neutrons)) + neutrons_opti(:,:) = 0.0_dp ! has to be intialized with zero because later summation + end if ! add other optionals... - + !------------------------------------------------------------------- - ! Initalize State variables either to the default value or + ! Initalize State variables either to the default value or ! from the restrat_files. - ! All variables to be initalized had been allocated to the required + ! All variables to be initalized had been allocated to the required ! space before this point (see, mo_startup: initialise) !------------------------------------------------------------------- if (.NOT. read_restart ) then - ! as default values, + ! as default values, ! all cells for all modeled basins are simultenously initalized ONLY ONCE call variables_default_init() #ifdef mrm2mhm @@ -292,8 +302,8 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) !------------------------------------------- call variables_default_init_routing() end if -#endif - else +#endif + else ! read from restart files, basin wise ... do ii = 1, nBasins call read_restart_states(ii, dirRestartIn(ii) ) @@ -308,28 +318,28 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) ! calculate NtimeSteps for this basin nTimeSteps = ( simPer(ii)%julEnd - simPer(ii)%julStart + 1 ) * NTSTEPDAY - + ! reinitialize time counter for LCover and MPR - ! -0.5 is due to the fact that dec2date routine + ! -0.5 is due to the fact that dec2date routine ! changes the day at 12:00 in NOON ! Whereas mHM needs day change at 00:00 h newTime = real(simPer(ii)%julStart,dp) yId = 0 ! get basin information - call get_basin_info ( ii, 0, nrows, ncols, iStart=s0, iEnd=e0, mask=mask0 ) - call get_basin_info ( ii, 1, nrows, ncols, ncells=nCells, iStart=s1, iEnd=e1, mask=mask1 ) + call get_basin_info ( ii, 0, nrows, ncols, iStart=s0, iEnd=e0, mask=mask0 ) + call get_basin_info ( ii, 1, nrows, ncols, ncells=nCells, iStart=s1, iEnd=e1, mask=mask1 ) -#ifdef mrm2mhm +#ifdef mrm2mhm if (processMatrix(8,1) .eq. 1) then ! read states from restart if (read_restart) call mrm_read_restart_states(ii, dirRestartIn(ii)) ! ! get basin information at L11 and L110 if routing is activated - call get_basin_info_mrm ( ii, 11, nrows, ncols, iStart=s11, iEnd=e11, mask=mask11 ) - call get_basin_info_mrm ( ii, 110, nrows, ncols, iStart=s110, iEnd=e110 ) + call get_basin_info_mrm ( ii, 11, nrows, ncols, iStart=s11, iEnd=e11, mask=mask11 ) + call get_basin_info_mrm ( ii, 110, nrows, ncols, iStart=s110, iEnd=e110 ) end if -#endif +#endif ! allocate space for local LAI grid allocate( LAI(s0:e0) ) LAI(:) = nodata_dp @@ -339,7 +349,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) allocate(TWS_field(s1:e1)) TWS_field(s1:e1) = nodata_dp end if - + ! Loop over time average_counter = 0 writeout_counter = 0 @@ -354,16 +364,16 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) s_meteo = s1 e_meteo = e1 ! time step for meteorological variable (daily values) - iMeteoTS = ceiling( real(tt,dp) / real(NTSTEPDAY,dp) ) + iMeteoTS = ceiling( real(tt,dp) / real(NTSTEPDAY,dp) ) else - ! read chunk of meteorological forcings data (reading, upscaling/downscaling) + ! read chunk of meteorological forcings data (reading, upscaling/downscaling) call prepare_meteo_forcings_data(ii, tt) ! set start and end of meteo position s_meteo = 1 e_meteo = e1 - s1 + 1 ! time step for meteorological variable (daily values) iMeteoTS = ceiling( real(tt,dp) / real(NTSTEPDAY,dp) ) & - - ( readPer%julStart - simPer(ii)%julStart ) + - ( readPer%julStart - simPer(ii)%julStart ) end if hour = mod(hour+timestep, 24) @@ -392,7 +402,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) ! customize iMeteoTS for process 5 - PET select case (processMatrix(5,1)) - ! (/ pet, tmin, tmax, netrad, absVapP,windspeed /) + ! (/ pet, tmin, tmax, netrad, absVapP,windspeed /) case(0) ! PET is input iMeteo_p5 = (/iMeteoTS, 1, 1, 1, 1, 1 /) case(1) ! HarSam @@ -402,7 +412,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) case(3) ! PenMon iMeteo_p5 = (/iMeteoTS, 1, 1, iMeteoTS, iMeteoTS, iMeteoTS /) end select - + !-------------------------------------------------------------------- ! call LAI function to get LAI fields for this timestep and basin !-------------------------------------------------------------------- @@ -412,18 +422,18 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) month_counter = month year_counter = year end if - ! + ! select case(timeStep_LAI_input) case(0) ! create gridded fields of LAI using long term monthly mean LAI values - ! and the corresponding LC file + ! and the corresponding LC file ! update LAI --> for 1st timestep and when month changes if( (tt .EQ. 1) .OR. (month_counter .NE. month) ) then do ll = 1, size(LAIUnitList) where( L0_LCover_LAI(s0:e0) .EQ. LAIUnitList(ll) ) LAI(:) = LAILUT(ll, month) end do end if - ! + ! case(-1) ! daily if ( (tt .EQ. 1) .OR. (day .NE. day_counter) ) then iGridLAI_TS = iGridLAI_TS + 1_i4 @@ -458,9 +468,9 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) ! X FLUXES (L1, L11 levels) ! -------------------------------------------------------------------------- call mhm(perform_mpr, read_restart, fracSealed_cityArea, & ! IN C - timeStep_LAI_input, year_counter, month_counter, day_counter, & ! IN C + timeStep_LAI_input, year_counter, month_counter, day_counter, & ! IN C tt, newTime-0.5_dp, processMatrix, c2TSTu, HorizonDepth_mHM, & ! IN C - nCells, nSoilHorizons_mHM, real(NTSTEPDAY,dp), mask0, & ! IN C + nCells, nSoilHorizons_mHM, real(NTSTEPDAY,dp), mask0, & ! IN C parameterset, & ! IN P LCyearId(year,ii), GeoUnitList, GeoUnitKar, LAIUnitList, LAILUT, & ! IN L0 L0_slope_emp(s0:e0), L0_Latitude(s0:e0), & ! IN L0 @@ -481,12 +491,13 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) L1_netrad(s_p5(4):e_p5(4), iMeteo_p5(4)), & ! IN F:PET L1_absvappress(s_p5(5):e_p5(5), iMeteo_p5(5)), & ! IN F:PET L1_windspeed(s_p5(6):e_p5(6), iMeteo_p5(6)), & ! IN F:PET - L1_pre(s_meteo:e_meteo,iMeteoTS), & ! IN F:Pre + L1_pre(s_meteo:e_meteo,iMeteoTS), & ! IN F:Pre L1_temp(s_meteo:e_meteo,iMeteoTS), & ! IN F:Temp yId, & ! INOUT C - L1_fForest(s1:e1), L1_fPerm(s1:e1), L1_fSealed(s1:e1), & ! INOUT L1 - L1_inter(s1:e1), L1_snowPack(s1:e1), L1_sealSTW(s1:e1), & ! INOUT S - L1_soilMoist(s1:e1,:), L1_unsatSTW(s1:e1), L1_satSTW(s1:e1), L1_neutrons, & ! INOUT S + L1_fForest(s1:e1), L1_fPerm(s1:e1), L1_fSealed(s1:e1), & ! INOUT L1 + L1_inter(s1:e1), L1_snowPack(s1:e1), L1_sealSTW(s1:e1), & ! INOUT S + L1_soilMoist(s1:e1,:), L1_unsatSTW(s1:e1), L1_satSTW(s1:e1), & ! INOUT S + L1_neutrons(s1:e1), & ! INOUT S L1_pet_calc(s1:e1), & ! INOUT X L1_aETSoil(s1:e1,:), L1_aETCanopy(s1:e1), L1_aETSealed(s1:e1), & ! INOUT X L1_baseflow(s1:e1), L1_infilSoil(s1:e1,:), L1_fastRunoff(s1:e1), & ! INOUT X @@ -560,7 +571,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) if (day_counter .NE. day ) day_counter = day if (month_counter .NE. month) month_counter = month if (year_counter .NE. year) year_counter = year - + ! increment of timestep newTime = julday(day,month,year) + real(hour+timestep,dp)/24._dp ! calculate new year, month and day @@ -586,7 +597,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) L11_qmod(s11:e11)) end if #endif - + ! output only for evaluation period tIndex_out = (tt-warmingDays(ii)*NTSTEPDAY) ! tt if write out of warming period @@ -598,7 +609,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) L1_fNotSealed = 1.0_dp - L1_fSealed nc = OutputDataset(ii, mask1) end if - + call nc%updateDataset( & s1 , & e1 , & @@ -625,7 +636,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) L1_infilSoil , & L1_preEffect & ) - + ! write data writeout = .false. if (timeStep_model_outputs .gt. 0) then @@ -705,21 +716,47 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) if( present(basin_avg_tws) ) then area_basin = sum(L1_areaCell(s1:e1) ) TWS_field(s1:e1) = L1_inter(s1:e1) + L1_snowPack(s1:e1) + L1_sealSTW(s1:e1) + & - L1_unsatSTW(s1:e1) + L1_satSTW(s1:e1) + L1_unsatSTW(s1:e1) + L1_satSTW(s1:e1) do gg = 1, nSoilHorizons_mHM - TWS_field(s1:e1) = TWS_field(s1:e1) + L1_soilMoist (s1:e1,gg) + TWS_field(s1:e1) = TWS_field(s1:e1) + L1_soilMoist (s1:e1,gg) end do - basin_avg_TWS_sim(tt,ii) = ( dot_product( TWS_field (s1:e1), L1_areaCell(s1:e1) ) / area_basin ) + basin_avg_TWS_sim(tt,ii) = ( dot_product( TWS_field (s1:e1), L1_areaCell(s1:e1) ) / area_basin ) end if !---------------------------------------------------------------------- - + + !---------------------------------------------------------------------- + ! FOR NEUTRONS + ! NOTE:: modeled neutrons are averaged daily + !---------------------------------------------------------------------- + if (present(neutrons_opti)) then + if ( tt .EQ. 1 ) writeout_counter = 1 + ! only for evaluation period - ignore warming days + if ( (tt-warmingDays(ii)*NTSTEPDAY) .GT. 0 ) then + ! decide for daily, monthly or yearly aggregation + ! daily + if (day .NE. day_counter) then + neutrons_opti(s1:e1,writeout_counter) = neutrons_opti(s1:e1,writeout_counter) / real(average_counter,dp) + writeout_counter = writeout_counter + 1 + average_counter = 0 + end if + + ! last timestep is already done - write_counter exceeds size(sm_opti, dim=2) + if (.not. (tt .eq. nTimeSteps) ) then + ! aggregate neutrons to needed time step for optimization + neutrons_opti(s1:e1,writeout_counter) = neutrons_opti(s1:e1,writeout_counter) + L1_neutrons(s1:e1) + end if + + average_counter = average_counter + 1 + end if + end if + end do !<< TIME STEPS LOOP ! deallocate space for temprory LAI fields deallocate(LAI) ! deallocate TWS field temporal variable if (allocated (TWS_field) ) deallocate(TWS_field) - + end do !<< BASIN LOOP #ifdef mrm2mhm @@ -732,7 +769,7 @@ SUBROUTINE mhm_eval(parameterset, runoff, sm_opti, basin_avg_tws) ! ========================================================================= ! SET TWS OUTPUT VARIABLE ! ========================================================================= - if( present(basin_avg_tws) ) basin_avg_tws = basin_avg_TWS_sim + if( present(basin_avg_tws) ) basin_avg_tws = basin_avg_TWS_sim end SUBROUTINE mhm_eval diff --git a/src/mHM/mo_objective_function.f90 b/src/mHM/mo_objective_function.f90 index 12d19a33..5f3c5792 100644 --- a/src/mHM/mo_objective_function.f90 +++ b/src/mHM/mo_objective_function.f90 @@ -45,7 +45,7 @@ MODULE mo_objective_function !> \brief Wrapper for objective functions. - !> \details The functions selects the objective function case defined in a namelist, + !> \details The functions selects the objective function case defined in a namelist, !> i.e. the global variable \e opti\_function.\n !> It return the objective function value for a specific parameter set. @@ -68,7 +68,7 @@ MODULE mo_objective_function ! None ! RETURN - !> \return real(dp) :: objective — objective function value + !> \return real(dp) :: objective — objective function value !> (which will be e.g. minimized by an optimization routine like DDS) ! RESTRICTIONS @@ -113,10 +113,13 @@ FUNCTION objective(parameterset) case (15) ! KGE for Q * RMSE for basin_avg TWS (standarized scored) objective = objective_kge_q_rmse_tws(parameterset) + case (17) + ! KGE of catchment average SM + objective = objective_neutrons_kge_catchment_avg(parameterset) case default stop "Error objective: opti_function not implemented yet." end select - + END FUNCTION objective ! ------------------------------------------------------------------ @@ -126,8 +129,8 @@ END FUNCTION objective !> \brief Objective function for soil moisture. - !> \details The objective function only depends on a parameter vector. - !> The model will be called with that parameter vector and + !> \details The objective function only depends on a parameter vector. + !> The model will be called with that parameter vector and !> the model output is subsequently compared to observed data.\n !> !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average @@ -139,13 +142,13 @@ END FUNCTION objective !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation \n !> is calculated and the objective function for a given basin \f$ i \f$ is !> \f[ \phi_{i} = 1.0 - KGE_{i} \f] - !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods. + !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods. !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.\n !> - !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 - !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. - !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f] \n - !> The observed data L1_sm, L1_sm_mask are global in this module. + !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 + !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. + !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f] \n + !> The observed data L1_sm, L1_sm_mask are global in this module. ! INTENT(IN) !> \param[in] "real(dp) :: parameterset(:)" 1D-array with parameters the model is run with @@ -166,7 +169,7 @@ END FUNCTION objective ! None ! RETURN - !> \return real(dp) :: objective_sm_kge_catchment_avg — objective function value + !> \return real(dp) :: objective_sm_kge_catchment_avg — objective function value !> (which will be e.g. minimized by an optimization routine like DDS) ! RESTRICTIONS @@ -184,7 +187,7 @@ END FUNCTION objective !> \date May 2015 FUNCTION objective_sm_kge_catchment_avg(parameterset) - + use mo_init_states, only : get_basin_info use mo_mhm_eval, only : mhm_eval use mo_message, only : message @@ -195,7 +198,7 @@ FUNCTION objective_sm_kge_catchment_avg(parameterset) use mo_global_variables, only: nBasins, & ! number of basins L1_sm, L1_sm_mask ! packed measured sm, sm-mask (dim1=ncells, dim2=time) use mo_mhm_constants, only: nodata_dp ! global nodata value - + implicit none real(dp), dimension(:), intent(in) :: parameterset @@ -223,7 +226,7 @@ FUNCTION objective_sm_kge_catchment_avg(parameterset) do iBasin=1, nBasins ! get basin information - call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) ! allocate allocate(mask_times (size(sm_opti, dim=2))) @@ -242,7 +245,7 @@ FUNCTION objective_sm_kge_catchment_avg(parameterset) if ( all(.NOT. L1_sm_mask(:,iTime)) .OR. (count(L1_sm_mask(:,iTime)) .LE. 10) ) then call message('WARNING: objective_sm_kge_catchment_avg: ignored currrent time step since less than') call message(' 10 valid cells available in soil moisture observation') - mask_times(iTime) = .FALSE. + mask_times(iTime) = .FALSE. cycle end if sm_catch_avg_basin(iTime) = average( L1_sm(s1:e1,iTime), mask=L1_sm_mask(s1:e1,iTime)) @@ -256,9 +259,9 @@ FUNCTION objective_sm_kge_catchment_avg(parameterset) end do objective_sm_kge_catchment_avg = objective_sm_kge_catchment_avg**onesixth - + call message(' objective_sm_kge_catchment_avg = ', num2str(objective_sm_kge_catchment_avg,'(F9.5)')) - + END FUNCTION objective_sm_kge_catchment_avg ! ------------------------------------------------------------------ @@ -268,11 +271,11 @@ END FUNCTION objective_sm_kge_catchment_avg !> \brief Objective function for soil moisture. - !> \details The objective function only depends on a parameter vector. - !> The model will be called with that parameter vector and + !> \details The objective function only depends on a parameter vector. + !> The model will be called with that parameter vector and !> the model output is subsequently compared to observed data.\n !> - !> Therefore the Pearson correlation between observed and modeled soil + !> Therefore the Pearson correlation between observed and modeled soil !> moisture on each grid cell \f$ j \f$ is compared !> \f[ r_j = r^2(SM_{obs}^j, SM_{sim}^j) \f] !> where \n @@ -281,17 +284,17 @@ END FUNCTION objective_sm_kge_catchment_avg !> \f$ SM_{sim} \f$ = simulated soil moisture.\n !> The observed data \f$ SM_{obs} \f$ are global in this module.\n !> - !> The the correlation is spatially averaged as - !> \f[ \phi_{i} = \frac{1}{K} \cdot \sum_{j=1}^K r_j \f] + !> The the correlation is spatially averaged as + !> \f[ \phi_{i} = \frac{1}{K} \cdot \sum_{j=1}^K r_j \f] ! !> where \f$ K \f$ denotes the number of valid cells in the study domain.\n ! - !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 - !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. - !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 }. \f] \n - !> The observed data L1_sm, L1_sm_mask are global in this module. + !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 + !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. + !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 }. \f] \n + !> The observed data L1_sm, L1_sm_mask are global in this module. + - ! INTENT(IN) !> \param[in] "real(dp) :: parameterset(:)" 1D-array with parameters the model is run with @@ -311,7 +314,7 @@ END FUNCTION objective_sm_kge_catchment_avg ! None ! RETURN - !> \return real(dp) :: objective_sm_corr — objective function value + !> \return real(dp) :: objective_sm_corr — objective function value !> (which will be e.g. minimized by an optimization routine like DDS) ! RESTRICTIONS @@ -329,7 +332,7 @@ END FUNCTION objective_sm_kge_catchment_avg !> \date March 2015 FUNCTION objective_sm_corr(parameterset) - + use mo_init_states, only : get_basin_info use mo_mhm_eval, only : mhm_eval use mo_message, only : message @@ -355,7 +358,7 @@ FUNCTION objective_sm_corr(parameterset) real(dp), parameter :: onesixth = 1.0_dp/6.0_dp real(dp), dimension(:,:), allocatable :: sm_opti ! simulated soil moisture ! ! (dim1=ncells, dim2=time) - + call mhm_eval(parameterset, sm_opti=sm_opti) ! initialize some variables @@ -364,10 +367,10 @@ FUNCTION objective_sm_corr(parameterset) ! loop over basin - for applying power law later on do iBasin=1, nBasins - ! init + ! init objective_sm_corr_basin = 0.0_dp ! get basin information - call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) ! temporal correlation is calculated on individual gridd cells do iCell = s1, e1 @@ -389,9 +392,9 @@ FUNCTION objective_sm_corr(parameterset) end do objective_sm_corr = objective_sm_corr**onesixth - + call message(' objective_sm_corr = ', num2str(objective_sm_corr,'(F9.5)')) - + END FUNCTION objective_sm_corr ! ------------------------------------------------------------------ @@ -401,11 +404,11 @@ END FUNCTION objective_sm_corr !> \brief Objective function for soil moisture. - !> \details The objective function only depends on a parameter vector. - !> The model will be called with that parameter vector and + !> \details The objective function only depends on a parameter vector. + !> The model will be called with that parameter vector and !> the model output is subsequently compared to observed data.\n !> - !> Therefore the Pattern Dissimilarity (PD) of observed and modeled soil + !> Therefore the Pattern Dissimilarity (PD) of observed and modeled soil !> moisture fields is calculated - aim: matching spatial patters !> \f[ E(t) = PD\left( SM_{obs}(t), SM_{sim}(t) \right) \f] !> where \n @@ -414,16 +417,16 @@ END FUNCTION objective_sm_corr !> \f$ SM_{sim} \f$ = simulated soil moisture.\n !> \f$ E(t) \f$ = pattern dissimilarity at timestep \f$ t \f$.\n ! - !> The the pattern dissimilaity (E) is spatially averaged as - !> \f[ \phi_{i} = \frac{1}{T} \cdot \sum_{t=1}^T E_t \f] + !> The the pattern dissimilaity (E) is spatially averaged as + !> \f[ \phi_{i} = \frac{1}{T} \cdot \sum_{t=1}^T E_t \f] ! !> where \f$ T \f$ denotes the number of time steps.\n ! - !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 - !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. - !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 } . \f] \n + !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 + !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. + !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 } . \f] \n !> The observed data L1_sm, L1_sm_mask are global in this module. - + ! !> The observed data L1_sm, L1_sm_mask are global in this module.\n @@ -446,7 +449,7 @@ END FUNCTION objective_sm_corr ! None ! RETURN - !> \return real(dp) :: objecive_sm_pd — objective function value + !> \return real(dp) :: objecive_sm_pd — objective function value !> (which will be e.g. minimized by an optimization routine like DDS) ! RESTRICTIONS @@ -464,7 +467,7 @@ END FUNCTION objective_sm_corr !> \date May 2015 FUNCTION objective_sm_pd(parameterset) - + use mo_init_states, only : get_basin_info use mo_mhm_eval, only : mhm_eval use mo_message, only : message @@ -474,12 +477,12 @@ FUNCTION objective_sm_pd(parameterset) use mo_global_variables, only: nBasins, & ! number of basins L1_sm, L1_sm_mask ! packed measured sm, sm-mask (dim1=ncells, dim2=time) use mo_mhm_constants, only: nodata_dp ! global nodata value - + implicit none real(dp), dimension(:), intent(in) :: parameterset real(dp) :: objective_sm_pd ! objective function value - + ! local integer(i4) :: iBasin ! basin loop counter integer(i4) :: iTime ! time loop counter @@ -504,7 +507,7 @@ FUNCTION objective_sm_pd(parameterset) do iBasin=1, nBasins ! get basin information - call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1, mask=mask1) + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1, mask=mask1) ! allocate allocate(mask_times (size(sm_opti, dim=2))) @@ -521,7 +524,7 @@ FUNCTION objective_sm_pd(parameterset) do iTime = 1, size(sm_opti, dim=2) mat1 = unpack( L1_sm(s1:e1,iTime), mask1, nodata_dp) mat2 = unpack( sm_opti(s1:e1,iTime), mask1, nodata_dp) - mask_sm = unpack(L1_sm_mask(s1:e1,iTime), mask1, .FALSE.) + mask_sm = unpack(L1_sm_mask(s1:e1,iTime), mask1, .FALSE.) pd_time_series = PD(mat1, mat2, mask=mask_sm, valid=mask_times(itime)) end do @@ -537,9 +540,9 @@ FUNCTION objective_sm_pd(parameterset) end do objective_sm_pd = objective_sm_pd**onesixth - + call message(' objective_sm_pd = ', num2str(objective_sm_pd,'(F9.5)')) - + END FUNCTION objective_sm_pd ! ------------------------------------------------------------------ @@ -549,11 +552,11 @@ END FUNCTION objective_sm_pd !> \brief Objective function for soil moisture. - !> \details The objective function only depends on a parameter vector. - !> The model will be called with that parameter vector and + !> \details The objective function only depends on a parameter vector. + !> The model will be called with that parameter vector and !> the model output is subsequently compared to observed data.\n !> - !> Therefore the sum of squared errors (SSE) of the standard score of observed and + !> Therefore the sum of squared errors (SSE) of the standard score of observed and !> modeled soil moisture is calculated. The standard score or normalization (anomaly) !> make the objctive function bias insensitive and basically the dynamics of the soil moisture !> is tried to capture by this obejective function. @@ -564,9 +567,9 @@ END FUNCTION objective_sm_pd !> \f$ SM_{sim} \f$ = simulated soil moisture.\n !> \f$ K \f$ = valid elements in study domain.\n ! - !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 - !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. - !> \f[ OF = \sqrt[6]{\sum(\phi_{i}/N)^6 }. \f] \n + !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 + !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. + !> \f[ OF = \sqrt[6]{\sum(\phi_{i}/N)^6 }. \f] \n !> The observed data L1_sm, L1_sm_mask are global in this module. ! INTENT(IN) @@ -588,7 +591,7 @@ END FUNCTION objective_sm_pd ! None ! RETURN - !> \return real(dp) :: objective_sm_sse_standard_score — objective function value + !> \return real(dp) :: objective_sm_sse_standard_score — objective function value !> (which will be e.g. minimized by an optimization routine like DDS) ! RESTRICTIONS @@ -606,7 +609,7 @@ END FUNCTION objective_sm_pd !> \date March 2015 FUNCTION objective_sm_sse_standard_score(parameterset) - + use mo_init_states, only : get_basin_info use mo_mhm_eval, only : mhm_eval use mo_message, only : message @@ -632,7 +635,7 @@ FUNCTION objective_sm_sse_standard_score(parameterset) real(dp), parameter :: onesixth = 1.0_dp/6.0_dp real(dp), dimension(:,:), allocatable :: sm_opti ! simulated soil moisture ! ! (dim1=ncells, dim2=time) - + call mhm_eval(parameterset, sm_opti=sm_opti) ! initialize some variables @@ -641,11 +644,11 @@ FUNCTION objective_sm_sse_standard_score(parameterset) ! loop over basin - for applying power law later on do iBasin=1, nBasins - ! init + ! init objective_sm_sse_standard_score_basin = 0.0_dp ! get basin information - call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) - + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) + ! standard_score signal is calculated on individual grid cells do iCell = s1, e1 @@ -658,7 +661,7 @@ FUNCTION objective_sm_sse_standard_score(parameterset) objective_sm_sse_standard_score_basin = objective_sm_sse_standard_score_basin + & SSE( standard_score(L1_sm(iCell,:), mask=L1_sm_mask(iCell,:)), & standard_score(sm_opti(iCell,:), mask=L1_sm_mask(iCell,:)), mask=L1_sm_mask(iCell,:)) - + end do print*, iBasin, objective_sm_sse_standard_score_basin ! calculate average soil moisture correlation over all basins with power law @@ -668,9 +671,9 @@ FUNCTION objective_sm_sse_standard_score(parameterset) end do objective_sm_sse_standard_score = objective_sm_sse_standard_score**onesixth - + call message(' objective_sm_sse_standard_score = ', num2str(objective_sm_sse_standard_score,'(E12.5)')) - + END FUNCTION objective_sm_sse_standard_score @@ -682,7 +685,7 @@ END FUNCTION objective_sm_sse_standard_score !> \brief Objective function of KGE for runoff and RMSE for basin_avg TWS (standarized scores) !> \details Objective function of KGE for runoff and RMSE for basin_avg TWS (standarized scores) - + ! INTENT(IN) !> \param[in] "real(dp) :: parameterset(:)" 1D-array with parameters the model is run with @@ -702,20 +705,20 @@ END FUNCTION objective_sm_sse_standard_score ! None ! RETURN - !> \return real(dp) :: objective_kge_q_rmse_tws — objective function value + !> \return real(dp) :: objective_kge_q_rmse_tws — objective function value !> (which will be e.g. minimized by an optimization routine like DDS) ! RESTRICTIONS !> \note Input values must be floating points. \n !> Actually, \f$ KGE \f$ will be returned such that it can be minimized. \n !> For TWS required at least 5 years of data! - + ! EXAMPLE ! para = (/ 1., 2, 3., -999., 5., 6. /) ! obj_value = objective_kge(para) ! LITERATURE - !> Gupta, Hoshin V., et al. "Decomposition of the mean squared error and NSE performance criteria: + !> Gupta, Hoshin V., et al. "Decomposition of the mean squared error and NSE performance criteria: !> Implications for improving hydrological modelling." Journal of Hydrology 377.1 (2009): 80-91. @@ -727,7 +730,7 @@ END FUNCTION objective_sm_sse_standard_score FUNCTION objective_kge_q_rmse_tws(parameterset) use mo_mhm_eval, only: mhm_eval - use mo_constants, only: eps_dp + use mo_constants, only: eps_dp use mo_errormeasures, only: kge, rmse use mo_global_variables, only: nBasins, evalPer use mo_julian, only: caldat @@ -739,8 +742,8 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) use mo_string_utils, only: num2str #ifdef mrm2mhm use mo_mrm_objective_function_runoff, only: extract_runoff -#endif - +#endif + implicit none real(dp), dimension(:), intent(in) :: parameterset @@ -763,21 +766,21 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) real(dp), dimension(:), allocatable :: tws_sim ! simulated tws real(dp), dimension(:), allocatable :: tws_obs ! measured tws - logical, dimension(:), allocatable :: tws_obs_mask ! mask for measured tws + logical, dimension(:), allocatable :: tws_obs_mask ! mask for measured tws integer(i4) :: nMonths ! total number of months integer(i4), dimension(:),allocatable :: month_classes ! vector with months' classes - + ! real(dp), DIMENSION(:), allocatable :: tws_sim_m, tws_obs_m ! monthly values original time series real(dp), DIMENSION(:), allocatable :: tws_sim_m_anom, tws_obs_m_anom ! monthly values anomaly time series - logical, DIMENSION(:), allocatable :: tws_obs_m_mask ! + logical, DIMENSION(:), allocatable :: tws_obs_m_mask ! real(dp), dimension(:), allocatable :: rmse_tws, kge_q ! rmse_tws(nBasins), kge_q(nGaugesTotal) real(dp) :: rmse_tws_avg, kge_q_avg ! obj. functions ! obtain hourly values of runoff and tws: call mHM_eval(parameterset, runoff=runoff, basin_avg_tws=tws) - + !-------------------------------------------- !! TWS !-------------------------------------------- @@ -785,20 +788,20 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) ! allocate allocate( rmse_tws(nBasins) ) rmse_tws(:) = nodata_dp - + do ii=1, nBasins - + ! extract tws the same way as runoff using mrm call extract_basin_avg_tws( ii, tws, tws_sim, tws_obs, tws_obs_mask ) ! check for whether to proceed with this basin or not ! potentially 5 years of data - if (count(tws_obs_mask) .lt. 365*5 ) then + if (count(tws_obs_mask) .lt. 365*5 ) then deallocate (tws_sim, tws_obs, tws_obs_mask) call message('objective_kge_q_rmse_tws: Length of TWS data of Basin ', trim(adjustl(num2str(ii))) , & ' less than 5 years: these data are not considered') - cycle - else + cycle + else ! get initial time of the evaluation period initTime(ii) = real(evalPer(ii)%julStart,dp) @@ -820,7 +823,7 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) allocate ( month_classes( nMonths ) ) allocate ( tws_obs_m_mask( nMonths ) ) allocate ( tws_obs_m_anom( nMonths ) ) - allocate ( tws_sim_m_anom( nMonths ) ) + allocate ( tws_sim_m_anom( nMonths ) ) month_classes(:) = 0 tws_obs_m_mask(:) = .TRUE. @@ -840,10 +843,10 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) ! define mask for missing data in observations (there are always data for simulations) where( abs( tws_obs_m - nodata_dp) .lt. eps_dp ) tws_obs_m_mask = .FALSE. - + ! calculate standard score - tws_obs_m_anom = classified_standard_score(tws_obs_m, month_classes, mask = tws_obs_m_mask) - tws_sim_m_anom = classified_standard_score(tws_sim_m, month_classes, mask = tws_obs_m_mask) + tws_obs_m_anom = classified_standard_score(tws_obs_m, month_classes, mask = tws_obs_m_mask) + tws_sim_m_anom = classified_standard_score(tws_sim_m, month_classes, mask = tws_obs_m_mask) rmse_tws(ii) = rmse(tws_sim_m_anom,tws_obs_m_anom, mask = tws_obs_m_mask) deallocate ( month_classes ) @@ -851,12 +854,12 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) deallocate ( tws_sim_m ) deallocate ( tws_obs_m_mask ) deallocate ( tws_sim_m_anom ) - deallocate ( tws_obs_m_anom ) + deallocate ( tws_obs_m_anom ) deallocate (tws_sim, tws_obs, tws_obs_mask) endif - + end do - + rmse_tws_avg = sum( rmse_tws(:), abs( rmse_tws - nodata_dp) .gt. eps_dp ) / & real(count( abs( rmse_tws - nodata_dp) .gt. eps_dp),dp) deallocate(rmse_tws) @@ -864,7 +867,7 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) !-------------------------------------------- !! RUNOFF !-------------------------------------------- -#ifdef mrm2mhm +#ifdef mrm2mhm nGaugesTotal = size(runoff, dim=2) allocate( kge_q(nGaugesTotal)) kge_q(:) = nodata_dp @@ -877,15 +880,15 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) ! check for whether to proceed with this basin or not ! potentially 5 years of data pp = count(runoff_agg .ge. 0.0_dp ) - if (pp .lt. 365*5 ) then + if (pp .lt. 365*5 ) then deallocate (runoff_agg, runoff_obs, runoff_obs_mask) - cycle - else + cycle + else ! calculate KGE for each basin: kge_q(gg) = kge( runoff_obs, runoff_agg, mask=runoff_obs_mask) - deallocate (runoff_agg, runoff_obs, runoff_obs_mask) + deallocate (runoff_agg, runoff_obs, runoff_obs_mask) end if - + end do ! calculate average KGE value for runoff @@ -895,11 +898,11 @@ FUNCTION objective_kge_q_rmse_tws(parameterset) #else call message('***ERROR: objective_kge_q_rmse_tws: missing routing module for optimization') stop -#endif +#endif - ! + ! objective_kge_q_rmse_tws = rmse_tws_avg * ( 1._dp - kge_q_avg ) - + call message(' objective_kge_q_rmse_tws = ', num2str(objective_kge_q_rmse_tws,'(F9.5)')) END FUNCTION objective_kge_q_rmse_tws @@ -966,7 +969,7 @@ subroutine extract_basin_avg_tws( basinId, tws, tws_sim, tws_obs, tws_obs_mask ) use mo_global_variables, only: basin_avg_TWS_obs, nMeasPerDay_TWS, evalPer, warmingDays, nTstepDay use mo_message, only: message - use mo_constants, only: eps_dp + use mo_constants, only: eps_dp use mo_mrm_constants, only: nodata_dp implicit none @@ -979,7 +982,7 @@ subroutine extract_basin_avg_tws( basinId, tws, tws_sim, tws_obs, tws_obs_mask ) real(dp), dimension(:), allocatable, intent(out) :: tws_sim ! aggregated simulated ! tws to the resolution ! of the measurement - real(dp), dimension(:), allocatable, intent(out) :: tws_obs ! extracted measured + real(dp), dimension(:), allocatable, intent(out) :: tws_obs ! extracted measured ! tws to exactly the ! evaluation period logical, dimension(:), allocatable, intent(out) :: tws_obs_mask ! mask of no data values @@ -1039,4 +1042,147 @@ subroutine extract_basin_avg_tws( basinId, tws, tws_sim, tws_obs, tws_obs_mask ) end subroutine extract_basin_avg_tws + ! ------------------------------------------------------------------ + + ! NAME + ! objective_neutrons_kge_catchment_avg + + !> \brief Objective function for neutrons. + + !> \details The objective function only depends on a parameter vector. + !> The model will be called with that parameter vector and + !> the model output is subsequently compared to observed data.\n + !> + !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average + !> neutrons (N) is calculated + !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f] + !> where \n + !> \f$ r \f$ = Pearson product-moment correlation coefficient \n + !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM \n + !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation \n + !> is calculated and the objective function for a given basin \f$ i \f$ is + !> \f[ \phi_{i} = 1.0 - KGE_{i} \f] + !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods. + !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.\n + !> + !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6 + !> norm to combine the \f$ \phi_{i} \f$ from all basins \f$ N \f$. + !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f] \n + !> The observed data L1_neutronsdata, L1_neutronsdata_mask are global in this module. + + ! INTENT(IN) + !> \param[in] "real(dp) :: parameterset(:)" 1D-array with parameters the model is run with + + ! INTENT(INOUT) + ! None + + ! INTENT(OUT) + ! None + + ! INTENT(IN), OPTIONAL + ! None + + ! INTENT(INOUT), OPTIONAL + ! None + + ! INTENT(OUT), OPTIONAL + ! None + + ! RETURN + !> \return real(dp) :: objective_neutrons_kge_catchment_avg — objective function value + !> (which will be e.g. minimized by an optimization routine) + + ! RESTRICTIONS + !> \note Input values must be floating points. \n + + ! EXAMPLE + ! para = (/ 1., 2, 3., -999., 5., 6. /) + ! obj_value = objective_neutrons_kge_catchment_avg(para) + + ! LITERATURE + ! none + + ! HISTORY + !> \author Martin Schroen + !> \date Jun 2015 + + FUNCTION objective_neutrons_kge_catchment_avg(parameterset) + + use mo_init_states, only : get_basin_info + use mo_mhm_eval, only : mhm_eval + use mo_message, only : message + use mo_moment, only : average + use mo_errormeasures, only : KGE + use mo_string_utils, only : num2str + ! + use mo_global_variables, only: nBasins, & ! number of basins + L1_neutronsdata, L1_neutronsdata_mask ! packed measured neutrons, neutrons-mask (dim1=ncells, dim2=time) + use mo_mhm_constants, only: nodata_dp ! global nodata value + + implicit none + + real(dp), dimension(:), intent(in) :: parameterset + real(dp) :: objective_neutrons_kge_catchment_avg + + ! local + integer(i4) :: iBasin ! basin loop counter + integer(i4) :: iTime ! time loop counter + integer(i4) :: nrows1, ncols1 ! level 1 number of culomns and rows + integer(i4) :: s1, e1 ! start and end index for the current basin + integer(i4) :: ncells1 ! ncells1 of level 1 + real(dp), parameter :: onesixth = 1.0_dp/6.0_dp ! for sixth root + real(dp), dimension(:), allocatable :: neutrons_catch_avg_basin ! spatial average of observed neutrons + real(dp), dimension(:), allocatable :: neutrons_opti_catch_avg_basin ! spatial avergae of modeled neutrons + real(dp), dimension(:,:), allocatable :: neutrons_opti ! simulated neutrons + ! ! (dim1=ncells, dim2=time) + logical, dimension(:), allocatable :: mask_times ! mask for valid neutrons catchment avg time steps + + call mhm_eval(parameterset, neutrons_opti=neutrons_opti) + + ! initialize some variables + objective_neutrons_kge_catchment_avg = 0.0_dp + + ! loop over basin - for applying power law later on + do iBasin=1, nBasins + + ! get basin information + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, iStart=s1, iEnd=e1 ) + + ! allocate + allocate(mask_times (size(neutrons_opti, dim=2))) + allocate(neutrons_catch_avg_basin (size(neutrons_opti, dim=2))) + allocate(neutrons_opti_catch_avg_basin(size(neutrons_opti, dim=2))) + + ! initalize + mask_times = .TRUE. + neutrons_catch_avg_basin = nodata_dp + neutrons_opti_catch_avg_basin = nodata_dp + + ! calculate catchment average soil moisture + do iTime = 1, size(neutrons_opti, dim=2) + + ! check for enough data points in time for correlation + if ( all(.NOT. L1_neutronsdata_mask(:,iTime))) then + write (*,*) 'WARNING: neutrons data at time ', iTime, ' is empty.' + !call message('WARNING: objective_neutrons_kge_catchment_avg: ignored currrent time step since less than') + !call message(' 10 valid cells available in soil moisture observation') + mask_times(iTime) = .FALSE. + cycle + end if + neutrons_catch_avg_basin(iTime) = average( L1_neutronsdata(s1:e1,iTime), mask=L1_neutronsdata_mask(s1:e1,iTime)) + neutrons_opti_catch_avg_basin(iTime) = average(neutrons_opti(s1:e1,iTime), mask=L1_neutronsdata_mask(s1:e1,iTime)) + end do + + ! calculate average neutrons KGE over all basins with power law + ! basins are weighted equally ( 1 / real(nBasin,dp))**6 + objective_neutrons_kge_catchment_avg = objective_neutrons_kge_catchment_avg + & + ( (1.0_dp-KGE(neutrons_catch_avg_basin, neutrons_opti_catch_avg_basin, mask=mask_times)) / real(nBasins,dp) )**6 + end do + + objective_neutrons_kge_catchment_avg = objective_neutrons_kge_catchment_avg**onesixth + + call message(' objective_neutrons_kge_catchment_avg = ', num2str(objective_neutrons_kge_catchment_avg,'(F9.5)')) + +END FUNCTION objective_neutrons_kge_catchment_avg + END MODULE mo_objective_function diff --git a/src/mHM/mo_read_config.f90 b/src/mHM/mo_read_config.f90 index 840ecc26..4ac76215 100644 --- a/src/mHM/mo_read_config.f90 +++ b/src/mHM/mo_read_config.f90 @@ -149,6 +149,7 @@ subroutine read_config() nSoilHorizons_sm_input, & ! No. of mhm soil horizons equivalent to soil moisture input basin_avg_TWS_obs, & ! basin avg TWS data fileTWS, & ! directory with basin average tws data + dirNeutrons, timeStep_neutrons_input, & ! directory where neutron data is located HorizonDepth_mHM, nSoilHorizons_mHM, tillageDepth, & ! soil horizons info for mHM fracSealed_cityArea, nLcoverScene, & ! land cover information LCfilename, LCyearId, & ! @@ -186,8 +187,8 @@ subroutine read_config() sce_ngs, sce_npg, sce_nps, & ! SCE: # complexes, # points per complex, ! ! # points per subcomplex mcmc_opti, & ! MCMC: if optimization mode of MCMC or only uncertainty estimation - mcmc_error_params ! parameters of error model used in likelihood - + mcmc_error_params ! parameters of error model used in likelihood + implicit none ! LOCAL variables @@ -286,11 +287,12 @@ subroutine read_config() ! ! used when timeStep_LAI_input<0 character(256), dimension(maxNoBasins) :: dir_soil_moisture ! soil moisture input character(256), dimension(maxNoBasins) :: file_TWS ! total water storage input file + character(256), dimension(maxNoBasins) :: dir_neutrons ! ground albedo neutron input ! integer(i4) :: nLCover_scene ! given number of land cover scenes - integer(i4), dimension(maxNLCovers) :: LCoverYearStart ! starting year LCover - integer(i4), dimension(maxNLCovers) :: LCoverYearEnd ! ending year LCover + integer(i4), dimension(maxNLCovers) :: LCoverYearStart ! starting year LCover + integer(i4), dimension(maxNLCovers) :: LCoverYearEnd ! ending year LCover character(256), dimension(maxNLCovers) :: LCoverfName ! filename of Lcover file real(dp), dimension(maxGeoUnit, nColPars) :: GeoParam ! geological parameters ! @@ -314,7 +316,12 @@ subroutine read_config() dir_MaxTemperature, dir_absVapPressure, dir_windspeed, & dir_NetRadiation, dir_gridded_LAI ! optional data used for optimization - namelist /optional_data/ dir_soil_moisture, nSoilHorizons_sm_input, timeStep_sm_input, file_TWS + namelist /optional_data/ & + dir_soil_moisture, & + nSoilHorizons_sm_input, & + timeStep_sm_input, & + file_TWS, & + dir_neutrons ! namelist spatial & temporal resolution, otmization information namelist /mainconfig/ timestep, iFlag_cordinate_sys, resolution_Hydrology, resolution_Routing, & L0Basin, optimize, optimize_restart, opti_method, opti_function, nBasins, read_restart, & @@ -544,19 +551,21 @@ subroutine read_config() call message(' defined number of soil horizions: ', trim(num2str(maxNoSoilHorizons)),'!') stop end if + dirNeutrons = dir_neutrons(1:nBasins) + timeStep_neutrons_input = -1 ! daily, hard-coded, to be flexibilized !=============================================================== ! Read evaluation basin average TWS data !=============================================================== fileTWS = file_TWS (1:nBasins) - + ! if (opti_function .EQ. 15) then !OROROR allocate(basin_avg_TWS_obs%basinId(nBasins)); basin_avg_TWS_obs%basinId = nodata_i4 allocate(basin_avg_TWS_obs%fName (nBasins)); basin_avg_TWS_obs%fName(:)= num2str(nodata_i4) - + do iBasin = 1, nBasins - if (trim(fileTWS(iBasin)) .EQ. trim(num2str(nodata_i4))) then + if (trim(fileTWS(iBasin)) .EQ. trim(num2str(nodata_i4))) then call message() call message('***ERROR: mhm.nml: Filename of evaluation TWS data ', & ' for subbasin ', trim(adjustl(num2str(iBasin))), & @@ -564,13 +573,13 @@ subroutine read_config() call message(' Error occured in namelist: evaluation_tws') stop end if - + basin_avg_TWS_obs%basinId(iBasin) = iBasin - basin_avg_TWS_obs%fname(iBasin) = trim(file_TWS(iBasin)) + basin_avg_TWS_obs%fname(iBasin) = trim(file_TWS(iBasin)) end do ! end if !OROROR - - + + !=============================================================== ! Read process selection list !=============================================================== @@ -1088,13 +1097,13 @@ subroutine read_config() #ifndef mrm2mhm call message('***ERROR processCase(8) equals 1, but mrm2mhm preprocessor flag is not given in Makefile') stop -#endif +#endif processMatrix(8, 1) = processCase(8) processMatrix(8, 2) = 5_i4 processMatrix(8, 3) = sum(processMatrix(1:8, 2)) call append(global_parameters, dummy_2d_dp) call append(global_parameters_name, (/'dummy', 'dummy', 'dummy', 'dummy', 'dummy'/)) - + case DEFAULT call message() call message('***ERROR: Process description for process "routing" does not exist!') diff --git a/src/mHM/mo_read_optional_data.f90 b/src/mHM/mo_read_optional_data.f90 index 8ce345a0..993f7391 100644 --- a/src/mHM/mo_read_optional_data.f90 +++ b/src/mHM/mo_read_optional_data.f90 @@ -15,8 +15,8 @@ MODULE mo_read_optional_data PRIVATE - PUBLIC :: read_soil_moisture, read_basin_avg_TWS - + PUBLIC :: read_soil_moisture, read_basin_avg_TWS, read_neutrons + ! ------------------------------------------------------------------ CONTAINS @@ -25,7 +25,7 @@ MODULE mo_read_optional_data ! NAME ! read_soil_moisture - + ! PURPOSE !> \brief Read soil moisture data from NetCDF file for calibration @@ -81,7 +81,7 @@ subroutine read_soil_moisture(iBasin) use mo_mhm_constants, only: nodata_dp use mo_global_variables, only: & dirSoil_moisture, & ! directory of meteo input - evalPer, & ! evaluation period + evalPer, & ! evaluation period L1_sm, L1_sm_mask, & ! soil mositure data and mask timeStep_sm_input, nTimeSteps_L1_sm ! input time step (d,m,y), number of time steps implicit none @@ -97,10 +97,10 @@ subroutine read_soil_moisture(iBasin) real(dp), dimension(:,:), allocatable :: L1_data_packed ! packed data at level-1 from 3D to 2D logical, dimension(:,:,:), allocatable :: L1_mask ! mask at level-1 logical, dimension(:,:), allocatable :: L1_mask_packed ! packed mask at level-1 from 3D to 2D - + ! get basic basin information at level-1 - call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, mask=mask1 ) - + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, mask=mask1 ) + ! basin characteristics and read meteo header call message(' Reading soil moisture for basin: ', trim(adjustl(num2str(iBasin))),' ...') call timer_start(1) @@ -112,24 +112,24 @@ subroutine read_soil_moisture(iBasin) allocate( L1_data_packed(nCells1, nTimeSteps_L1_sm)) allocate( L1_mask_packed(nCells1, nTimeSteps_L1_sm)) do t = 1, nTimeSteps_L1_sm - L1_data_packed(:,t) = pack( L1_data(:,:,t), MASK=mask1(:,:) ) - L1_mask_packed(:,t) = pack( L1_mask(:,:,t), MASK=mask1(:,:) ) + L1_data_packed(:,t) = pack( L1_data(:,:,t), MASK=mask1(:,:) ) + L1_mask_packed(:,t) = pack( L1_mask(:,:,t), MASK=mask1(:,:) ) end do - + ! append call append( L1_sm, L1_data_packed(:,:), fill_value=nodata_dp ) call append( L1_sm_mask, L1_mask_packed(:,:), fill_value=.FALSE. ) ! for multi basin calibration number of time steps may vary for different basins if (iBasin .GT. 1) nTimeSteps_L1_sm = size(L1_sm, 2) - + !free space - deallocate(L1_data, L1_data_packed) + deallocate(L1_data, L1_data_packed) call timer_stop(1) call message(' in ', trim(num2str(timer_get(1),'(F9.3)')), ' seconds.') call timer_clear(1) - + end subroutine read_soil_moisture ! --------------------------------------------------------------------------- @@ -139,7 +139,7 @@ end subroutine read_soil_moisture !> \brief Read basin average TWS timeseries from file, the same way runoff is read - !> \details Read basin average TWS timeseries + !> \details Read basin average TWS timeseries !> Allocate global basin_avg_TWS variable that contains the simulated values after the simulation. ! INTENT(IN) @@ -173,7 +173,7 @@ end subroutine read_soil_moisture ! None ! HISTORY - ! \author Oldrich Rakovec, based on modification of mrm_read_discharge by S. Thober + ! \author Oldrich Rakovec, based on modification of mrm_read_discharge by S. Thober ! \date Oct 2015 ! subroutine read_basin_avg_TWS() @@ -210,7 +210,7 @@ subroutine read_basin_avg_TWS() maxTimeSteps = maxval( simPer(1:nBasins)%julEnd - simPer(1:nBasins)%julStart + 1 ) * nTstepDay allocate( basin_avg_TWS_sim(maxTimeSteps, nBasins) ) basin_avg_TWS_sim = nodata_dp - + ! ************************************************ ! READ basin average TWS TIME SERIES ! ************************************************ @@ -229,7 +229,118 @@ subroutine read_basin_avg_TWS() call paste(basin_avg_TWS_obs%TWS, data_dp_1d, nodata_dp) deallocate (data_dp_1d) end do - + end subroutine read_basin_avg_TWS - + + ! ------------------------------------------------------------------ + + ! NAME + ! read_neutrons + + ! PURPOSE + !> \brief Read neutrons data from NetCDF file for calibration + + !> \details This routine reads oberved neutron fields which are used for model + !> calibration. The neutrons file is expected to be called "neutrons.nc" with + !> a variable "neutrons" inside. The data are read only for the evaluation period + !> they are intended to be used for calibration. Neutrons data are only + !> read if one of the corresponding objective functions is chosen. + + ! CALLING SEQUENCE + + ! INTENT(IN) + !> \param[in] "integer(i4) :: iBasin" Basin Id + + ! INTENT(INOUT) + ! None + + ! INTENT(OUT) + ! None + + ! INTENT(IN), OPTIONAL + ! None + + + ! INTENT(INOUT), OPTIONAL + ! None + + ! INTENT(OUT), OPTIONAL + ! None + + ! RETURN + ! None + + ! RESTRICTIONS + + ! EXAMPLE + + ! LITERATURE + ! None + + ! HISTORY + !> \author Martin Schroen + !> \date Jul 2015 + + subroutine read_neutrons(iBasin) + use mo_message, only: message + use mo_string_utils, only: num2str + use mo_init_states, only: get_basin_info + use mo_read_forcing_nc, only: read_forcing_nc + use mo_timer, only: & + timer_start, timer_stop, timer_get, timer_clear ! Timing of processes + use mo_append, only: append ! append data + use mo_mhm_constants, only: nodata_dp + use mo_global_variables, only: & + dirNeutrons, & ! directory of meteo input + evalPer, & ! evaluation period + L1_neutronsdata, L1_neutronsdata_mask, & ! soil mositure data and mask + timeStep_neutrons_input, nTimeSteps_L1_neutrons ! input time step (d,m,y), number of time steps + implicit none + + integer(i4), intent(in) :: iBasin ! Basin Id + + ! local variables + integer(i4) :: t ! loop vars packing L1_data to L1_data_packed + integer(i4) :: nrows1, ncols1 ! level 1 number of culomns and rows + logical, dimension(:,:), allocatable :: mask1 ! mask of level 1 for packing + integer(i4) :: ncells1 ! ncells1 of level 1 + real(dp), dimension(:,:,:), allocatable :: L1_data ! data at level-1 + real(dp), dimension(:,:), allocatable :: L1_data_packed ! packed data at level-1 from 3D to 2D + logical, dimension(:,:,:), allocatable :: L1_mask ! mask at level-1 + logical, dimension(:,:), allocatable :: L1_mask_packed ! packed mask at level-1 from 3D to 2D + + ! get basic basin information at level-1 + call get_basin_info( iBasin, 1, nrows1, ncols1, nCells=nCells1, mask=mask1 ) + + ! basin characteristics and read meteo header + call message(' Reading neutrons for basin: ', trim(adjustl(num2str(iBasin))),' ...') + call timer_start(1) + call read_forcing_nc( dirNeutrons(iBasin), nRows1, nCols1, evalPer(iBasin), trim('neutrons'), L1_data, mask1, & + nctimestep=timeStep_neutrons_input, nocheck=.TRUE., maskout=L1_mask) + + ! pack variables + nTimeSteps_L1_neutrons = size(L1_data, 3) + allocate( L1_data_packed(nCells1, nTimeSteps_L1_neutrons)) + allocate( L1_mask_packed(nCells1, nTimeSteps_L1_neutrons)) + do t = 1, nTimeSteps_L1_neutrons + L1_data_packed(:,t) = pack( L1_data(:,:,t), MASK=mask1(:,:) ) + L1_mask_packed(:,t) = pack( L1_mask(:,:,t), MASK=mask1(:,:) ) + end do + + ! append + call append( L1_neutronsdata, L1_data_packed(:,:), fill_value=nodata_dp ) + call append( L1_neutronsdata_mask, L1_mask_packed(:,:), fill_value=.FALSE. ) + + ! for multi basin calibration number of time steps may vary for different basins + if (iBasin .GT. 1) nTimeSteps_L1_neutrons = size(L1_neutronsdata, 2) + + !free space + deallocate(L1_data, L1_data_packed) + + call timer_stop(1) + call message(' in ', trim(num2str(timer_get(1),'(F9.3)')), ' seconds.') + call timer_clear(1) + + end subroutine read_neutrons + END MODULE mo_read_optional_data diff --git a/src/mHM/mo_write_ascii.f90 b/src/mHM/mo_write_ascii.f90 index edc02f8f..69647356 100644 --- a/src/mHM/mo_write_ascii.f90 +++ b/src/mHM/mo_write_ascii.f90 @@ -19,7 +19,7 @@ MODULE mo_write_ascii ! Written Christoph Schneider, May 2013 ! Modified, Juliane Mai, May 2013 - module version and documentation - ! Modified, Luis Samaniego, Nov 2013 - improving all formats + ! Modified, Luis Samaniego, Nov 2013 - improving all formats ! Modified, Luis Samaniego, Mar 2014 - added inflow gauge information write out ! Modified, Stephan Thober, Jun 2014 - bug fixed: in writing network properties ! Modified, Rohini Kumar, Jun 2014 - bug fixed: writing of max and min value of discharge @@ -39,7 +39,7 @@ MODULE mo_write_ascii ! PURPOSE !> \brief This modules writes the results of the configuration into an ASCII-file - !> \details + !> \details ! INTENT(IN) ! None @@ -66,16 +66,16 @@ MODULE mo_write_ascii ! None ! EXAMPLE - ! + ! ! LITERATURE - ! + ! ! HISTORY !> \author Christoph Schneider !> \date May 2013 ! Modified, Juliane Mai, May 2013 - module version and documentation - ! Stephan Thober, Jun 2014 - bug fix in L11 config print out + ! Stephan Thober, Jun 2014 - bug fix in L11 config print out ! Stephan Thober, Jun 2014 - updated read_restart ! Rohini, Luis , Jul 2015 - updated version, L1 level prints @@ -110,7 +110,7 @@ Subroutine write_configfile() dirTemperature, & dirReferenceET, & dirOut, & - dirRestartOut, & + dirRestartOut, & warmPer, & evalPer, & SimPer, & @@ -119,7 +119,7 @@ Subroutine write_configfile() use mo_common_variables, only: & global_parameters, & global_parameters_name -#ifdef mrm2mhm +#ifdef mrm2mhm use mo_mrm_global_variables, only: & basin_mrm, & gauge, & @@ -130,17 +130,17 @@ Subroutine write_configfile() L11_toN, & L11_rOrder, & L11_label, & - L11_length, & + L11_length, & L11_slope, & L11_ID, & L1_L11_ID, & L1_areaCell, & nGaugesTotal, & nInflowGaugesTotal, & - resolutionRouting, & + resolutionRouting, & dirGauges -#endif - +#endif + implicit none ! @@ -159,26 +159,26 @@ Subroutine write_configfile() call message(' Error-Code', num2str(err) ) stop end if - write(uconfig, 200) + write(uconfig, 200) write(uconfig, 100) 'mHM-UFZ v-'//trim(version) write(uconfig, 100) 'L. Samaniego & R. Kumar, UFZ' - write(uconfig, 200) + write(uconfig, 200) write(uconfig, 100) write(uconfig, 201) ' M A I N mHM C O N F I G U R A T I O N I N F O R M A T I O N ' write(uconfig, 100) write(uconfig, 103) 'Number of basins ', nBasins -#ifdef mrm2mhm +#ifdef mrm2mhm write(uconfig, 103) 'Total No. of nodes ', L11_nCells write(uconfig, 103) 'Total No. of reaches ', L11_nCells-1 -#endif +#endif write(uconfig, 103) 'No. of cells L0 ', L0_nCells write(uconfig, 103) 'No. of cells L1 ', L1_nCells -#ifdef mrm2mhm +#ifdef mrm2mhm if ( processMatrix(8,1) .ne. 0 ) then write(uconfig, 103) 'No. of cells L11 ', L11_nCells write(uconfig, 103) 'Total No. of gauges ', nGaugesTotal end if -#endif +#endif write(uconfig, 103) 'Time Step [h] ', timeStep do i=1, nBasins select case (iFlag_cordinate_sys) @@ -188,22 +188,22 @@ Subroutine write_configfile() if ( processMatrix(8,1) .ne. 0 ) then write(uconfig, 301) 'Basin ',i, ' Routing Resolution [m] ', resolutionRouting(i) end if -#endif +#endif case(1) write(uconfig, 302) 'Basin ',i, ' Hydrology Resolution [o] ', resolutionHydrology(i) #ifdef mrm2mhm if ( processMatrix(8,1) .ne. 0 ) then write(uconfig, 302) 'Basin ',i, ' Routing Resolution [o] ', resolutionRouting(i) end if -#endif - +#endif + end select end do write(uconfig, 126) 'Flag READ restart ', read_restart write(uconfig, 126) 'Flag WRITE restart ', write_restart ! !****************** - ! Model Run period + ! Model Run period !****************** do j = 1, nBasins write(uconfig, 115) ' Model Run Periods for Basin ', num2str(j) @@ -212,10 +212,10 @@ Subroutine write_configfile() ' Day Month Year Day Month Year' write(uconfig,117) & 'Warming Period (1) ',& - warmPer(j)%dStart, warmPer(j)%mStart, warmPer(j)%yStart ,& + warmPer(j)%dStart, warmPer(j)%mStart, warmPer(j)%yStart ,& warmPer(j)%dEnd , warmPer(j)%mEnd , warmPer(j)%yEnd ,& 'Evaluation Period (2) ',& - evalPer(j)%dStart ,evalPer(j)%mStart , evalPer(j)%yStart ,& + evalPer(j)%dStart ,evalPer(j)%mStart , evalPer(j)%yStart ,& evalPer(j)%dEnd ,evalPer(j)%mEnd , evalPer(j)%yEnd ,& 'Simulation Period (1)+(2) ',& SimPer(j)%dStart , SimPer(j)%mStart , SimPer(j)%yStart ,& @@ -223,7 +223,7 @@ Subroutine write_configfile() end do !********************************* - ! Model Land Cover Observations + ! Model Land Cover Observations !********************************* do j = 1, nBasins write(uconfig,118) ' Land Cover Observations for Basin ', num2str(i) @@ -247,11 +247,11 @@ Subroutine write_configfile() i, global_parameters(i,1), global_parameters(i,2), global_parameters(i,3), & trim(adjustl(global_parameters_name(i))) end do -#ifdef mrm2mhm +#ifdef mrm2mhm ! basin runoff data if ( processMatrix(8,1) .ne. 0 ) then write(uconfig, 202) ' Basin Runoff Data ' - write(uconfig, 107) ' Gauge No.', ' Basin Id', ' Qmax[m3/s]', ' Qmin[m3/s]' + write(uconfig, 107) ' Gauge No.', ' Basin Id', ' Qmax[m3/s]', ' Qmin[m3/s]' do i=1, nGaugesTotal if( any(gauge%Q(:,i) > nodata_dp) ) then write(uconfig,108) i, gauge%basinId(i), maxval(gauge%Q(:,i), gauge%Q(:,i) > nodata_dp), & @@ -264,7 +264,7 @@ Subroutine write_configfile() ! inflow gauge data if ( nInflowGaugesTotal .GT. 0 ) then write(uconfig, 202) ' Basin Inflow Data ' - write(uconfig, 107) ' Gauge No.', ' Basin Id', ' Qmax[m3/s]', ' Qmin[m3/s]' + write(uconfig, 107) ' Gauge No.', ' Basin Id', ' Qmax[m3/s]', ' Qmin[m3/s]' do i=1, nInflowGaugesTotal if( all(InflowGauge%Q(:,i) > nodata_dp) ) then write(uconfig,108) i, InflowGauge%basinId(i), maxval(InflowGauge%Q(:,i), InflowGauge%Q(:,i) > nodata_dp), & @@ -274,7 +274,7 @@ Subroutine write_configfile() end if end do end if -#endif +#endif ! basin config write(uconfig,218) 'Basin-wise Configuration' do n=1,nBasins @@ -288,18 +288,18 @@ Subroutine write_configfile() write(uconfig, 224) 'Directory to morphological input ', dirMorpho(n) write(uconfig, 224) 'Directory to land cover input ', dirLCover(n) -#ifdef mrm2mhm +#ifdef mrm2mhm if ( processMatrix(8,1) .ne. 0 ) then write(uconfig, 224) 'Directory to gauging station input ', dirGauges(n) end if -#endif +#endif write(uconfig, 224) 'Directory to precipitation input ', dirPrecipitation(n) write(uconfig, 224) 'Directory to temperature input ', dirTemperature(n) write(uconfig, 224) 'Directory to reference ET input ', dirReferenceET(n) write(uconfig, 224) 'Directory to write output by default ', dirOut(n) write(uconfig, 224) 'Directory to write output when restarted ', dirRestartOut(n) -#ifdef mrm2mhm +#ifdef mrm2mhm if ( processMatrix(8,1) .ne. 0 ) then write(uconfig, 102) 'River Network (Routing level)' write(uconfig, 100) 'Label 0 = intermediate draining cell ' @@ -352,12 +352,12 @@ Subroutine write_configfile() end do write(uconfig,114) ' Total[km2]', sum(L1_areaCell(basin%L1_iStart(n): basin%L1_iEnd(n))) end if -#endif +#endif ! end do write(uconfig,*) - close(uconfig) + close(uconfig) !! Formats 100 format (a80) @@ -416,7 +416,7 @@ end Subroutine write_configfile ! CALLING SEQUENCE ! INTENT(IN) - !> \param[in] "real(dp) :: best_OF" best objective function value as returned + !> \param[in] "real(dp) :: best_OF" best objective function value as returned !> by the optimization routine !> \param[in] "real(dp), dimension(:) :: best_paramSet" best associated global parameter set @@ -464,7 +464,7 @@ subroutine write_optifile(best_OF, best_paramSet, param_names) implicit none - real(dp), intent(in) :: best_OF + real(dp), intent(in) :: best_OF real(dp), dimension(:), intent(in) :: best_paramSet character(len=*), dimension(:), intent(in) :: param_names @@ -484,14 +484,14 @@ subroutine write_optifile(best_OF, best_paramSet, param_names) stop end if - ! header + ! header write(formHeader, *) '(a40,',n_params,'a40)' ! len(param_names(1))=256 but only 39 characters taken here ! write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii))), ii=1, n_params) write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii)(1:39))), ii=1, n_params) ! output - write(formParams, *) '( es40.14, ', n_params,'(es40.14) )' + write(formParams, *) '( es40.14, ', n_params,'(es40.14) )' write(uopti, formParams) best_OF, (best_paramSet(ii), ii=1, n_params) ! close file @@ -511,8 +511,8 @@ end subroutine write_optifile ! PURPOSE !> \brief Write final, optimized parameter set in a namelist format. - !> \details Write final, optimized parameter set in a namelist format. - !> Only parameters of processes which were switched on are written to the namelist. + !> \details Write final, optimized parameter set in a namelist format. + !> Only parameters of processes which were switched on are written to the namelist. !> All others are discarded. ! CALLING SEQUENCE @@ -555,7 +555,7 @@ end subroutine write_optifile !> \author Juliane Mai !> \date Dec 2013 - ! Modified, + ! Modified, subroutine write_optinamelist(processMatrix, parameters, maskpara, parameters_name) @@ -566,7 +566,7 @@ subroutine write_optinamelist(processMatrix, parameters, maskpara, parameters_na implicit none - integer(i4), dimension(nProcesses, 3), intent(in) :: processMatrix ! information about which process + integer(i4), dimension(nProcesses, 3), intent(in) :: processMatrix ! information about which process ! ! case was used real(dp), dimension(:,:), intent(in) :: parameters ! (min, max, opti) logical, dimension(size(parameters,1)), intent(in) :: maskpara ! .true. if parameter was calibrated @@ -652,17 +652,19 @@ subroutine write_optinamelist(processMatrix, parameters, maskpara, parameters_na write(uopti_nml,*) '&geoparameter' end if case(10) - if (processMatrix(iProc,1) .eq. 1) then + if (processMatrix(iProc,1) .ge. 1) then write(uopti_nml,*) '&neutrons1' end if end select do iPar=iPar_Start, processMatrix(iProc, 3) + if (maskpara(iPar)) then flag=' 1 ' else flag=' 0 ' end if + write(uopti_nml,*) trim(adjustl(parameters_name(iPar))), ' = ', & parameters(iPar,1), ' , ', & parameters(iPar,2), ' , ', & @@ -687,6 +689,3 @@ subroutine write_optinamelist(processMatrix, parameters, maskpara, parameters_na end subroutine write_optinamelist END MODULE mo_write_ascii - - - diff --git a/src/mHM/mo_write_fluxes_states.f90 b/src/mHM/mo_write_fluxes_states.f90 index 42a5cc77..0eb17a3c 100755 --- a/src/mHM/mo_write_fluxes_states.f90 +++ b/src/mHM/mo_write_fluxes_states.f90 @@ -9,7 +9,7 @@ !> \date Apr 2013 ! Modified: ! David Schaefer, Aug 2015 - major rewrite -! +! module mo_write_fluxes_states use mo_kind , only: i4, dp @@ -26,7 +26,7 @@ module mo_write_fluxes_states logical, pointer :: mask(:,:) !> mask to reconstruct data real(dp), allocatable :: data(:) !> store the data between writes integer(i4) :: counter = 0 !> count the number of updateVariable calls - + contains procedure, public :: updateVariable procedure, public :: writeVariableTimestep @@ -56,13 +56,13 @@ module mo_write_fluxes_states procedure newOutputDataset end interface OutputDataset - + private - + public :: OutputDataset contains - + !------------------------------------------------------------------ ! NAME ! newOutputVariable @@ -129,12 +129,12 @@ end function newOutputVariable ! ! PURPOSE !> \brief Update OutputVariable - !> \details Add the array given as actual argument + !> \details Add the array given as actual argument !> to the derived type's component 'data' ! ! CALLING SEQUENCE ! -> with nc of type(OutputVariable): - ! call var%updateVariable(data) + ! call var%updateVariable(data) ! ! INTENT(IN) !> \param[in] "type(NcDataset) :: nc" -> NcDataset which contains the variable @@ -191,10 +191,10 @@ end subroutine updateVariable ! ! CALLING SEQUENCE ! -> with var of type(OutputVariable): - ! call var%updateVariable(data) + ! call var%updateVariable(data) ! ! INTENT(IN) - !> \param[in] "integer(i4) :: timestep" + !> \param[in] "integer(i4) :: timestep" !> -> index along the time dimension of the netcdf variable ! ! INTENT(INOUT) @@ -241,9 +241,9 @@ subroutine writeVariableTimestep(self, timestep) (/1,1,timestep/)) self%data = 0 self%counter = 0 - + end subroutine writeVariableTimestep - + !------------------------------------------------------------------ ! NAME ! newOutputDataset @@ -258,7 +258,7 @@ end subroutine writeVariableTimestep ! nc = OutputDataset(ibasin, mask1) ! ! INTENT(IN) - !> \param[in] "integer(i4) :: ibasin" -> basin id + !> \param[in] "integer(i4) :: ibasin" -> basin id !> \param[in] "logical :: mask1" -> L1 mask to reconstruct the data ! ! INTENT(INOUT) @@ -304,7 +304,7 @@ function newOutputDataset(ibasin, mask1) result(out) integer(i4), intent(in) :: ibasin logical, target :: mask1(:,:) - type(OutputDataset) :: out + type(OutputDataset) :: out ! local integer(i4) :: ii, nn, ncols, nrows, ncells character(3) :: dtype @@ -318,9 +318,9 @@ function newOutputDataset(ibasin, mask1) result(out) unit = fluxesUnit(ibasin) dims1 = (/"easting ", "northing","time "/) nc = createOutputFile(ibasin) - + ii = 0 - + if (outputFlxState(1)) then ii = ii + 1 tmpvars(ii) = OutputVariable( & @@ -375,7 +375,7 @@ function newOutputDataset(ibasin, mask1) result(out) ncells, mask1, .true.) call writeVariableAttributes(tmpvars(ii), "reservoir of sealed areas (sealedSTW)", "mm") end if - + if (outputFlxState(7)) then ii = ii + 1 tmpvars(ii) = OutputVariable( & @@ -395,7 +395,7 @@ function newOutputDataset(ibasin, mask1) result(out) if (outputFlxState(18)) then ii = ii + 1 tmpvars(ii) = OutputVariable( & - nc%setVariable("Neutrons", dtype, dims1), & + nc%setVariable("neutrons", dtype, dims1), & ncells, mask1, .true.) call writeVariableAttributes(tmpvars(ii), "ground albedo neutrons", "cph") end if @@ -506,7 +506,7 @@ function newOutputDataset(ibasin, mask1) result(out) out = OutputDataset(ibasin, nc, tmpvars(1:ii)) end function newOutputDataset - + !------------------------------------------------------------------ ! NAME ! updateDataset @@ -530,33 +530,33 @@ end function newOutputDataset ! L1_runoffSeal, L1_fastRunoff, L1_slowRunoff , & ! L1_baseflow , L1_percol , L1_infilSoil , & ! L1_preEffect ) - ! + ! ! ! INTENT(IN) !> \param[in] "sidx" -> start index of the basin related data in L1_* arguments !> \param[in] "eidx" -> end index of the basin related data in L1_* arguments - !> \param[in] "L1_fSealed" - !> \param[in] "L1_fNotSealed" - !> \param[in] "L1_inter" - !> \param[in] "L1_snowPack" - !> \param[in] "L1_soilMoist" - !> \param[in] "L1_soilMoistSat" - !> \param[in] "L1_sealSTW" - !> \param[in] "L1_unsatSTW" - !> \param[in] "L1_satSTW" - !> \param[in] "L1_neutrons" - !> \param[in] "L1_pet" - !> \param[in] "L1_aETSoil" - !> \param[in] "L1_aETCanopy" - !> \param[in] "L1_aETSealed" - !> \param[in] "L1_total_runoff" - !> \param[in] "L1_runoffSeal" - !> \param[in] "L1_fastRunoff" - !> \param[in] "L1_slowRunoff" - !> \param[in] "L1_baseflow" - !> \param[in] "L1_percol" - !> \param[in] "L1_infilSoil" - !> \param[in] "L1_preEffect" + !> \param[in] "L1_fSealed" + !> \param[in] "L1_fNotSealed" + !> \param[in] "L1_inter" + !> \param[in] "L1_snowPack" + !> \param[in] "L1_soilMoist" + !> \param[in] "L1_soilMoistSat" + !> \param[in] "L1_sealSTW" + !> \param[in] "L1_unsatSTW" + !> \param[in] "L1_satSTW" + !> \param[in] "L1_neutrons" + !> \param[in] "L1_pet" + !> \param[in] "L1_aETSoil" + !> \param[in] "L1_aETCanopy" + !> \param[in] "L1_aETSealed" + !> \param[in] "L1_total_runoff" + !> \param[in] "L1_runoffSeal" + !> \param[in] "L1_fastRunoff" + !> \param[in] "L1_slowRunoff" + !> \param[in] "L1_baseflow" + !> \param[in] "L1_percol" + !> \param[in] "L1_infilSoil" + !> \param[in] "L1_preEffect" ! ! INTENT(INOUT) ! None @@ -606,7 +606,7 @@ subroutine updateDataset(& L1_runoffSeal, L1_fastRunoff, L1_slowRunoff , & L1_baseflow , L1_percol , L1_infilSoil , & L1_preEffect ) - + use mo_global_variables, only : outputFlxState, nSoilHorizons_mHM class(OutputDataset), intent(inout), target :: self @@ -642,7 +642,7 @@ subroutine updateDataset(& ii = 0 vars => self%vars - + if (outputFlxState(1)) then ii = ii + 1 call vars(ii)%updateVariable(L1_inter(sidx:eidx)) @@ -678,7 +678,7 @@ subroutine updateDataset(& ii = ii + 1 call vars(ii)%updateVariable(L1_sealSTW(sidx:eidx)) end if - + if (outputFlxState(7)) then ii = ii + 1 call vars(ii)%updateVariable(L1_unsatSTW(sidx:eidx)) @@ -753,7 +753,7 @@ subroutine updateDataset(& ii = ii + 1 call vars(ii)%updateVariable(L1_preEffect(sidx:eidx)) end if - + end subroutine updateDataset !------------------------------------------------------------------ @@ -767,7 +767,7 @@ end subroutine updateDataset ! ! CALLING SEQUENCE ! -> with nc of type(OutputDataset) - ! call nc%writeTimestep(timestep) + ! call nc%writeTimestep(timestep) ! ! INTENT(IN) !> \param[in] "integer(i4) :: timestep" The model timestep to write @@ -817,7 +817,7 @@ subroutine writeTimestep(self, timestep) do ii = 1, size(self%vars) call self%vars(ii)%writeVariableTimestep(self%counter) end do - + end subroutine writeTimestep !------------------------------------------------------------------ @@ -831,7 +831,7 @@ end subroutine writeTimestep ! ! CALLING SEQUENCE ! -> with nc of type(OutputDataset): - ! call nc%close() + ! call nc%close() ! ! INTENT(IN) ! None @@ -872,7 +872,7 @@ subroutine close(self) use mo_String_utils, only: num2str use mo_message, only: message - use mo_global_variables, only: dirOut + use mo_global_variables, only: dirOut class(OutputDataset) :: self call self%nc%close() @@ -1011,10 +1011,10 @@ end function createOutputFile !> \brief Write output variable attributes ! ! CALLING SEQUENCE - ! call writeVariableAttributes(var, long_name, unit) + ! call writeVariableAttributes(var, long_name, unit) ! ! INTENT(IN) - !> \param[in] "type(OutputVariable) :: var" + !> \param[in] "type(OutputVariable) :: var" !> \param[in] "character(*) :: long_name" -> variable name !> \param[in] "character(*) :: unit" -> physical unit ! @@ -1060,7 +1060,7 @@ subroutine writeVariableAttributes(var, long_name, unit) call var%nc%setAttribute("coordinates","lat lon") end subroutine writeVariableAttributes - + !------------------------------------------------------------------ ! NAME @@ -1153,7 +1153,7 @@ end subroutine mapCoordinates ! ! INTENT(IN) !> \param[in] "integer(i4) :: iBasin" -> basin number - !> \param[in] "type(gridGeoRef) :: level" -> grid reference + !> \param[in] "type(gridGeoRef) :: level" -> grid reference ! ! INTENT(INOUT) ! None @@ -1201,21 +1201,21 @@ subroutine geoCoordinates(ibasin, level, lat, lon) type(gridGeoRef), intent(in) :: level real(dp), intent(out), allocatable :: lat(:,:), lon(:,:) integer(i4) :: ncols, nrows - integer(i4) :: ii, pos - - nrows = level%nrows(ibasin) + integer(i4) :: ii, pos + + nrows = level%nrows(ibasin) ncols = level%ncols(ibasin) - - pos = 1 - if ( ibasin .gt. 1 ) then - do ii = 1, ibasin -1 - pos = pos + level%ncols(ii) * level%nrows(ii) + + pos = 1 + if ( ibasin .gt. 1 ) then + do ii = 1, ibasin -1 + pos = pos + level%ncols(ii) * level%nrows(ii) end do end if - - lat = reshape(L1_rect_latitude(pos:pos+nrows*ncols-1), (/nrows, ncols/)) - lon = reshape(L1_rect_longitude(pos:pos+nrows*ncols-1), (/nrows, ncols/)) - + + lat = reshape(L1_rect_latitude(pos:pos+nrows*ncols-1), (/nrows, ncols/)) + lon = reshape(L1_rect_longitude(pos:pos+nrows*ncols-1), (/nrows, ncols/)) + end subroutine geoCoordinates !------------------------------------------------------------------