Merge pull request #452 from NCAR/develop

Bugfixes and minor improvements

build/Makefile

@@ -91,12 +91,12 @@ ifeq "$(FC)" "gfortran"
 	endif
 
 # Production runs
-FLAGS_NOAH = -O3 -ffree-form -fdefault-real-8 -ffree-line-length-none -fmax-errors=0 $(FLAGS_OMP)
+FLAGS_NOAH = -O3 -ffree-form -ffree-line-length-none -fmax-errors=0 $(FLAGS_OMP)
 FLAGS_COMM = -O3 -ffree-line-length-none -fmax-errors=0 $(FLAGS_OMP)
 FLAGS_SUMMA = -O3 -ffree-line-length-none -fmax-errors=0 $(FLAGS_OMP)
 
 # Debug runs
-#FLAGS_NOAH = -p -g -ffree-form -fdefault-real-8 -ffree-line-length-none -fmax-errors=0 -fbacktrace -Wno-unused -Wno-unused-dummy-argument
+#FLAGS_NOAH = -p -g -ffree-form -ffree-line-length-none -fmax-errors=0 -fbacktrace -Wno-unused -Wno-unused-dummy-argument
 #FLAGS_COMM = -p -g -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds
 #FLAGS_SUMMA = -p -g -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds
 
@@ -110,12 +110,12 @@ ifeq "$(FC)" "ifort"
 	endif
 
 # Production runs
-FLAGS_NOAH = -O3 -autodouble -noerror_limit -FR -auto -fltconsistency $(FLAGS_OMP)
+FLAGS_NOAH = -O3 -noerror_limit -FR -auto -fltconsistency $(FLAGS_OMP)
 FLAGS_COMM = -O3 -FR -auto -fltconsistency -fpe0 $(FLAGS_OMP)
 FLAGS_SUMMA = -O3 -FR -auto -fltconsistency -fpe0 $(FLAGS_OMP)
 
 # Debug runs
-#FLAGS_NOAH = -O0 -p -g -warn nounused -autodouble -noerror_limit -FR -auto -WB -traceback -fltconsistency
+#FLAGS_NOAH = -O0 -p -g -warn nounused -noerror_limit -FR -auto -WB -traceback -fltconsistency
 #FLAGS_COMM = -O0 -p -g -debug -warn all -check all -FR -auto -WB -traceback -fltconsistency -fpe0
 #FLAGS_SUMMA = -O0 -p -g -debug -warn all -check all -FR -auto -WB -traceback -fltconsistency -fpe0
 endif
@@ -357,7 +357,7 @@ update_version:
 
 # compile Noah-MP routines
 compile_noah:
-	$(FC_EXE) $(FLAGS_NOAH) -c $(NOAHMP)
+	$(FC_EXE) $(FLAGS_NOAH) -c $(NRUTIL) $(NOAHMP)
 
 # compile common routines
 compile_comm:

build/source/driver/summa_defineOutput.f90

@@ -100,11 +100,11 @@ subroutine summa_defineOutputFiles(modelTimeStep, summa1_struc, err, message)
  ! *** define the name of the model output file
  ! *****************************************************************************
 
- ! define full name of output file 
+ ! define full name of output file
  if(modelTimeStep==1)then
   select case(newOutputFile)
    case(noNewFiles);          ! do nothing, just ensure validity of outputfile option
-   case(newFileEveryOct1);    
+   case(newFileEveryOct1);
    case default; err=20; message=trim(message)//'unable to identify the option to define new output files'; return
   end select
 
@@ -141,7 +141,7 @@ subroutine summa_defineOutputFiles(modelTimeStep, summa1_struc, err, message)
   end do  ! (looping through HRUs)
 
   ! write GRU parameters
-  call writeParm(integerMissing,bparStruct%gru(iGRU),bpar_meta,err,cmessage)
+  call writeParm(iGRU,bparStruct%gru(iGRU),bpar_meta,err,cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
 
  end do  ! (looping through GRUs)

build/source/driver/summa_globalData.f90

@@ -107,7 +107,7 @@ subroutine summa_defineGlobalData(err, message)
  doJacobian=.false.        ! initialize the Jacobian flag
 
  ! define double precision NaNs (shared in globalData)
- dNaN = ieee_value(1._dp, ieee_quiet_nan)
+ dNaN = ieee_value(1._rkind, ieee_quiet_nan)
 
  ! populate metadata for all model variables
  call popMetadat(err,cmessage)

build/source/driver/summa_init.f90

@@ -175,9 +175,9 @@ subroutine summa_initialize(summa1_struc, err, message)
  ncid(:) = integerMissing
 
  ! initialize the elapsed time for cumulative quantities
- elapsedRead=0._dp
- elapsedWrite=0._dp
- elapsedPhysics=0._dp
+ elapsedRead=0._rkind
+ elapsedWrite=0._rkind
+ elapsedPhysics=0._rkind
 
  ! get the command line arguments
  call getCommandArguments(summa1_struc,err,cmessage)

build/source/driver/summa_modelRun.f90

@@ -72,16 +72,16 @@ subroutine summa_runPhysics(modelTimeStep, summa1_struc, err, message)
  integer(i4b)                          :: iGRU,jGRU,kGRU        ! GRU indices
  ! local variables: veg phenology
  logical(lgt)                          :: computeVegFluxFlag    ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
- real(dp)                              :: notUsed_canopyDepth   ! NOT USED: canopy depth (m)
- real(dp)                              :: notUsed_exposedVAI    ! NOT USED: exposed vegetation area index (m2 m-2)
+ real(rkind)                              :: notUsed_canopyDepth   ! NOT USED: canopy depth (m)
+ real(rkind)                              :: notUsed_exposedVAI    ! NOT USED: exposed vegetation area index (m2 m-2)
  ! local variables: parallelize the model run
  integer(i4b), allocatable             :: ixExpense(:)          ! ranked index GRU w.r.t. computational expense
  integer(i4b), allocatable             :: totalFluxCalls(:)     ! total number of flux calls for each GRU
  ! local variables: timing information
  integer*8                             :: openMPstart,openMPend ! time for the start of the parallelization section
  integer*8, allocatable                :: timeGRUstart(:)       ! time GRUs start
- real(dp),  allocatable                :: timeGRUcompleted(:)   ! time required to complete each GRU
- real(dp),  allocatable                :: timeGRU(:)            ! time spent on each GRU
+ real(rkind),  allocatable                :: timeGRUcompleted(:)   ! time required to complete each GRU
+ real(rkind),  allocatable                :: timeGRU(:)            ! time spent on each GRU
  ! ---------------------------------------------------------------------------------------
  ! associate to elements in the data structure
  summaVars: associate(&
@@ -171,7 +171,7 @@ subroutine summa_runPhysics(modelTimeStep, summa1_struc, err, message)
 
  ! compute the total number of flux calls from the previous time step
  do jGRU=1,nGRU
-  totalFluxCalls(jGRU) = 0._dp
+  totalFluxCalls(jGRU) = 0._rkind
   do iHRU=1,gru_struc(jGRU)%hruCount
    totalFluxCalls(jGRU) = totalFluxCalls(jGRU) + indxStruct%gru(jGRU)%hru(iHRU)%var(iLookINDEX%numberFluxCalc)%dat(1)
   end do
@@ -268,8 +268,8 @@ subroutine summa_runPhysics(modelTimeStep, summa1_struc, err, message)
   !$omp critical(saveTiming)
   ! save timing information
   call system_clock(openMPend)
-  timeGRU(iGRU)          = real(openMPend - timeGRUstart(iGRU), kind(dp))
-  timeGRUcompleted(iGRU) = real(openMPend - openMPstart       , kind(dp))
+  timeGRU(iGRU)          = real(openMPend - timeGRUstart(iGRU), kind(rkind))
+  timeGRUcompleted(iGRU) = real(openMPend - openMPstart       , kind(rkind))
   !$omp end critical(saveTiming)
 
  end do  ! (looping through GRUs)

build/source/driver/summa_restart.f90

@@ -178,7 +178,7 @@ subroutine summa_readRestart(summa1_struc, err, message)
 
    ! initialize canopy drip
    ! NOTE: canopy drip from the previous time step is used to compute throughfall for the current time step
-   fluxStruct%gru(iGRU)%hru(iHRU)%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp  ! not used
+   fluxStruct%gru(iGRU)%hru(iHRU)%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._rkind  ! not used
 
   end do  ! end looping through HRUs
 
@@ -201,14 +201,14 @@ subroutine summa_readRestart(summa1_struc, err, message)
 
    ! the basin-average aquifer storage is not used if the groundwater is included in the local column
    case(localColumn)
-    bvarStruct%gru(iGRU)%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 0._dp ! set to zero to be clear that there is no basin-average aquifer storage in this configuration
+    bvarStruct%gru(iGRU)%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 0._rkind ! set to zero to be clear that there is no basin-average aquifer storage in this configuration
 
    ! the local column aquifer storage is not used if the groundwater is basin-average
    ! (i.e., where multiple HRUs drain to a basin-average aquifer)
    case(singleBasin)
-    bvarStruct%gru(iGRU)%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 1._dp
+    bvarStruct%gru(iGRU)%var(iLookBVAR%basin__AquiferStorage)%dat(1) = 1._rkind
     do iHRU=1,gru_struc(iGRU)%hruCount
-     progStruct%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarAquiferStorage)%dat(1) = 0._dp  ! set to zero to be clear that there is no local aquifer storage in this configuration
+     progStruct%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarAquiferStorage)%dat(1) = 0._rkind  ! set to zero to be clear that there is no local aquifer storage in this configuration
     end do
 
    ! error check

build/source/driver/summa_setup.f90

@@ -191,7 +191,7 @@ subroutine summa_paramSetup(summa1_struc, err, message)
  ! *****************************************************************************
 
  ! define monthly fraction of green vegetation
- greenVegFrac_monthly = (/0.01_dp, 0.02_dp, 0.03_dp, 0.07_dp, 0.50_dp, 0.90_dp, 0.95_dp, 0.96_dp, 0.65_dp, 0.24_dp, 0.11_dp, 0.02_dp/)
+ greenVegFrac_monthly = (/0.01_rkind, 0.02_rkind, 0.03_rkind, 0.07_rkind, 0.50_rkind, 0.90_rkind, 0.95_rkind, 0.96_rkind, 0.65_rkind, 0.24_rkind, 0.11_rkind, 0.02_rkind/)
 
  ! read Noah soil and vegetation tables
  call soil_veg_gen_parm(trim(SETTINGS_PATH)//trim(VEGPARM),                            & ! filename for vegetation table
@@ -298,7 +298,7 @@ subroutine summa_paramSetup(summa1_struc, err, message)
 
   ! compute total area of the upstream HRUS that flow into each HRU
   do iHRU=1,gru_struc(iGRU)%hruCount
-   upArea%gru(iGRU)%hru(iHRU) = 0._dp
+   upArea%gru(iGRU)%hru(iHRU) = 0._rkind
    do jHRU=1,gru_struc(iGRU)%hruCount
     ! check if jHRU flows into iHRU; assume no exchange between GRUs
     if(typeStruct%gru(iGRU)%hru(jHRU)%var(iLookTYPE%downHRUindex)==typeStruct%gru(iGRU)%hru(iHRU)%var(iLookID%hruId))then
@@ -309,7 +309,7 @@ subroutine summa_paramSetup(summa1_struc, err, message)
 
   ! identify the total basin area for a GRU (m2)
   associate(totalArea => bvarStruct%gru(iGRU)%var(iLookBVAR%basin__totalArea)%dat(1) )
-  totalArea = 0._dp
+  totalArea = 0._rkind
   do iHRU=1,gru_struc(iGRU)%hruCount
    totalArea = totalArea + attrStruct%gru(iGRU)%hru(iHRU)%var(iLookATTR%HRUarea)
   end do

build/source/driver/summa_type.f90

@@ -91,11 +91,11 @@ MODULE summa_type
     ! define miscellaneous variables
     integer(i4b)                     :: summa1open                 ! flag to define if the summa file is open??
     integer(i4b)                     :: numout                     ! number of output variables??
-    real(dp)                         :: ts                         ! model time step ??
+    real(rkind)                         :: ts                         ! model time step ??
     integer(i4b)                     :: nGRU                       ! number of grouped response units
     integer(i4b)                     :: nHRU                       ! number of global hydrologic response units
     integer(i4b)                     :: hruCount                   ! number of local hydrologic response units
-    real(dp),dimension(12)           :: greenVegFrac_monthly       ! fraction of green vegetation in each month (0-1)
+    real(rkind),dimension(12)           :: greenVegFrac_monthly       ! fraction of green vegetation in each month (0-1)
     character(len=256)               :: summaFileManagerFile       ! path/name of file defining directories and files
 
 end type summa1_type_dec

build/source/driver/summa_util.f90

@@ -350,7 +350,7 @@ subroutine stop_program(err,message)
  integer(i4b)                       :: endModelRun(8)  ! final time
  integer(i4b)                       :: localErr        ! local error code
  integer(i4b)                       :: iFreq           ! loop through output frequencies
- real(dp)                           :: elpSec          ! elapsed seconds
+ real(rkind)                           :: elpSec          ! elapsed seconds
 
  ! close any remaining output files
  ! NOTE: use the direct NetCDF call with no error checking since the file may already be closed
@@ -392,9 +392,9 @@ subroutine stop_program(err,message)
 
  ! print total elapsed time
  write(outunit,"(/,A,1PG15.7,A)")                                             '     elapsed time = ', elpSec,                ' s'
- write(outunit,"(A,1PG15.7,A)")                                               '       or           ', elpSec/60_dp,          ' m'
- write(outunit,"(A,1PG15.7,A)")                                               '       or           ', elpSec/3600_dp,        ' h'
- write(outunit,"(A,1PG15.7,A/)")                                              '       or           ', elpSec/86400_dp,       ' d'
+ write(outunit,"(A,1PG15.7,A)")                                               '       or           ', elpSec/60_rkind,          ' m'
+ write(outunit,"(A,1PG15.7,A)")                                               '       or           ', elpSec/3600_rkind,        ' h'
+ write(outunit,"(A,1PG15.7,A/)")                                              '       or           ', elpSec/86400_rkind,       ' d'
 
  ! print the number of threads
  write(outunit,"(A,i10,/)")                                                   '   number threads = ', nThreads

build/source/dshare/data_types.f90

@@ -48,18 +48,18 @@ MODULE data_types
   integer(i4b),allocatable               :: var_ix(:)                ! index of each forcing data variable in the data structure
   integer(i4b),allocatable               :: data_id(:)               ! netcdf variable id for each forcing data variable
   character(len=256),allocatable         :: varName(:)               ! netcdf variable name for each forcing data variable
-  real(dp)                               :: firstJulDay              ! first julian day in forcing file
-  real(dp)                               :: convTime2Days            ! factor to convert time to days
+  real(rkind)                               :: firstJulDay              ! first julian day in forcing file
+  real(rkind)                               :: convTime2Days            ! factor to convert time to days
  end type file_info
 
  ! ***********************************************************************************************************
  ! Define metadata on model parameters
  ! ***********************************************************************************************************
  ! define a data type to store model parameter information
  type,public  :: par_info
-  real(dp)                               :: default_val              ! default parameter value
-  real(dp)                               :: lower_limit              ! lower bound
-  real(dp)                               :: upper_limit              ! upper bound
+  real(rkind)                               :: default_val              ! default parameter value
+  real(rkind)                               :: lower_limit              ! lower bound
+  real(rkind)                               :: upper_limit              ! upper bound
  endtype par_info
 
  ! ***********************************************************************************************************
@@ -131,7 +131,7 @@ MODULE data_types
  ! NOTE: use derived types here to facilitate adding the "variable" dimension
  ! ** double precision type
  type, public :: dlength
-  real(dp),allocatable                :: dat(:)    ! dat(:)
+  real(rkind),allocatable                :: dat(:)    ! dat(:)
  endtype dlength
  ! ** integer type (4 byte)
  type, public :: ilength
@@ -168,7 +168,7 @@ MODULE data_types
 
  ! ** double precision type of fixed length
  type, public :: var_d
-  real(dp),allocatable                :: var(:)    ! var(:)
+  real(rkind),allocatable                :: var(:)    ! var(:)
  endtype var_d
  ! ** integer type of fixed length (4 byte)
  type, public :: var_i
@@ -181,7 +181,7 @@ MODULE data_types
 
  ! ** double precision type of fixed length
  type, public :: hru_d
-  real(dp),allocatable                :: hru(:)    ! hru(:)
+  real(rkind),allocatable                :: hru(:)    ! hru(:)
  endtype hru_d
  ! ** integer type of fixed length (4 byte)
  type, public :: hru_i

build/source/dshare/get_ixname.f90

@@ -171,6 +171,7 @@ function get_ixAttr(varName)
   case('contourLength' ); get_ixAttr = iLookATTR%contourLength  ! length of contour at downslope edge of HRU (m)
   case('HRUarea'       ); get_ixAttr = iLookATTR%HRUarea        ! area of each HRU (m2)
   case('mHeight'       ); get_ixAttr = iLookATTR%mHeight        ! measurement height above bare ground (m)
+  case('aspect'        ); get_ixAttr = iLookATTR%aspect         ! azimuth in degrees East of North (degrees)
   ! get to here if cannot find the variable
   case default
    get_ixAttr = integerMissing
@@ -886,6 +887,8 @@ function get_ixbvar(varName)
   case('basin__AquiferRecharge'        ); get_ixbvar = iLookBVAR%basin__AquiferRecharge          ! recharge to the aquifer (m s-1)
   case('basin__AquiferBaseflow'        ); get_ixbvar = iLookBVAR%basin__AquiferBaseflow          ! baseflow from the aquifer (m s-1)
   case('basin__AquiferTranspire'       ); get_ixbvar = iLookBVAR%basin__AquiferTranspire         ! transpiration from the aquifer (m s-1)
+  case('basin__TotalRunoff'            ); get_ixbvar = iLookBVAR%basin__TotalRunoff              ! total runoff to channel from all active components (m s-1)
+  case('basin__SoilDrainage'           ); get_ixbvar = iLookBVAR%basin__SoilDrainage             ! soil drainage (m s-1)
   ! variables to compute runoff
   case('routingRunoffFuture'           ); get_ixbvar = iLookBVAR%routingRunoffFuture             ! runoff in future time steps (m s-1)
   case('routingFractionFuture'         ); get_ixbvar = iLookBVAR%routingFractionFuture           ! fraction of runoff in future time steps (-)