From fb46a3f569c448125cad0b0008896f2c1814d5d0 Mon Sep 17 00:00:00 2001 From: Kyle <kyle.c.klenk@gmail.com> Date: Wed, 28 Sep 2022 17:34:00 +0000 Subject: [PATCH] copied mDecisions from SUMMA --- build/source/engine/mDecisions.f90 | 276 ++++++++++++----------------- 1 file changed, 118 insertions(+), 158 deletions(-) diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90 index 4a863a7..72258c0 100755 --- a/build/source/engine/mDecisions.f90 +++ b/build/source/engine/mDecisions.f90 @@ -19,7 +19,6 @@ ! along with this program. If not, see <http://www.gnu.org/licenses/>. module mDecisions_module -USE data_types,only: var_i USE nrtype USE var_lookup, only: maxvarDecisions ! maximum number of decisions implicit none @@ -60,9 +59,6 @@ integer(i4b),parameter,public :: laiScaling = 72 ! exponential fun ! look-up values for the choice of numerical method integer(i4b),parameter,public :: bEuler = 81 ! home-grown backward Euler solution with long time steps integer(i4b),parameter,public :: sundials = 82 ! SUNDIALS/IDA solution -integer(i4b),parameter,public :: iterative = 81 ! iterative -integer(i4b),parameter,public :: nonIterative = 82 ! non-iterative -integer(i4b),parameter,public :: iterSurfEnergyBal = 83 ! iterate only on the surface energy balance ! look-up values for method used to compute derivative integer(i4b),parameter,public :: numerical = 91 ! numerical solution integer(i4b),parameter,public :: analytical = 92 ! analytical solution @@ -153,27 +149,24 @@ integer(i4b),parameter,public :: windUnload = 322 ! Roesch et al 20 ! look-up values for the choice of energy equation integer(i4b),parameter,public :: enthalpyFD = 323 ! enthalpyFD integer(i4b),parameter,public :: closedForm = 324 ! closedForm -! look-up values for the choice of DAE solver -integer(i4b),parameter,public :: sundialIDA = 325 ! IDA solver form Sundials package -integer(i4b),parameter,public :: backwEuler = 326 ! backward Euler method implemented by Martyn ! ----------------------------------------------------------------------------------------------------------- contains - ! ************************************************************************************************ - ! public subroutine mDecisions: save model decisions as named integers - ! ************************************************************************************************ -subroutine mDecisions(numSteps,err,message) +! ************************************************************************************************ +! public subroutine mDecisions: save model decisions as named integers +! ************************************************************************************************ +subroutine mDecisions(err,message) ! model time structures USE multiconst,only:secprday ! number of seconds in a day USE var_lookup,only:iLookTIME ! named variables that identify indices in the time structures + USE globalData,only:refTime,refJulday ! reference time + USE globalData,only:oldTime ! time from the previous time step + USE globalData,only:startTime,finshTime ! start/end time of simulation USE globalData,only:dJulianStart ! julian day of start time of simulation USE globalData,only:dJulianFinsh ! julian day of end time of simulation - USE globalData,only:startTime,finshTime,refTime - USE globalData,only:refJulday ! julian day of start time of simulation - USE globalData,only:refJulday ! julian day of end time of simulation USE globalData,only:data_step ! length of data step (s) - USE globalData,only:numtim ! the number of timesteps in the model + USE globalData,only:numtim ! number of time steps in the simulation ! model decision structures USE globaldata,only:model_decisions ! model decision structure USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure @@ -188,31 +181,25 @@ subroutine mDecisions(numSteps,err,message) USE time_utils_module,only:extractTime ! extract time info from units string USE time_utils_module,only:compjulday ! compute the julian day USE time_utils_module,only:fracDay ! compute fractional day - USE summaActors_FileManager,only: SIM_START_TM, SIM_END_TM ! time info from control file module - + USE summaFileManager,only: SIM_START_TM, SIM_END_TM ! time info from control file module implicit none ! define output - integer(i4b),intent(inout) :: numSteps integer(i4b),intent(out) :: err ! error code character(*),intent(out) :: message ! error message ! define local variables character(len=256) :: cmessage ! error message for downwind routine - real(dp) :: dsec,dsec_tz ! second + real(rkind) :: dsec,dsec_tz ! second ! initialize error control - err=0;message='mDecisions/' - - ! ------------------------------------------------------------------------------------------------- - ! ------------------------------------------------------------------------------------------------- + err=0; message='mDecisions/' ! read information from model decisions file, and populate model decisions structure call readoption(err,cmessage) if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if - ! ------------------------------------------------------------------------------------------------- ! put reference time information into the time structures - call extractTime(forc_meta(iLookFORCE%time)%varunit, & ! date-time string + call extractTime(forc_meta(iLookFORCE%time)%varunit, & ! date-time string refTime%var(iLookTIME%iyyy), & ! year refTime%var(iLookTIME%im), & ! month refTime%var(iLookTIME%id), & ! day @@ -227,14 +214,14 @@ subroutine mDecisions(numSteps,err,message) ! compute the julian date (fraction of day) for the reference time call compjulday(& - refTime%var(iLookTIME%iyyy), & ! year - refTime%var(iLookTIME%im), & ! month - refTime%var(iLookTIME%id), & ! day - refTime%var(iLookTIME%ih), & ! hour - refTime%var(iLookTIME%imin), & ! minute - 0._dp, & ! second - refJulday, & ! julian date for the start of the simulation - err, cmessage) ! error control + refTime%var(iLookTIME%iyyy), & ! year + refTime%var(iLookTIME%im), & ! month + refTime%var(iLookTIME%id), & ! day + refTime%var(iLookTIME%ih), & ! hour + refTime%var(iLookTIME%imin), & ! minute + 0._rkind, & ! second + refJulday, & ! julian date for the start of the simulation + err, cmessage) ! error control if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if ! put simulation start time information into the time structures @@ -253,18 +240,18 @@ subroutine mDecisions(numSteps,err,message) ! compute the julian date (fraction of day) for the start of the simulation call compjulday(& - startTime%var(iLookTIME%iyyy), & ! year - startTime%var(iLookTIME%im), & ! month - startTime%var(iLookTIME%id), & ! day - startTime%var(iLookTIME%ih), & ! hour - startTime%var(iLookTIME%imin), & ! minute - 0._dp, & ! second - dJulianStart, & ! julian date for the start of the simulation - err, cmessage) ! error control + startTime%var(iLookTIME%iyyy), & ! year + startTime%var(iLookTIME%im), & ! month + startTime%var(iLookTIME%id), & ! day + startTime%var(iLookTIME%ih), & ! hour + startTime%var(iLookTIME%imin), & ! minute + 0._rkind, & ! second + dJulianStart, & ! julian date for the start of the simulation + err, cmessage) ! error control if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if ! put simulation end time information into the time structures - call extractTime(trim(SIM_END_TM), & ! date-time string + call extractTime(trim(SIM_END_TM), & ! date-time string finshTime%var(iLookTIME%iyyy), & ! year finshTime%var(iLookTIME%im), & ! month finshTime%var(iLookTIME%id), & ! day @@ -284,27 +271,25 @@ subroutine mDecisions(numSteps,err,message) finshTime%var(iLookTIME%id), & ! day finshTime%var(iLookTIME%ih), & ! hour finshTime%var(iLookTIME%imin), & ! minute - 0._dp, & ! second + 0._rkind, & ! second dJulianFinsh, & ! julian date for the end of the simulation err, cmessage) ! error control if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if ! check start and finish time - ! write(*,'(a,i4,1x,4(i2,1x))') 'startTime: iyyy, im, id, ih, imin = ', startTime%var(1:5) - ! write(*,'(a,i4,1x,4(i2,1x))') 'finshTime: iyyy, im, id, ih, imin = ', finshTime%var(1:5) + write(*,'(a,i4,1x,4(i2,1x))') 'startTime: iyyy, im, id, ih, imin = ', startTime%var(1:5) + write(*,'(a,i4,1x,4(i2,1x))') 'finshTime: iyyy, im, id, ih, imin = ', finshTime%var(1:5) ! check that simulation end time is > start time if(dJulianFinsh < dJulianStart)then; err=20; message=trim(message)//'end time of simulation occurs before start time'; return; end if ! initialize the old time vector (time from the previous time step) + oldTime%var(:) = startTime%var(:) ! compute the number of time steps - numSteps = nint( (dJulianFinsh - dJulianStart)*secprday/data_step ) + 1 - numTim = numSteps - - ! write(*,'(a,1x,i10)') 'number of time steps = ', numSteps + numtim = nint( (dJulianFinsh - dJulianStart)*secprday/data_step ) + 1 + write(*,'(a,1x,i10)') 'number of time steps = ', numtim - ! ------------------------------------------------------------------------------------------------- ! set Noah-MP options DVEG=3 ! option for dynamic vegetation @@ -322,7 +307,7 @@ subroutine mDecisions(numSteps,err,message) case('CLM_Type'); model_decisions(iLookDECISIONS%soilStress)%iDecision = CLM_Type ! thresholded linear function of matric head case('SiB_Type'); model_decisions(iLookDECISIONS%soilStress)%iDecision = SiB_Type ! exponential of the log of matric head case default - err=10; message=trim(message)//"unknown soil moisture function [option="//trim(model_decisions(iLookDECISIONS%soilStress)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown soil moisture function [option="//trim(model_decisions(iLookDECISIONS%soilStress)%cDecision)//"]"; return end select ! identify the choice of function for stomatal resistance @@ -333,129 +318,103 @@ subroutine mDecisions(numSteps,err,message) case('BallBerryFlex' ); model_decisions(iLookDECISIONS%stomResist)%iDecision = BallBerryFlex ! flexible Ball-Berry scheme case('BallBerryTest' ); model_decisions(iLookDECISIONS%stomResist)%iDecision = BallBerryTest ! flexible Ball-Berry scheme (testing) case default - err=10; message=trim(message)//"unknown stomatal resistance function [option="//trim(model_decisions(iLookDECISIONS%stomResist)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown stomatal resistance function [option="//trim(model_decisions(iLookDECISIONS%stomResist)%cDecision)//"]"; return end select ! identify the leaf temperature controls on photosynthesis + stomatal resistance if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbTempFunc)%cDecision)) - case('q10Func' ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = q10Func - case('Arrhenius' ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = Arrhenius - case default - err=10; message=trim(message)//"unknown leaf temperature function [option="//trim(model_decisions(iLookDECISIONS%bbTempFunc)%cDecision)//"]"; return + case('q10Func' ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = q10Func + case('Arrhenius' ); model_decisions(iLookDECISIONS%bbTempFunc)%iDecision = Arrhenius + case default + err=10; message=trim(message)//"unknown leaf temperature function [option="//trim(model_decisions(iLookDECISIONS%bbTempFunc)%cDecision)//"]"; return end select end if ! identify the humidity controls on stomatal resistance if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbHumdFunc)%cDecision)) - case('humidLeafSurface' ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = humidLeafSurface - case('scaledHyperbolic' ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = scaledHyperbolic - case default - err=10; message=trim(message)//"unknown humidity function [option="//trim(model_decisions(iLookDECISIONS%bbHumdFunc)%cDecision)//"]"; return + case('humidLeafSurface' ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = humidLeafSurface + case('scaledHyperbolic' ); model_decisions(iLookDECISIONS%bbHumdFunc)%iDecision = scaledHyperbolic + case default + err=10; message=trim(message)//"unknown humidity function [option="//trim(model_decisions(iLookDECISIONS%bbHumdFunc)%cDecision)//"]"; return end select end if ! identify functions for electron transport function (dependence of photosynthesis on PAR) if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbElecFunc)%cDecision)) - case('linear' ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linear - case('linearJmax' ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linearJmax - case('quadraticJmax' ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = quadraticJmax - case default - err=10; message=trim(message)//"unknown electron transport function [option="//trim(model_decisions(iLookDECISIONS%bbElecFunc)%cDecision)//"]"; return + case('linear' ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linear + case('linearJmax' ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = linearJmax + case('quadraticJmax' ); model_decisions(iLookDECISIONS%bbElecFunc)%iDecision = quadraticJmax + case default + err=10; message=trim(message)//"unknown electron transport function [option="//trim(model_decisions(iLookDECISIONS%bbElecFunc)%cDecision)//"]"; return end select end if ! identify the use of the co2 compensation point in the stomatal conductance calaculations if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbCO2point)%cDecision)) - case('origBWB' ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = origBWB - case('Leuning' ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = Leuning - case default - err=10; message=trim(message)//"unknown option for the co2 compensation point [option="//trim(model_decisions(iLookDECISIONS%bbCO2point)%cDecision)//"]"; return + case('origBWB' ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = origBWB + case('Leuning' ); model_decisions(iLookDECISIONS%bbCO2point)%iDecision = Leuning + case default + err=10; message=trim(message)//"unknown option for the co2 compensation point [option="//trim(model_decisions(iLookDECISIONS%bbCO2point)%cDecision)//"]"; return end select end if ! identify the iterative numerical solution method used in the Ball-Berry stomatal resistance parameterization if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbNumerics)%cDecision)) - case('NoahMPsolution' ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = NoahMPsolution ! the NoahMP solution (and CLM4): fixed point iteration; max 3 iterations - case('newtonRaphson' ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = newtonRaphson ! full Newton-Raphson iterative solution to convergence - case default - err=10; message=trim(message)//"unknown option for the Ball-Berry numerical solution [option="//trim(model_decisions(iLookDECISIONS%bbNumerics)%cDecision)//"]"; return + case('NoahMPsolution' ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = NoahMPsolution ! the NoahMP solution (and CLM4): fixed point iteration; max 3 iterations + case('newtonRaphson' ); model_decisions(iLookDECISIONS%bbNumerics)%iDecision = newtonRaphson ! full Newton-Raphson iterative solution to convergence + case default + err=10; message=trim(message)//"unknown option for the Ball-Berry numerical solution [option="//trim(model_decisions(iLookDECISIONS%bbNumerics)%cDecision)//"]"; return end select end if ! identify the controls on carbon assimilation if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbAssimFnc)%cDecision)) - case('colimitation' ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = colimitation ! enable colimitation, as described by Collatz et al. (1991) and Sellers et al. (1996) - case('minFunc' ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = minFunc ! do not enable colimitation: use minimum of the three controls on carbon assimilation - case default - err=10; message=trim(message)//"unknown option for the controls on carbon assimilation [option="//trim(model_decisions(iLookDECISIONS%bbAssimFnc)%cDecision)//"]"; return + case('colimitation' ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = colimitation ! enable colimitation, as described by Collatz et al. (1991) and Sellers et al. (1996) + case('minFunc' ); model_decisions(iLookDECISIONS%bbAssimFnc)%iDecision = minFunc ! do not enable colimitation: use minimum of the three controls on carbon assimilation + case default + err=10; message=trim(message)//"unknown option for the controls on carbon assimilation [option="//trim(model_decisions(iLookDECISIONS%bbAssimFnc)%cDecision)//"]"; return end select end if ! identify the scaling of photosynthesis from the leaf to the canopy if(model_decisions(iLookDECISIONS%stomResist)%iDecision >= BallBerryFlex)then select case(trim(model_decisions(iLookDECISIONS%bbCanIntg8)%cDecision)) - case('constantScaling' ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = constantScaling ! constant scaling factor - case('laiScaling' ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = laiScaling ! exponential function of LAI (Leuning, Plant Cell Env 1995: "Scaling from..." [eq 9]) - case default - err=10; message=trim(message)//"unknown option for scaling of photosynthesis from the leaf to the canopy [option="//trim(model_decisions(iLookDECISIONS%bbCanIntg8)%cDecision)//"]"; return + case('constantScaling' ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = constantScaling ! constant scaling factor + case('laiScaling' ); model_decisions(iLookDECISIONS%bbCanIntg8)%iDecision = laiScaling ! exponential function of LAI (Leuning, Plant Cell Env 1995: "Scaling from..." [eq 9]) + case default + err=10; message=trim(message)//"unknown option for scaling of photosynthesis from the leaf to the canopy [option="//trim(model_decisions(iLookDECISIONS%bbCanIntg8)%cDecision)//"]"; return end select end if - ! ************************************ - ! ************************************ - ! ************************************ - ! ************************************ - ! SUNDIALS ADDITIONS - ! identify the numerical method select case(trim(model_decisions(iLookDECISIONS%num_method)%cDecision)) - case('bEuler'); model_decisions(iLookDECISIONS%num_method)%iDecision = bEuler - case('sundials'); model_decisions(iLookDECISIONS%num_method)%iDecision = sundials - case('itertive'); model_decisions(iLookDECISIONS%num_method)%iDecision = iterative ! iterative - case('non_iter'); model_decisions(iLookDECISIONS%num_method)%iDecision = nonIterative ! non-iterative - case('itersurf'); model_decisions(iLookDECISIONS%num_method)%iDecision = iterSurfEnergyBal ! iterate only on the surface energy balance + case('bEuler' ); model_decisions(iLookDECISIONS%num_method)%iDecision = bEuler ! home-grown backward Euler solution with long time steps + case('itertive' ); model_decisions(iLookDECISIONS%num_method)%iDecision = bEuler ! home-grown backward Euler solution (included for backwards compatibility) + case('sundials' ); model_decisions(iLookDECISIONS%num_method)%iDecision = sundials ! SUNDIALS/IDA solution case default - err=10; message=trim(message)//"unknown numerical method [option="//trim(model_decisions(iLookDECISIONS%num_method)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown numerical method [option="//trim(model_decisions(iLookDECISIONS%num_method)%cDecision)//"]"; return end select - ! how to compute heat capacity in energy equation + ! how to compute heat capacity in energy equation select case(trim(model_decisions(iLookDECISIONS%howHeatCap)%cDecision)) case('enthalpyFD'); model_decisions(iLookDECISIONS%howHeatCap)%iDecision = enthalpyFD ! enthalpyFD case('closedForm'); model_decisions(iLookDECISIONS%howHeatCap)%iDecision = closedForm ! closedForm case default - ! TODO: after adding howHeatCap decision in corresponding file we should delete the next line - model_decisions(iLookDECISIONS%howHeatCap)%iDecision = closedForm - ! err=10; message=trim(message)//"unknown Cp computation [option="//trim(model_decisions(iLookDECISIONS%howHeatCap)%cDecision)//"]"; return - end select - - ! how to solve the system of differential equations - select case(trim(model_decisions(iLookDECISIONS%diffEqSolv)%cDecision)) - case('sundialIDA'); model_decisions(iLookDECISIONS%diffEqSolv)%iDecision = sundialIDA ! enthalpyFD - case('backwEuler'); model_decisions(iLookDECISIONS%diffEqSolv)%iDecision = backwEuler ! closedForm - case default - ! TODO: after adding diffEqSolv decision in corresponding file we should delete the next line - model_decisions(iLookDECISIONS%diffEqSolv)%iDecision = sundialIDA - ! err=10; message=trim(message)//"unknown DAE solver [option="//trim(model_decisions(iLookDECISIONS%diffEqSolv)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown Cp computation [option="//trim(model_decisions(iLookDECISIONS%howHeatCap)%cDecision)//"]"; return end select - ! ************************************ - ! ************************************ - ! ************************************ - ! ************************************ - ! END SUNDIALS ADDITTION - ! identify the method used to calculate flux derivatives select case(trim(model_decisions(iLookDECISIONS%fDerivMeth)%cDecision)) case('numericl'); model_decisions(iLookDECISIONS%fDerivMeth)%iDecision = numerical ! numerical case('analytic'); model_decisions(iLookDECISIONS%fDerivMeth)%iDecision = analytical ! analytical case default - err=10; message=trim(message)//"unknown method used to calculate flux derivatives [option="//trim(model_decisions(iLookDECISIONS%fDerivMeth)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown method used to calculate flux derivatives [option="//trim(model_decisions(iLookDECISIONS%fDerivMeth)%cDecision)//"]"; return end select ! identify the method used to determine LAI and SAI @@ -463,7 +422,7 @@ subroutine mDecisions(numSteps,err,message) case('monTable'); model_decisions(iLookDECISIONS%LAI_method)%iDecision = monthlyTable ! LAI/SAI taken directly from a monthly table for different vegetation classes case('specified'); model_decisions(iLookDECISIONS%LAI_method)%iDecision = specified ! LAI/SAI computed from green vegetation fraction and winterSAI and summerLAI parameters case default - err=10; message=trim(message)//"unknown method to determine LAI and SAI [option="//trim(model_decisions(iLookDECISIONS%LAI_method)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown method to determine LAI and SAI [option="//trim(model_decisions(iLookDECISIONS%LAI_method)%cDecision)//"]"; return end select ! identify the canopy interception parameterization @@ -472,7 +431,7 @@ subroutine mDecisions(numSteps,err,message) case('sparseCanopy'); model_decisions(iLookDECISIONS%cIntercept)%iDecision = sparseCanopy case('storageFunc'); model_decisions(iLookDECISIONS%cIntercept)%iDecision = storageFunc case default - err=10; message=trim(message)//"unknown canopy interception parameterization [option="//trim(model_decisions(iLookDECISIONS%cIntercept)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown canopy interception parameterization [option="//trim(model_decisions(iLookDECISIONS%cIntercept)%cDecision)//"]"; return end select ! identify the form of Richards' equation @@ -480,7 +439,7 @@ subroutine mDecisions(numSteps,err,message) case('moisture'); model_decisions(iLookDECISIONS%f_Richards)%iDecision = moisture ! moisture-based form case('mixdform'); model_decisions(iLookDECISIONS%f_Richards)%iDecision = mixdform ! mixed form case default - err=10; message=trim(message)//"unknown form of Richards' equation [option="//trim(model_decisions(iLookDECISIONS%f_Richards)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown form of Richards' equation [option="//trim(model_decisions(iLookDECISIONS%f_Richards)%cDecision)//"]"; return end select ! identify the groundwater parameterization @@ -489,7 +448,7 @@ subroutine mDecisions(numSteps,err,message) case('bigBuckt'); model_decisions(iLookDECISIONS%groundwatr)%iDecision = bigBucket ! a big bucket (lumped aquifer model) case('noXplict'); model_decisions(iLookDECISIONS%groundwatr)%iDecision = noExplicit ! no explicit groundwater parameterization case default - err=10; message=trim(message)//"unknown groundwater parameterization [option="//trim(model_decisions(iLookDECISIONS%groundwatr)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown groundwater parameterization [option="//trim(model_decisions(iLookDECISIONS%groundwatr)%cDecision)//"]"; return end select ! identify the hydraulic conductivity profile @@ -497,7 +456,7 @@ subroutine mDecisions(numSteps,err,message) case('constant'); model_decisions(iLookDECISIONS%hc_profile)%iDecision = constant ! constant hydraulic conductivity with depth case('pow_prof'); model_decisions(iLookDECISIONS%hc_profile)%iDecision = powerLaw_profile ! power-law profile case default - err=10; message=trim(message)//"unknown hydraulic conductivity profile [option="//trim(model_decisions(iLookDECISIONS%hc_profile)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown hydraulic conductivity profile [option="//trim(model_decisions(iLookDECISIONS%hc_profile)%cDecision)//"]"; return end select ! identify the upper boundary conditions for thermodynamics @@ -506,7 +465,7 @@ subroutine mDecisions(numSteps,err,message) case('nrg_flux'); model_decisions(iLookDECISIONS%bcUpprTdyn)%iDecision = energyFlux ! energy flux case('zeroFlux'); model_decisions(iLookDECISIONS%bcUpprTdyn)%iDecision = zeroFlux ! zero flux case default - err=10; message=trim(message)//"unknown upper boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcUpprTdyn)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown upper boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcUpprTdyn)%cDecision)//"]"; return end select ! identify the lower boundary conditions for thermodynamics @@ -514,7 +473,7 @@ subroutine mDecisions(numSteps,err,message) case('presTemp'); model_decisions(iLookDECISIONS%bcLowrTdyn)%iDecision = prescribedTemp ! prescribed temperature case('zeroFlux'); model_decisions(iLookDECISIONS%bcLowrTdyn)%iDecision = zeroFlux ! zero flux case default - err=10; message=trim(message)//"unknown lower boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcLowrTdyn)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown lower boundary conditions for thermodynamics [option="//trim(model_decisions(iLookDECISIONS%bcLowrTdyn)%cDecision)//"]"; return end select ! identify the upper boundary conditions for soil hydrology @@ -522,7 +481,7 @@ subroutine mDecisions(numSteps,err,message) case('presHead'); model_decisions(iLookDECISIONS%bcUpprSoiH)%iDecision = prescribedHead ! prescribed head (volumetric liquid water content for mixed form of Richards' eqn) case('liq_flux'); model_decisions(iLookDECISIONS%bcUpprSoiH)%iDecision = liquidFlux ! liquid water flux case default - err=10; message=trim(message)//"unknown upper boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcUpprSoiH)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown upper boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcUpprSoiH)%cDecision)//"]"; return end select ! identify the lower boundary conditions for soil hydrology @@ -532,7 +491,7 @@ subroutine mDecisions(numSteps,err,message) case('drainage'); model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision = freeDrainage ! free drainage case('zeroFlux'); model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision = zeroFlux ! zero flux case default - err=10; message=trim(message)//"unknown lower boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcLowrSoiH)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown lower boundary conditions for soil hydrology [option="//trim(model_decisions(iLookDECISIONS%bcLowrSoiH)%cDecision)//"]"; return end select ! identify the choice of parameterization for vegetation roughness length and displacement height @@ -541,7 +500,7 @@ subroutine mDecisions(numSteps,err,message) case('CM_QJRMS1988' ); model_decisions(iLookDECISIONS%veg_traits)%iDecision = CM_QJRMS1988 ! Choudhury and Monteith (QJRMS 1998) "A four layer model for the heat budget..." case('vegTypeTable' ); model_decisions(iLookDECISIONS%veg_traits)%iDecision = vegTypeTable ! constant parameters dependent on the vegetation type case default - err=10; message=trim(message)//"unknown parameterization for vegetation roughness length and displacement height [option="//trim(model_decisions(iLookDECISIONS%veg_traits)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown parameterization for vegetation roughness length and displacement height [option="//trim(model_decisions(iLookDECISIONS%veg_traits)%cDecision)//"]"; return end select ! identify the choice of parameterization for the rooting profile @@ -550,7 +509,7 @@ subroutine mDecisions(numSteps,err,message) case('powerLaw','notPopulatedYet'); model_decisions(iLookDECISIONS%rootProfil)%iDecision = powerLaw ! simple power-law rooting profile case('doubleExp'); model_decisions(iLookDECISIONS%rootProfil)%iDecision = doubleExp ! the double exponential function of Xeng et al. (JHM 2001) case default - err=10; message=trim(message)//"unknown parameterization for rooting profile [option="//trim(model_decisions(iLookDECISIONS%rootProfil)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown parameterization for rooting profile [option="//trim(model_decisions(iLookDECISIONS%rootProfil)%cDecision)//"]"; return end select ! identify the choice of parameterization for canopy emissivity @@ -558,7 +517,7 @@ subroutine mDecisions(numSteps,err,message) case('simplExp'); model_decisions(iLookDECISIONS%canopyEmis)%iDecision = simplExp ! simple exponential function case('difTrans'); model_decisions(iLookDECISIONS%canopyEmis)%iDecision = difTrans ! parameterized as a function of diffuse transmissivity case default - err=10; message=trim(message)//"unknown parameterization for canopy emissivity [option="//trim(model_decisions(iLookDECISIONS%canopyEmis)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown parameterization for canopy emissivity [option="//trim(model_decisions(iLookDECISIONS%canopyEmis)%cDecision)//"]"; return end select ! choice of parameterization for snow interception @@ -566,7 +525,7 @@ subroutine mDecisions(numSteps,err,message) case('stickySnow'); model_decisions(iLookDECISIONS%snowIncept)%iDecision = stickySnow ! maximum interception capacity an increasing function of temerature case('lightSnow' ); model_decisions(iLookDECISIONS%snowIncept)%iDecision = lightSnow ! maximum interception capacity an inverse function of new snow density case default - err=10; message=trim(message)//"unknown option for snow interception capacity[option="//trim(model_decisions(iLookDECISIONS%snowIncept)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for snow interception capacity[option="//trim(model_decisions(iLookDECISIONS%snowIncept)%cDecision)//"]"; return end select ! identify the choice of wind profile @@ -574,7 +533,7 @@ subroutine mDecisions(numSteps,err,message) case('exponential' ); model_decisions(iLookDECISIONS%windPrfile)%iDecision = exponential ! exponential wind profile extends to the surface case('logBelowCanopy'); model_decisions(iLookDECISIONS%windPrfile)%iDecision = logBelowCanopy ! logarithmic profile below the vegetation canopy case default - err=10; message=trim(message)//"unknown option for choice of wind profile[option="//trim(model_decisions(iLookDECISIONS%windPrfile)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for choice of wind profile[option="//trim(model_decisions(iLookDECISIONS%windPrfile)%cDecision)//"]"; return end select ! identify the choice of atmospheric stability function @@ -583,7 +542,7 @@ subroutine mDecisions(numSteps,err,message) case('louisinv'); model_decisions(iLookDECISIONS%astability)%iDecision = louisInversePower ! Louis (1979) inverse power function case('mahrtexp'); model_decisions(iLookDECISIONS%astability)%iDecision = mahrtExponential ! Mahrt (1987) exponential case default - err=10; message=trim(message)//"unknown stability function [option="//trim(model_decisions(iLookDECISIONS%astability)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown stability function [option="//trim(model_decisions(iLookDECISIONS%astability)%cDecision)//"]"; return end select ! choice of canopy shortwave radiation method @@ -594,7 +553,7 @@ subroutine mDecisions(numSteps,err,message) case('NL_scatter' ); model_decisions(iLookDECISIONS%canopySrad)%iDecision = NL_scatter ! Simplified method Nijssen and Lettenmaier (JGR 1999) case('BeersLaw' ); model_decisions(iLookDECISIONS%canopySrad)%iDecision = BeersLaw ! Beer's Law (as implemented in VIC) case default - err=10; message=trim(message)//"unknown canopy radiation method [option="//trim(model_decisions(iLookDECISIONS%canopySrad)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown canopy radiation method [option="//trim(model_decisions(iLookDECISIONS%canopySrad)%cDecision)//"]"; return end select ! choice of albedo representation @@ -602,7 +561,7 @@ subroutine mDecisions(numSteps,err,message) case('conDecay'); model_decisions(iLookDECISIONS%alb_method)%iDecision = constantDecay ! constant decay (e.g., VIC, CLASS) case('varDecay'); model_decisions(iLookDECISIONS%alb_method)%iDecision = variableDecay ! variable decay (e.g., BATS approach, with destructive metamorphism + soot content) case default - err=10; message=trim(message)//"unknown option for snow albedo [option="//trim(model_decisions(iLookDECISIONS%alb_method)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for snow albedo [option="//trim(model_decisions(iLookDECISIONS%alb_method)%cDecision)//"]"; return end select ! choice of snow compaction routine @@ -610,7 +569,7 @@ subroutine mDecisions(numSteps,err,message) case('consettl'); model_decisions(iLookDECISIONS%compaction)%iDecision = constantSettlement ! constant settlement rate case('anderson'); model_decisions(iLookDECISIONS%compaction)%iDecision = andersonEmpirical ! semi-empirical method of Anderson (1976) case default - err=10; message=trim(message)//"unknown option for snow compaction [option="//trim(model_decisions(iLookDECISIONS%compaction)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for snow compaction [option="//trim(model_decisions(iLookDECISIONS%compaction)%cDecision)//"]"; return end select ! choice of method to combine and sub-divide snow layers @@ -618,7 +577,7 @@ subroutine mDecisions(numSteps,err,message) case('jrdn1991'); model_decisions(iLookDECISIONS%snowLayers)%iDecision = sameRulesAllLayers ! SNTHERM option: same combination/sub-dividion rules applied to all layers case('CLM_2010'); model_decisions(iLookDECISIONS%snowLayers)%iDecision = rulesDependLayerIndex ! CLM option: combination/sub-dividion rules depend on layer index case default - err=10; message=trim(message)//"unknown option for combination/sub-division of snow layers [option="//trim(model_decisions(iLookDECISIONS%snowLayers)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for combination/sub-division of snow layers [option="//trim(model_decisions(iLookDECISIONS%snowLayers)%cDecision)//"]"; return end select ! choice of thermal conductivity representation for snow @@ -628,7 +587,7 @@ subroutine mDecisions(numSteps,err,message) case('jrdn1991'); model_decisions(iLookDECISIONS%thCondSnow)%iDecision = Jordan1991 ! Jordan (1991) case('smnv2000'); model_decisions(iLookDECISIONS%thCondSnow)%iDecision = Smirnova2000 ! Smirnova et al. (2000) case default - err=10; message=trim(message)//"unknown option for thermal conductivity of snow [option="//trim(model_decisions(iLookDECISIONS%thCondSnow)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for thermal conductivity of snow [option="//trim(model_decisions(iLookDECISIONS%thCondSnow)%cDecision)//"]"; return end select ! choice of thermal conductivity representation for soil @@ -637,7 +596,7 @@ subroutine mDecisions(numSteps,err,message) case('mixConstit' ); model_decisions(iLookDECISIONS%thCondSoil)%iDecision = mixConstit ! mixture of constituents case('hanssonVZJ' ); model_decisions(iLookDECISIONS%thCondSoil)%iDecision = hanssonVZJ ! test case for the mizoguchi lab experiment, Hansson et al. VZJ 2004 case default - err=10; message=trim(message)//"unknown option for thermal conductivity of soil [option="//trim(model_decisions(iLookDECISIONS%thCondSoil)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for thermal conductivity of soil [option="//trim(model_decisions(iLookDECISIONS%thCondSoil)%cDecision)//"]"; return end select ! choice of method for the spatial representation of groundwater @@ -645,7 +604,7 @@ subroutine mDecisions(numSteps,err,message) case('localColumn'); model_decisions(iLookDECISIONS%spatial_gw)%iDecision = localColumn ! separate groundwater in each local soil column case('singleBasin'); model_decisions(iLookDECISIONS%spatial_gw)%iDecision = singleBasin ! single groundwater store over the entire basin case default - err=10; message=trim(message)//"unknown option for spatial representation of groundwater [option="//trim(model_decisions(iLookDECISIONS%spatial_gw)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for spatial representation of groundwater [option="//trim(model_decisions(iLookDECISIONS%spatial_gw)%cDecision)//"]"; return end select ! choice of routing method @@ -653,7 +612,7 @@ subroutine mDecisions(numSteps,err,message) case('timeDlay'); model_decisions(iLookDECISIONS%subRouting)%iDecision = timeDelay ! time-delay histogram case('qInstant'); model_decisions(iLookDECISIONS%subRouting)%iDecision = qInstant ! instantaneous routing case default - err=10; message=trim(message)//"unknown option for sub-grid routing [option="//trim(model_decisions(iLookDECISIONS%subRouting)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for sub-grid routing [option="//trim(model_decisions(iLookDECISIONS%subRouting)%cDecision)//"]"; return end select ! choice of new snow density @@ -664,7 +623,7 @@ subroutine mDecisions(numSteps,err,message) case('pahaut_76'); model_decisions(iLookDECISIONS%snowDenNew)%iDecision = pahaut_76 ! Pahaut 1976, wind speed dependent (derived from Col de Porte, French Alps) case('constDens'); model_decisions(iLookDECISIONS%snowDenNew)%iDecision = constDens ! Constant new snow density case default - err=10; message=trim(message)//"unknown option for new snow density [option="//trim(model_decisions(iLookDECISIONS%snowDenNew)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for new snow density [option="//trim(model_decisions(iLookDECISIONS%snowDenNew)%cDecision)//"]"; return end select ! choice of snow unloading from canopy @@ -672,49 +631,52 @@ subroutine mDecisions(numSteps,err,message) case('meltDripUnload','notPopulatedYet'); model_decisions(iLookDECISIONS%snowUnload)%iDecision = meltDripUnload ! Hedstrom and Pomeroy (1998), Storck et al 2002 (snowUnloadingCoeff & ratioDrip2Unloading) case('windUnload'); model_decisions(iLookDECISIONS%snowUnload)%iDecision = windUnload ! Roesch et al 2001, formulate unloading based on wind and temperature case default - err=10; message=trim(message)//"unknown option for snow unloading [option="//trim(model_decisions(iLookDECISIONS%snowUnload)%cDecision)//"]"; return + err=10; message=trim(message)//"unknown option for snow unloading [option="//trim(model_decisions(iLookDECISIONS%snowUnload)%cDecision)//"]"; return end select + ! ----------------------------------------------------------------------------------------------------------------------------------------------- + ! check for consistency among options + ! ----------------------------------------------------------------------------------------------------------------------------------------------- ! check zero flux lower boundary for topmodel baseflow option select case(model_decisions(iLookDECISIONS%groundwatr)%iDecision) case(qbaseTopmodel) - if(model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision /= zeroFlux)then - message=trim(message)//'lower boundary condition for soil hydology must be zeroFlux with qbaseTopmodel option for groundwater' - err=20; return - end if + if(model_decisions(iLookDECISIONS%bcLowrSoiH)%iDecision /= zeroFlux)then + message=trim(message)//'lower boundary condition for soil hydology must be zeroFlux with qbaseTopmodel option for groundwater' + err=20; return + end if end select ! check power-law profile is selected when using topmodel baseflow option select case(model_decisions(iLookDECISIONS%groundwatr)%iDecision) case(qbaseTopmodel) - if(model_decisions(iLookDECISIONS%hc_profile)%iDecision /= powerLaw_profile)then - message=trim(message)//'power-law transmissivity profile must be selected when using topmodel baseflow option' - err=20; return - end if + if(model_decisions(iLookDECISIONS%hc_profile)%iDecision /= powerLaw_profile)then + message=trim(message)//'power-law transmissivity profile must be selected when using topmodel baseflow option' + err=20; return + end if end select ! check bigBucket groundwater option is used when for spatial groundwater is singleBasin if(model_decisions(iLookDECISIONS%spatial_gw)%iDecision == singleBasin)then if(model_decisions(iLookDECISIONS%groundwatr)%iDecision /= bigBucket)then - message=trim(message)//'groundwater parameterization must be bigBucket when using singleBasin for spatial_gw' - err=20; return + message=trim(message)//'groundwater parameterization must be bigBucket when using singleBasin for spatial_gw' + err=20; return end if end if end subroutine mDecisions - ! ************************************************************************************************ - ! private subroutine readoption: read information from model decisions file - ! ************************************************************************************************ +! ************************************************************************************************ +! private subroutine readoption: read information from model decisions file +! ************************************************************************************************ subroutine readoption(err,message) ! used to read information from model decisions file USE ascii_util_module,only:file_open ! open file USE ascii_util_module,only:linewidth ! max character number for one line USE ascii_util_module,only:get_vlines ! get a vector of non-comment lines - USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files - USE summaActors_FileManager,only:M_DECISIONS ! definition of modeling options + USE summaFileManager,only:SETTINGS_PATH ! path for metadata files + USE summaFileManager,only:M_DECISIONS ! definition of modeling options USE get_ixname_module,only:get_ixdecisions ! identify index of named variable USE globalData,only:model_decisions ! model decision structure implicit none @@ -735,7 +697,7 @@ subroutine readoption(err,message) err=0; message='readoption/' ! build filename infile = trim(SETTINGS_PATH)//trim(M_DECISIONS) - ! write(*,'(2(a,1x))') 'decisions file = ', trim(infile) + write(*,'(2(a,1x))') 'decisions file = ', trim(infile) ! open file call file_open(trim(infile),unt,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; end if @@ -753,13 +715,11 @@ subroutine readoption(err,message) if (err/=0) then; err=30; message=trim(message)//"errorReadLine"; return; end if ! get the index of the decision in the data structure iVar = get_ixdecisions(trim(option)) - ! write(*,'(i4,1x,a)') iDecision, trim(option)//': '//trim(decision) + write(*,'(i4,1x,a)') iDecision, trim(option)//': '//trim(decision) if(iVar<=0)then; err=40; message=trim(message)//"cannotFindDecisionIndex[name='"//trim(option)//"']"; return; end if ! populate the model decisions structure model_decisions(iVar)%cOption = trim(option) model_decisions(iVar)%cDecision = trim(decision) end do end subroutine readoption - - end module mDecisions_module -- GitLab