Upadte summa Merge branch 'develop' of https://github.com/CH-Earth/summa

 into develop

build/Makefile

@@ -300,10 +300,10 @@ DRIVER__EX = summa.exe
 
 # Define version number
 VERSIONFILE = $(DRIVER_DIR)/summaversion.inc
-VERSION = $(shell git tag | tail -n 1)
 BULTTIM = $(shell date)
 GITBRCH = $(shell git describe --long --all --always | sed -e's/heads\///')
 GITHASH = $(shell git rev-parse HEAD)
+VERSION = $(shell git show-ref --tags | grep $GITHASH | sed 's/.*tags\///' | grep . || echo "undefined")
 
 #========================================================================
 # PART 3: Checks

build/source/dshare/flxMapping.f90

@@ -164,9 +164,9 @@ subroutine flxMapping(err,message)
  flux2state_orig(iLookFLUX%scalarAquiferBaseflow)           = flux2state(state1=iname_matLayer,  state2=integerMissing)
 
  ! derived variables
- flux2state_orig(iLookFLUX%scalarTotalET)                   = flux2state(state1=iname_nrgCanopy, state2=integerMissing)
+ flux2state_orig(iLookFLUX%scalarTotalET)                   = flux2state(state1=iname_nrgCanopy, state2=iname_nrgLayer)
  flux2state_orig(iLookFLUX%scalarTotalRunoff)               = flux2state(state1=iname_matLayer,  state2=integerMissing)
- flux2state_orig(iLookFLUX%scalarNetRadiation)              = flux2state(state1=iname_nrgCanopy, state2=integerMissing)
+ flux2state_orig(iLookFLUX%scalarNetRadiation)              = flux2state(state1=iname_nrgCanopy, state2=iname_nrgLayer)
 
  ! ** copy across flux metadata
  do iVar=1,nFlux

build/source/dshare/get_ixname.f90

@@ -558,8 +558,9 @@ function get_ixdiag(varName)
   case('scalarVGn_m'                    ); get_ixdiag = iLookDIAG%scalarVGn_m                      ! van Genuchten "m" parameter (-)
   case('scalarKappa'                    ); get_ixdiag = iLookDIAG%scalarKappa                      ! constant in the freezing curve function (m K-1)
   case('scalarVolLatHt_fus'             ); get_ixdiag = iLookDIAG%scalarVolLatHt_fus               ! volumetric latent heat of fusion     (J m-3)
-  ! number of function evaluations
+  ! timing information
   case('numFluxCalls'                   ); get_ixdiag = iLookDIAG%numFluxCalls                     ! number of flux calls (-)
+  case('wallClockTime'                  ); get_ixdiag = iLookDIAG%wallClockTime                    ! wall clock time (s)
   ! get to here if cannot find the variable
   case default
    get_ixdiag = integerMissing

build/source/dshare/globalData.f90

@@ -305,6 +305,7 @@ MODULE globalData
  logical(lgt),save,public                    :: globalPrintFlag=.false.     ! flag to compute the Jacobian
  integer(i4b),save,public                    :: chunksize=1024              ! chunk size for the netcdf read/write
  integer(i4b),save,public                    :: outputPrecision=nf90_double ! variable type
+ integer(i4b),save,public                    :: outputCompressionLevel=4             ! output netcdf file deflate level: 0-9. 0 is no compression.
 
  ! define result from the time calls
  integer(i4b),dimension(8),save,public       :: startInit,endInit       ! date/time for the start and end of the initialization
@@ -331,14 +332,14 @@ MODULE globalData
  ! output file information
  logical(lgt),dimension(maxvarFreq),save,public :: outFreq              ! true if the output frequency is desired
  integer(i4b),dimension(maxvarFreq),save,public :: ncid                 ! netcdf output file id
- 
+
  ! look-up values for the choice of the time zone information (formerly in modelDecisions module)
  integer(i4b),parameter,public               :: ncTime=1                ! time zone information from NetCDF file (timeOffset = longitude/15. - ncTimeOffset)
  integer(i4b),parameter,public               :: utcTime=2               ! all times in UTC (timeOffset = longitude/15. hours)
  integer(i4b),parameter,public               :: localTime=3             ! all times local (timeOffset = 0)
- 
+
  ! define fixed dimensions
  integer(i4b),parameter,public               :: nBand=2          ! number of spectral bands
- integer(i4b),parameter,public               :: nTimeDelay=2000  ! number of hours in the time delay histogram (default: ~1 season = 24*365/4)
+ integer(i4b),parameter,public               :: nTimeDelay=2000  ! number of time steps in the time delay histogram (default: ~1 season = 24*365/4)
 
 END MODULE globalData

build/source/dshare/popMetadat.f90

@@ -420,8 +420,9 @@ subroutine popMetadat(err,message)
  diag_meta(iLookDIAG%scalarVGn_m)                     = var_info('scalarVGn_m'                    , 'van Genuchten "m" parameter'                                      , '-'               , get_ixVarType('midSoil'), iMissVec, iMissVec, .false.)
  diag_meta(iLookDIAG%scalarKappa)                     = var_info('scalarKappa'                    , 'constant in the freezing curve function'                          , 'm K-1'           , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  diag_meta(iLookDIAG%scalarVolLatHt_fus)              = var_info('scalarVolLatHt_fus'             , 'volumetric latent heat of fusion'                                 , 'J m-3'           , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
- ! number of function evaluations
+ ! timing information
  diag_meta(iLookDIAG%numFluxCalls)                    = var_info('numFluxCalls'                   , 'number of flux calls'                                             , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
+ diag_meta(iLookDIAG%wallClockTime)                   = var_info('wallClockTime'                  , 'wall clock time'                                                  , 's'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
 
  ! -----
  ! * local model fluxes...
@@ -708,6 +709,7 @@ subroutine read_output_file(err,message)
  USE globalData, only: flux_meta               ! data structure for local flux variables
  USE globalData, only: deriv_meta              ! data structure for local flux derivatives
  USE globalData, only: outputPrecision         ! data structure for output precision
+ USE globalData, only: outputCompressionLevel  ! data structure for output netcdf deflate level 
 
  ! structures of named variables
  USE var_lookup, only: iLookTYPE               ! named variables for categorical data
@@ -814,6 +816,19 @@ subroutine read_output_file(err,message)
     end if
     cycle
   end if
+
+  ! set output netcdf file compression level if given. default is level 4.
+  if (trim(varName)=='outputCompressionLevel') then
+    statName = trim(lineWords(nWords))
+    read(statName, *) outputCompressionLevel
+    if ((outputCompressionLevel .LT. 0) .or. (outputCompressionLevel .GT. 9)) then
+      err=20
+      cmessage='outputCompressionLevel must be between 0 and 9.'
+      message=trim(message)//trim(cmessage)//trim(varName);
+      return
+    end if
+    cycle
+  end if
 
   ! --- variables with multiple statistics options --------------------------
 

build/source/dshare/var_lookup.f90

@@ -435,8 +435,9 @@ MODULE var_lookup
   integer(i4b)    :: scalarVGn_m                     = integerMissing ! van Genuchten "m" parameter (-)
   integer(i4b)    :: scalarKappa                     = integerMissing ! constant in the freezing curve function (m K-1)
   integer(i4b)    :: scalarVolLatHt_fus              = integerMissing ! volumetric latent heat of fusion     (J m-3)
-  ! number of function evaluations
+  ! timing information
   integer(i4b)    :: numFluxCalls                    = integerMissing ! number of flux calls (-)
+  integer(i4b)    :: wallClockTime                   = integerMissing ! wall clock time (s)
  endtype iLook_diag
 
  ! ***********************************************************************************************************
@@ -811,7 +812,7 @@ MODULE var_lookup
                                                                          51, 52, 53, 54, 55, 56, 57, 58, 59, 60,&
                                                                          61, 62, 63, 64, 65, 66, 67, 68, 69, 70,&
                                                                          71, 72, 73, 74, 75, 76, 77, 78, 79, 80,&
-                                                                         81, 82, 83)
+                                                                         81, 82, 83, 84)
  ! named variables: model fluxes
  type(iLook_flux),    public,parameter :: iLookFLUX     =iLook_flux    (  1,  2,  3,  4,  5,  6,  7,  8,  9, 10,&
                                                                          11, 12, 13, 14, 15, 16, 17, 18, 19, 20,&

build/source/engine/check_icond.f90

@@ -80,17 +80,18 @@ subroutine check_icond(nGRU,                          & ! number of GRUs and HRU
  integer(i4b)                           :: iHRU          ! loop index
 
  ! temporary variables for realism checks
- integer(i4b)                      :: iLayer             ! index of model layer
- integer(i4b)                      :: iSoil              ! index of soil layer
- real(rkind)                          :: fLiq               ! fraction of liquid water on the vegetation canopy (-)
- real(rkind)                          :: vGn_m              ! van Genutchen "m" parameter (-)
- real(rkind)                          :: tWat               ! total water on the vegetation canopy (kg m-2)
- real(rkind)                          :: scalarTheta        ! liquid water equivalent of total water [liquid water + ice] (-)
- real(rkind)                          :: h1,h2              ! used to check depth and height are consistent
- integer(i4b)                      :: nLayers            ! total number of layers
- real(rkind)                          :: kappa              ! constant in the freezing curve function (m K-1)
- integer(i4b)                      :: nSnow              ! number of snow layers
- real(rkind),parameter                :: xTol=1.e-10_rkind     ! small tolerance to address precision issues
+ integer(i4b)                   :: iLayer                ! index of model layer
+ integer(i4b)                   :: iSoil                 ! index of soil layer
+ real(rkind)                    :: fLiq                  ! fraction of liquid water on the vegetation canopy (-)
+ real(rkind)                    :: vGn_m                 ! van Genutchen "m" parameter (-)
+ real(rkind)                    :: tWat                  ! total water on the vegetation canopy (kg m-2)
+ real(rkind)                    :: scalarTheta           ! liquid water equivalent of total water [liquid water + ice] (-)
+ real(rkind)                    :: h1,h2                 ! used to check depth and height are consistent
+ integer(i4b)                   :: nLayers               ! total number of layers
+ real(rkind)                    :: kappa                 ! constant in the freezing curve function (m K-1)
+ integer(i4b)                   :: nSnow                 ! number of snow layers
+ real(rkind),parameter          :: xTol=1.e-10_rkind     ! small tolerance to address precision issues
+ real(rkind),parameter          :: canIceTol=1.e-3_rkind ! small tolerance to allow existence of canopy ice for above-freezing temperatures (kg m-2)
  ! --------------------------------------------------------------------------------------------------------
 
  ! Start procedure here
@@ -148,15 +149,15 @@ subroutine check_icond(nGRU,                          & ! number of GRUs and HRU
    ! compute the constant in the freezing curve function (m K-1)
    kappa  = (iden_ice/iden_water)*(LH_fus/(gravity*Tfreeze))  ! NOTE: J = kg m2 s-2
 
-   ! modify the liquid water and ice in the canopy
-   if(scalarCanopyIce > 0._rkind .and. scalarCanopyTemp > Tfreeze)then
-    message=trim(message)//'canopy ice > 0 when canopy temperature > Tfreeze'
+   ! check canopy ice content for unrealistic situations
+   if(scalarCanopyIce > canIceTol .and. scalarCanopyTemp > Tfreeze)then
+    ! ice content > threshold, terminate run
+    write(message,'(A,E22.16,A,E9.3,A,F7.3,A,F7.3,A)') trim(message)//'canopy ice (=',scalarCanopyIce,') > canIceTol (=',canIceTol,') when canopy temperature (=',scalarCanopyTemp,') > Tfreeze (=',Tfreeze,')'
     err=20; return
+   else if(scalarCanopyIce > 0._rkind .and. scalarCanopyTemp > Tfreeze)then
+    ! if here, ice content < threshold. Could be sublimation on previous timestep or simply wrong input. Print a warning
+	write(*,'(A,E22.16,A,F7.3,A,F7.3,A)') 'Warning: canopy ice content in restart file (=',scalarCanopyIce,') > 0 when canopy temperature (=',scalarCanopyTemp,') > Tfreeze (=',Tfreeze,'). Continuing.',NEW_LINE('a')
    end if
-   fLiq = fracliquid(scalarCanopyTemp,snowfrz_scale)  ! fraction of liquid water (-)
-   tWat = scalarCanopyLiq + scalarCanopyIce           ! total water (kg m-2)
-   scalarCanopyLiq = fLiq*tWat                        ! mass of liquid water on the canopy (kg m-2)
-   scalarCanopyIce = (1._rkind - fLiq)*tWat              ! mass of ice on the canopy (kg m-2)
 
    ! number of layers
    nLayers = gru_struc(iGRU)%hruInfo(iHRU)%nSnow + gru_struc(iGRU)%hruInfo(iHRU)%nSoil

build/source/engine/coupled_em.f90

@@ -64,9 +64,8 @@ module coupled_em_module
 
 ! look-up values for the numerical method
 USE mDecisions_module,only:         &
- iterative,                         &      ! iterative
- nonIterative,                      &      ! non-iterative
- iterSurfEnergyBal                         ! iterate only on the surface energy balance
+ bEuler,                            &      ! home-grown backward Euler solution with long time steps
+ sundials                                  ! SUNDIALS/IDA solution
 
 ! look-up values for the maximum interception capacity
 USE mDecisions_module,only:         &
@@ -132,15 +131,11 @@ subroutine coupled_em(&
  ! the model solver
  USE indexState_module,only:indexState      ! define indices for all model state variables and layers
  USE opSplittin_module,only:opSplittin      ! solve the system of thermodynamic and hydrology equations for a given substep
+ USE time_utils_module,only:elapsedSec      ! calculate the elapsed time
  ! additional subroutines
  USE tempAdjust_module,only:tempAdjust      ! adjust snow temperature associated with new snowfall
  USE snwDensify_module,only:snwDensify      ! snow densification (compaction and cavitation)
  USE var_derive_module,only:calcHeight      ! module to calculate height at layer interfaces and layer mid-point
- ! look-up values for the numerical method
- USE mDecisions_module,only:         &
-  iterative,                         &      ! iterative
-  nonIterative,                      &      ! non-iterative
-  iterSurfEnergyBal                         ! iterate only on the surface energy balance
  ! look-up values for the maximum interception capacity
  USE mDecisions_module,only:          &
                        stickySnow,    &      ! maximum interception capacity an increasing function of temerature
@@ -234,12 +229,24 @@ subroutine coupled_em(&
  logical(lgt), parameter              :: printBalance=.false.   ! flag to print the balance checks
  real(rkind), allocatable                :: liqSnowInit(:)         ! volumetric liquid water conetnt of snow at the start of the time step
  real(rkind), allocatable                :: liqSoilInit(:)         ! soil moisture at the start of the time step
+ ! timing information
+ real(rkind)                             :: startTime              ! start time (used to compute wall clock time)
+ real(rkind)                             :: endTime                ! end time (used to compute wall clock time)
  ! ----------------------------------------------------------------------------------------------------------------------------------------------
  ! initialize error control
  err=0; message="coupled_em/"
 
  ! This is the start of a data step for a local HRU
 
+ ! get the start time
+ call cpu_time(startTime)
+
+ ! check the sundials decision
+ if(model_decisions(iLookDECISIONS%num_method)%iDecision==sundials)then
+  message=trim(message)//'still need to implement the sundials solver'
+  err=20; return
+ endif
+
  ! check that the decision is supported
  if(model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket .and. &
     model_decisions(iLookDECISIONS%spatial_gw)%iDecision/=localColumn)then
@@ -1213,6 +1220,12 @@ subroutine coupled_em(&
   err=20; return
  end if
 
+ ! get the end time
+ call cpu_time(endTime)
+
+ ! get the elapsed time
+ diag_data%var(iLookDIAG%wallClockTime)%dat(1) = endTime - startTime
+
  end subroutine coupled_em
 
 

build/source/engine/mDecisions.f90

@@ -57,9 +57,8 @@ module mDecisions_module
 integer(i4b),parameter,public :: constantScaling      =  71    ! constant scaling factor
 integer(i4b),parameter,public :: laiScaling           =  72    ! exponential function of LAI (Leuning, Plant Cell Env 1995: "Scaling from..." [eq 9])
 ! look-up values for the choice of numerical method
-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
+integer(i4b),parameter,public :: bEuler               =  81    ! home-grown backward Euler solution with long time steps
+integer(i4b),parameter,public :: sundials             =  82    ! SUNDIALS/IDA solution
 ! 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
@@ -397,9 +396,9 @@ subroutine mDecisions(err,message)
 
  ! identify the numerical method
  select case(trim(model_decisions(iLookDECISIONS%num_method)%cDecision))
-  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
  end select

build/source/engine/ssdNrgFlux.f90

@@ -64,10 +64,6 @@ module ssdNrgFlux_module
 
 ! provide access to look-up values for model decisions
 USE mDecisions_module,only:      &
- ! look-up values for the numerical method
- iterative,                      & ! iterative
- nonIterative,                   & ! non-iterative
- iterSurfEnergyBal,              & ! iterate only on the surface energy balance
  ! look-up values for method used to compute derivative
  numerical,                      & ! numerical solution
  analytical,                     & ! analytical solution

build/source/engine/sunGeomtry.f90

@@ -108,8 +108,10 @@ SUBROUTINE CLRSKY_RAD(MONTH,DAY,HOUR,DT,SLOPE,AZI,LAT,HRI,COSZEN)
  TPI=-TAN(LP)*TAN(D)
  IF(ABS(TPI).LT.1.0) THEN
   TP=ACOS(TPI)
- ELSE
-  TP=0.0
+ ELSE IF(TPI.LT.-1.0) THEN ! 24h daylight
+  TP=ACOS(-1.0)
+ ELSE IF(TPI.GT.1.0) THEN ! 24h dark
+  TP=ACOS(1.0)
  ENDIF
  ! Calculate time adjustment for ground slope, aspect and latitude (DDT = 0 for level surface)
  DDT=ATAN(SIN(AZI1)*SIN(SLOPE1)/(COS(SLOPE1)*COS(LAT1)-COS(AZI1)*SIN(SLOPE1)*SIN(LAT1)))
@@ -150,8 +152,10 @@ SUBROUTINE CLRSKY_RAD(MONTH,DAY,HOUR,DT,SLOPE,AZI,LAT,HRI,COSZEN)
   TPI=-TAN(LP)*TAN(D)
   IF(ABS(TPI).LT.1.0) THEN
    TP=ACOS(TPI)
-  ELSE
-   TP=0.0
+  ELSE IF(TPI.LT.-1.0) THEN ! 24h daylight
+   TP=ACOS(-1.0)
+  ELSE IF(TPI.GT.1.0) THEN ! 24h dark
+   TP=ACOS(1.0)
   ENDIF
   ! Set beginning time to sunrise
   T1=MAX(-TP-DDT,-TD)