Merge pull request #537 from CH-Earth/develop

Develop into master v3.2.0

.github/PULL_REQUEST_TEMPLATE.md

@@ -5,4 +5,4 @@ Make sure all the relevant boxes are checked (and only check the box if you actu
 - [ ] New tests added (describe which tests were performed to test the changes)
 - [ ] Science test figures (add figures to PR comment and describe the tests)
 - [ ] Checked that the new code conforms to the [SUMMA coding conventions](https://github.com/NCAR/summa/blob/master/docs/howto/summa_coding_conventions.md)
-- [ ] Describe the change in the release notes (use either `./summa/docs/whats-new.md` or `./summa/docs/minor-changes.md` depending on what changed)
+- [ ] Describe the change in the release notes (use  `./summa/docs/whats-new.md`)

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

@@ -143,7 +143,7 @@ subroutine flxMapping(err,message)
  flux2state_orig(iLookFLUX%mLayerLiqFluxSnow)               = flux2state(state1=iname_watLayer,  state2=integerMissing)
 
  ! liquid water fluxes for the soil domain
- flux2state_orig(iLookFLUX%scalarRainPlusMelt)              = flux2state(state1=iname_watLayer,  state2=integerMissing)
+ flux2state_orig(iLookFLUX%scalarRainPlusMelt)              = flux2state(state1=iname_matLayer,  state2=integerMissing)
  flux2state_orig(iLookFLUX%scalarMaxInfilRate)              = flux2state(state1=iname_matLayer,  state2=integerMissing)
  flux2state_orig(iLookFLUX%scalarInfiltration)              = flux2state(state1=iname_matLayer,  state2=integerMissing)
  flux2state_orig(iLookFLUX%scalarExfiltration)              = flux2state(state1=iname_matLayer,  state2=integerMissing)
@@ -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

@@ -197,8 +197,8 @@ MODULE globalData
  ! * part 2: globally constant variables/structures that require initialization
  ! ----------------------------------------------------------------------------------------------------------------
 
- ! define Indian bread (NaN)
- real(rkind),save,public                        :: dNaN
+ ! define Not-a-Number (NaN)
+ real(rkind),save,public                     :: dNaN
 
  ! define default parameter values and parameter bounds
  type(par_info),save,public                  :: localParFallback(maxvarMpar) ! local column default parameters
@@ -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

@@ -204,8 +204,8 @@ subroutine popMetadat(err,message)
  mpar_meta(iLookPARAM%rootDistExp)           = var_info('rootDistExp'           , 'exponent for the vertical distribution of root density'           , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%plantWiltPsi)          = var_info('plantWiltPsi'          , 'matric head at wilting point'                                     , 'm'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%soilStressParam)       = var_info('soilStressParam'       , 'parameter in the exponential soil stress function'                , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
- mpar_meta(iLookPARAM%critSoilWilting)       = var_info('critSoilWilting'       , 'critical vol. liq. water content when plants are wilting'         , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
- mpar_meta(iLookPARAM%critSoilTranspire)     = var_info('critSoilTranspire'     , 'critical vol. liq. water content when transpiration is limited'   , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
+ mpar_meta(iLookPARAM%critSoilWilting)       = var_info('critSoilWilting'       , 'critical vol. liq. water content when plants are wilting'         , '-'               , get_ixVarType('parSoil'), iMissVec, iMissVec, .false.)
+ mpar_meta(iLookPARAM%critSoilTranspire)     = var_info('critSoilTranspire'     , 'critical vol. liq. water content when transpiration is limited'   , '-'               , get_ixVarType('parSoil'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%critAquiferTranspire)  = var_info('critAquiferTranspire'  , 'critical aquifer storage value when transpiration is limited'     , 'm'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%minStomatalResistance) = var_info('minStomatalResistance' , 'minimum stomatal resistance'                                      , 's m-1'           , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%leafDimension)         = var_info('leafDimension'         , 'characteristic leaf dimension'                                    , 'm'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
@@ -241,7 +241,7 @@ subroutine popMetadat(err,message)
  ! scalar soil properties
  mpar_meta(iLookPARAM%fieldCapacity)         = var_info('fieldCapacity'         , 'soil field capacity (vol liq water content when baseflow begins)' , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%wettingFrontSuction)   = var_info('wettingFrontSuction'   , 'Green-Ampt wetting front suction'                                 , 'm'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
- mpar_meta(iLookPARAM%theta_mp)              = var_info('theta_mp'              , 'volumetric liquid water content when macropore flow begins'       , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
+ mpar_meta(iLookPARAM%theta_mp)              = var_info('theta_mp'              , 'volumetric liquid water content when macropore flow begins'       , '-'               , get_ixVarType('parSoil'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%mpExp)                 = var_info('mpExp'                 , 'empirical exponent in macropore flow equation'                    , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%kAnisotropic)          = var_info('kAnisotropic'          , 'anisotropy factor for lateral hydraulic conductivity'             , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
  mpar_meta(iLookPARAM%zScale_TOPMODEL)       = var_info('zScale_TOPMODEL'       , 'TOPMODEL scaling factor used in lower boundary condition for soil', 'm'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
@@ -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
@@ -274,7 +275,7 @@ subroutine check_icond(nGRU,                          & ! number of GRUs and HRU
     h1 = sum(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerDepth)%dat(1:iLayer)) ! sum of the depths up to the current layer
     h2 = progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(iLayer) - progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(0)  ! difference between snow-atm interface and bottom of layer
     if(abs(h1 - h2) > 1.e-6_rkind)then
-     write(message,'(a,1x,i0)') trim(message)//'mis-match between layer depth and layer height; layer = ', iLayer, '; sum depths = ',h1,'; height = ',h2
+     write(message,'(a,1x,i0,a,f5.3,a,f5.3)') trim(message)//'mis-match between layer depth and layer height; layer = ', iLayer, '; sum depths = ',h1,'; height = ',h2
      err=20; return
     end if
    end do