merged in prms 5.3.0

GSFLOW/src/prms/gwflow.f90

@@ -17,7 +17,7 @@ MODULE PRMS_GWFLOW
 !   Local Variables
       character(len=*), parameter :: MODDESC = 'Groundwater'
       character(len=6), parameter :: MODNAME = 'gwflow'
-      character(len=*), parameter :: Version_gwflow = '2023-11-13'
+      character(len=*), parameter :: Version_gwflow = '2024-01-22'
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gwstor_minarea(:), Gwin_dprst(:)
       DOUBLE PRECISION, SAVE :: Basin_gw_upslope
       INTEGER, SAVE :: Gwminarea_flag
@@ -29,7 +29,7 @@ MODULE PRMS_GWFLOW
       DOUBLE PRECISION, SAVE :: Basin_gwstor_minarea_wb
       REAL, SAVE, ALLOCATABLE :: Gwres_flow(:), Gwres_sink(:)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gw_upslope(:), Gwres_in(:)
-      REAL, SAVE, ALLOCATABLE :: Hru_gw_cascadeflow(:)
+      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Hru_gw_cascadeflow(:)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gw_in_soil(:), Gw_in_ssr(:), Hru_lateral_flow(:)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Gwstor_minarea_wb(:), Hru_streamflow_out(:), Lakein_gwflow(:)
       DOUBLE PRECISION, SAVE, ALLOCATABLE :: Lake_seepage(:), Gw_seep_lakein(:), Lake_seepage_gwr(:)
@@ -93,7 +93,7 @@ INTEGER FUNCTION gwflowdecl()
      &       'acre-inches', Gw_upslope)
 
         ALLOCATE ( Hru_gw_cascadeflow(Ngw) )
-        CALL declvar_real(MODNAME, 'hru_gw_cascadeflow', 'ngw', Ngw, &
+        CALL declvar_dble(MODNAME, 'hru_gw_cascadeflow', 'ngw', Ngw, &
      &       'Cascading groundwater flow from each GWR', &
      &       'inches', Hru_gw_cascadeflow)
 
@@ -269,7 +269,7 @@ INTEGER FUNCTION gwflowinit()
       USE PRMS_MODULE, ONLY: Ngw, Nlake, Print_debug, Init_vars_from_file, &
      &    Dprst_flag, Inputerror_flag, Gwr_swale_flag
       USE PRMS_GWFLOW
-      USE PRMS_BASIN, ONLY: Gwr_type, Hru_area, Basin_area_inv, Active_gwrs, Gwr_route_order, &
+      USE PRMS_BASIN, ONLY: Gwr_type, Hru_area_dble, Basin_area_inv, Active_gwrs, Gwr_route_order, &
      &                      Lake_hru_id, Weir_gate_flag, Hru_storage
       USE PRMS_FLOWVARS, ONLY: Gwres_stor
       use prms_utils, only: read_error
@@ -296,10 +296,10 @@ INTEGER FUNCTION gwflowinit()
       Basin_gwstor = 0.0D0
       DO j = 1, Active_gwrs
         i = Gwr_route_order(j)
-        Basin_gwstor = Basin_gwstor + Gwres_stor(i)*DBLE(Hru_area(i))
+        Basin_gwstor = Basin_gwstor + Gwres_stor(i) * Hru_area_dble(i)
         IF ( Gwstor_min(i)>0.0 ) THEN
           Gwminarea_flag = 1
-          Gwstor_minarea(i) = DBLE( Gwstor_min(i)*Hru_area(i) )
+          Gwstor_minarea(i) = DBLE( Gwstor_min(i) ) * Hru_area_dble(i)
         ENDIF
         IF ( Gwflow_coef(i)>1.0 ) THEN
           IF ( Print_debug>DEBUG_less ) PRINT *, 'WARNING, gwflow_coef value > 1.0 for GWR:', i, Gwflow_coef(i)
@@ -377,7 +377,7 @@ INTEGER FUNCTION gwflowrun()
      &    Gwr_add_water_use, Gwr_transfer_water_use, Nowyear, Nowmonth, Nowday
       USE PRMS_GWFLOW
       USE PRMS_BASIN, ONLY: Active_gwrs, Gwr_route_order, Lake_type, &
-     &    Basin_area_inv, Hru_area, Gwr_type, Lake_hru_id, Weir_gate_flag, Hru_area_dble, Hru_storage
+     &    Basin_area_inv, Gwr_type, Lake_hru_id, Weir_gate_flag, Hru_area_dble, Hru_storage
       USE PRMS_FLOWVARS, ONLY: Soil_to_gw, Ssr_to_gw, Sroff, Ssres_flow, Gwres_stor, Lake_vol
       USE PRMS_CASCADE, ONLY: Ncascade_gwr
       USE PRMS_SET_TIME, ONLY: Cfs_conv
@@ -390,14 +390,14 @@ INTEGER FUNCTION gwflowrun()
       INTRINSIC :: DBLE, DABS, SNGL, MIN
 ! Local Variables
       INTEGER :: i, j, jj, jjj
-      REAL :: dnflow
+      DOUBLE PRECISION :: dnflow
       DOUBLE PRECISION :: gwin, gwstor, gwsink, gwflow, gwstor_last, seepage, gwarea, inch2acre_feet
 !***********************************************************************
       gwflowrun = 0
 
       IF ( Cascadegw_flag>CASCADEGW_OFF ) THEN
         Gw_upslope = 0.0D0
-        Hru_gw_cascadeflow = 0.0
+        Hru_gw_cascadeflow = 0.0D0
         Basin_dnflow = 0.0D0
         Basin_gw_upslope = 0.0D0
         IF ( Nlake>0 ) Lakein_gwflow = 0.0D0
@@ -466,8 +466,8 @@ INTEGER FUNCTION gwflowrun()
         gwarea = Hru_area_dble(i)
         gwstor = Gwres_stor(i)*gwarea ! acre-inches
         ! soil_to_gw is for whole HRU, not just perv
-        Gw_in_soil(i) = Soil_to_gw(i)*Hru_area(i)
-        Gw_in_ssr(i) = Ssr_to_gw(i)*Hru_area(i)
+        Gw_in_soil(i) = DBLE( Soil_to_gw(i) ) * gwarea
+        Gw_in_ssr(i) = DBLE( Ssr_to_gw(i) ) * gwarea
         gwin = Gw_in_soil(i) + Gw_in_ssr(i)
         IF ( Cascadegw_flag>CASCADEGW_OFF ) THEN
           gwin = gwin + Gw_upslope(i)
@@ -518,29 +518,32 @@ INTEGER FUNCTION gwflowrun()
 
         gwsink = 0.0D0
         IF ( gwstor<0.0D0 ) THEN ! could happen with water use
-          IF ( Print_debug>DEBUG_less ) PRINT *, 'Warning, groundwater reservoir for HRU:', i, ' is < 0.0, set to 0.0', gwstor
-          gwstor = 0.0D0
-        ENDIF
+          IF ( Print_debug>DEBUG_less ) PRINT *, 'Warning, groundwater reservoir for HRU:', i, ' is < 0.0', gwstor
+          gwflow = 0.0D0
+          Gwres_sink(i) = 0.0
+        ELSE
+
 ! Compute groundwater discharge
-        gwflow = gwstor*DBLE( Gwflow_coef(i) )
+          gwflow = gwstor*DBLE( Gwflow_coef(i) )
 
 ! Reduce storage by outflow
-        gwstor = gwstor - gwflow
+          gwstor = gwstor - gwflow
 
-        IF ( Gwsink_coef(i)>0.0 ) THEN
-          gwsink = MIN( gwstor*DBLE( Gwsink_coef(i) ), gwstor ) ! if gwsink_coef > 1, could have had negative gwstor
-          gwstor = gwstor - gwsink
-        ENDIF
+          IF ( Gwsink_coef(i)>0.0 ) THEN
+            gwsink = MIN( gwstor*DBLE( Gwsink_coef(i) ), gwstor ) ! if gwsink_coef > 1, could have had negative gwstor
+            gwstor = gwstor - gwsink
+          ENDIF
 ! if gwr_swale_flag = 1 swale GWR flow goes to sink, 2 included in stream network and cascades
 ! maybe gwr_swale_flag = 3 abs(hru_segment) so hru_segment could be changed from 0 to allow HRU swales
-        IF ( Gwr_swale_flag==ACTIVE ) THEN
-          IF ( Gwr_type(i)==SWALE ) THEN
-            gwsink = gwsink + gwflow
-            gwflow = 0.0D0
+          IF ( Gwr_swale_flag==ACTIVE ) THEN
+            IF ( Gwr_type(i)==SWALE ) THEN
+              gwsink = gwsink + gwflow
+              gwflow = 0.0D0
+            ENDIF
           ENDIF
+          Gwres_sink(i) = SNGL( gwsink/gwarea )
+          Basin_gwsink = Basin_gwsink + gwsink
         ENDIF
-        Gwres_sink(i) = SNGL( gwsink/gwarea )
-        Basin_gwsink = Basin_gwsink + gwsink
         Basin_gwstor = Basin_gwstor + gwstor
 
         Gwres_flow(i) = SNGL( gwflow/gwarea )
@@ -589,19 +592,19 @@ SUBROUTINE rungw_cascade(Igwr, Ncascade_gwr, Gwres_flow, Dnflow)
 ! Arguments
       INTEGER, INTENT(IN) :: Igwr, Ncascade_gwr
       REAL, INTENT(INOUT) :: Gwres_flow
-      REAL, INTENT(OUT) :: Dnflow
+      DOUBLE PRECISION, INTENT(OUT) :: Dnflow
 ! Local variables
       INTEGER :: j, k
       DOUBLE PRECISION :: gwflow_in
 !***********************************************************************
-      Dnflow = 0.0
+      Dnflow = 0.0D0
       DO k = 1, Ncascade_gwr
         j = Gwr_down(k, Igwr)
         ! Gwres_flow is in inches
 ! if gwr_down(k, Igwr) > 0, cascade contributes to a downslope GWR
         IF ( j>0 ) THEN
-          Gw_upslope(j) = Gw_upslope(j) + Gwres_flow*Cascade_gwr_area(k, Igwr)
-          Dnflow = Dnflow + Gwres_flow*Gwr_down_frac(k, Igwr)
+          Gw_upslope(j) = Gw_upslope(j) + DBLE( Gwres_flow*Cascade_gwr_area(k, Igwr) )
+          Dnflow = Dnflow + DBLE( Gwres_flow*Gwr_down_frac(k, Igwr) )
 ! if gwr_down(k, Igwr) < 0, cascade contributes to a stream
         ELSEIF ( j<0 ) THEN
           j = IABS( j )
@@ -612,7 +615,7 @@ SUBROUTINE rungw_cascade(Igwr, Ncascade_gwr, Gwres_flow, Dnflow)
       ENDDO
 
       ! gwres_flow reduced by cascading flow to HRUs
-      Gwres_flow = Gwres_flow - Dnflow
+      Gwres_flow = Gwres_flow - SNGL( Dnflow )
       IF ( Gwres_flow<0.0 ) Gwres_flow = 0.0
 
       END SUBROUTINE rungw_cascade

GSFLOW/src/prms/nhru_summary.f90

@@ -132,6 +132,7 @@ SUBROUTINE nhru_summaryinit()
 ! Local Variables
       INTEGER :: ios, ierr, size, jj, values_lastrow
       CHARACTER(LEN=MAXFILE_LENGTH) :: fileName
+      CHARACTER(LEN=4) ::file_suffix
 !***********************************************************************
       Begin_results = ACTIVE
       IF ( Prms_warmup>0 ) Begin_results = OFF
@@ -244,20 +245,23 @@ SUBROUTINE nhru_summaryinit()
       WRITE ( Output_fmt2, 9002 ) Nhru
       ALLOCATE ( Nhru_var_daily(Nhru, NhruOutVars) )
       Nhru_var_daily = 0.0
+      file_suffix = '.csv'
+      IF ( write_binary_cbh == ACTIVE ) file_suffix = '.bin'
+      print *, write_binary_cbh, file_suffix
       DO jj = 1, NhruOutVars
         IF ( Daily_flag==ACTIVE ) THEN
-          fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'.csv'
+          fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//file_suffix
           CALL PRMS_open_output_file(Dailyunit(jj), fileName, 'xxx', write_binary_cbh, ios)
           IF ( ios/=0 ) CALL error_stop('in nhru_summary, daily', ERROR_open_out)
           CALL write_header_date( Dailyunit(jj) )
         ENDIF
         IF ( NhruOut_freq>MEAN_MONTHLY ) THEN
           IF ( NhruOut_freq==MEAN_YEARLY ) THEN
-            fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'_meanyearly.csv'
+            fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'_meanyearly'//file_suffix
             CALL PRMS_open_output_file(Yearlyunit(jj), fileName, 'xxx', write_binary_cbh, ios)
             IF ( ios/=0 ) CALL error_stop('in nhru_summary, mean yearly', ERROR_open_out)
           ELSE  !IF ( NhruOut_freq==YEARLY ) THEN
-            fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'_yearly.csv'
+            fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'_yearly'//file_suffix
             CALL PRMS_open_output_file(Yearlyunit(jj), fileName, 'xxx', write_binary_cbh, ios)
             IF ( ios/=0 ) CALL error_stop('in nhru_summary, yearly', ERROR_open_out)
           ENDIF
@@ -270,7 +274,7 @@ SUBROUTINE nhru_summaryinit()
             CALL PRMS_open_output_file(Monthlyunit(jj), fileName, 'xxx', write_binary_cbh, ios)
             IF ( ios/=0 ) CALL error_stop('in nhru_summary, mean monthly', ERROR_open_out)
           ELSE
-            fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'_monthly.csv'
+            fileName = NhruOutBaseFileName(:numchars(NhruOutBaseFileName))//NhruOutVar_names(jj)(:Nc_vars(jj))//'_monthly'//file_suffix
             CALL PRMS_open_output_file(Monthlyunit(jj), fileName, 'xxx', write_binary_cbh, ios)
             IF ( ios/=0 ) CALL error_stop('in nhru_summary, monthly', ERROR_open_out)
           ENDIF
@@ -472,13 +476,13 @@ SUBROUTINE write_header_date(Iunit)
         IF ( write_binary_cbh == OFF ) THEN
           WRITE ( Iunit, Output_fmt2 ) (j, j=1,Nhru)
         ELSE
-          WRITE ( Iunit ) 'Date', (', ', j, j=1,Nhru)
+          WRITE ( Iunit ) 'Date', (j, j=1,Nhru)
         ENDIF
       ELSE
         IF ( write_binary_cbh == OFF ) THEN
           WRITE ( Iunit, Output_fmt2 ) (Nhm_id(j), j=1,Nhru)
         ELSE
-          WRITE ( Iunit ) 'Date', (', ', Nhm_id(j), j=1,Nhru)
+          WRITE ( Iunit ) 'Date', (Nhm_id(j), j=1,Nhru)
         ENDIF
       ENDIF
 
@@ -515,12 +519,11 @@ SUBROUTINE write_CBH_yearly_values( ivar )
           ELSE
             WRITE ( Yearlyunit(ivar) ) Lastyear
             first = 1
-            last = 0
             DO jj = 1, nrows
-              last = last + NhruOutNcol
+              last = first + NhruOutNcol
               IF ( last > Nhru ) last = Nhru
-              WRITE ( Yearlyunit(ivar) ) Nhru_var_yearly(first,ivar), (',', Nhru_var_yearly(j,ivar), j=first+1,last)
-              first = first + NhruOutNcol
+              WRITE ( Yearlyunit(ivar) ) (Nhru_var_yearly(j,ivar), j=first,last)
+              first = last + 1
             ENDDO
           ENDIF
         ENDIF
@@ -532,7 +535,7 @@ SUBROUTINE write_CBH_yearly_values( ivar )
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Yearlyunit(ivar), Output_fmt3int ) Lastyear, (Nhru_var_int(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Yearlyunit(ivar) ) Lastyear, (',', Nhru_var_int(j,ivar), j=1,Nhru)
+            WRITE ( Yearlyunit(ivar) ) Lastyear, (Nhru_var_int(j,ivar), j=1,Nhru)
           ENDIF
         ELSE
           IF ( write_binary_cbh == OFF ) THEN
@@ -541,12 +544,11 @@ SUBROUTINE write_CBH_yearly_values( ivar )
           ELSE
             WRITE ( Yearlyunit(ivar) ) Lastyear
             first = 1
-            last = 0
             DO jj = 1, nrows
-              last = last + NhruOutNcol
+              last = first + NhruOutNcol
               IF ( last > Nhru ) last = Nhru
-              WRITE ( Yearlyunit(ivar) ) Nhru_var_int(first,ivar), (',', Nhru_var_int(j,ivar), j=first+1, last)
-              first = first + NhruOutNcol
+              WRITE ( Yearlyunit(ivar) ) (Nhru_var_int(j,ivar), j=first, last)
+              first = last + 1
             ENDDO
           ENDIF
         ENDIF
@@ -576,21 +578,20 @@ SUBROUTINE write_cbh_monthly_values( ivar )
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Monthlyunit(ivar), Output_fmt) Nowyear, Nowmonth, Nowday, (Nhru_var_monthly(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Monthlyunit(ivar) ) Nowyear, '-', Nowmonth, '-', Nowday, (',', Nhru_var_monthly(j,ivar), j=1,Nhru)
+            WRITE ( Monthlyunit(ivar) ) Nowyear, Nowmonth, Nowday, (Nhru_var_monthly(j,ivar), j=1,Nhru)
           ENDIF
         ELSE
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Monthlyunit(ivar), Output_date_fmt) Nowyear, Nowmonth, Nowday
             WRITE ( Monthlyunit(ivar), Output_grid_fmt) (Nhru_var_monthly(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Monthlyunit(ivar) ) Nowyear, '-', Nowmonth, '-', Nowday
+            WRITE ( Monthlyunit(ivar) ) Nowyear, Nowmonth, Nowday
             first = 1
-            last = 0
             DO jj = 1, nrows
-              last = last + NhruOutNcol
+              last = first + NhruOutNcol
               IF ( last > Nhru ) last = Nhru
-              WRITE ( Monthlyunit(ivar) ) Nhru_var_monthly(first,ivar), (',', Nhru_var_monthly(j,ivar), j=first+1,last)
-              first = first + NhruOutNcol
+              WRITE ( Monthlyunit(ivar) ) (Nhru_var_monthly(j,ivar), j=first,last)
+              first = last + 1
             ENDDO
           ENDIF
         ENDIF
@@ -602,7 +603,7 @@ SUBROUTINE write_cbh_monthly_values( ivar )
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Monthlyunit(ivar), Output_fmtint) Nowyear, Nowmonth, Nowday, (Nhru_var_int(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Monthlyunit(ivar) ) Nowyear, '-', Nowmonth, '-', Nowday, (',', Nhru_var_int(j,ivar), j=1,Nhru)
+            WRITE ( Monthlyunit(ivar) ) Nowyear, Nowmonth, Nowday, (Nhru_var_int(j,ivar), j=1,Nhru)
           ENDIF
         ELSE
           IF ( write_binary_cbh == OFF ) THEN
@@ -611,12 +612,11 @@ SUBROUTINE write_cbh_monthly_values( ivar )
           ELSE
             WRITE ( Monthlyunit(jj) ) Nowyear, Nowmonth, Nowday
             first = 1
-            last = 0
             DO jj = 1, nrows
-              last = last + NhruOutNcol
+              last = first + NhruOutNcol
               IF ( last > Nhru ) last = Nhru
-              WRITE ( Monthlyunit(ivar) ) Nhru_var_int(first,ivar), (',', Nhru_var_int(j,ivar), j=first+1,last)
-              first = first + NhruOutNcol
+              WRITE ( Monthlyunit(ivar) ) (Nhru_var_int(j,ivar), j=first,last)
+              first = last + 1
             ENDDO
           ENDIF
         ENDIF
@@ -645,21 +645,20 @@ SUBROUTINE write_cbh_daily_values( ivar )
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Dailyunit(ivar), Output_fmt) Nowyear, Nowmonth, Nowday, (Nhru_var_daily(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Dailyunit(ivar) ) Nowyear, '-', Nowmonth, '-', Nowday, (',', Nhru_var_daily(j,ivar), j=1,Nhru)
+            WRITE ( Dailyunit(ivar) ) Nowyear, Nowmonth, Nowday, (Nhru_var_daily(j,ivar), j=1,Nhru)
           ENDIF
         ELSE
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Dailyunit(ivar), Output_date_fmt) Nowyear, Nowmonth, Nowday
             WRITE ( Dailyunit(ivar), Output_grid_fmt) (Nhru_var_daily(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Dailyunit(ivar) ) Nowyear, '-', Nowmonth, '-', Nowday
+            WRITE ( Dailyunit(ivar) ) Nowyear, Nowmonth, Nowday
             first = 1
-            last = 0
             DO jj = 1, nrows
-              last = last + NhruOutNcol
+              last = first + NhruOutNcol
               IF ( last > Nhru ) last = Nhru
-              WRITE ( Dailyunit(ivar) ) Nhru_var_daily(first,ivar), (',', Nhru_var_daily(j,ivar), j=first+1,last)
-              first = first + NhruOutNcol
+              WRITE ( Dailyunit(ivar) ) (Nhru_var_daily(j,ivar), j=first,last)
+              first = last + 1
             ENDDO
           ENDIF
         ENDIF
@@ -668,7 +667,7 @@ SUBROUTINE write_cbh_daily_values( ivar )
           IF ( write_binary_cbh == OFF ) THEN
             WRITE ( Dailyunit(ivar), Output_fmtint) Nowyear, Nowmonth, Nowday, (Nhru_var_int(j,ivar), j=1,Nhru)
           ELSE
-            WRITE ( Dailyunit(ivar) ) Nowyear, '-', Nowmonth, '-', Nowday, (',', Nhru_var_int(j,ivar), j=1,Nhru)
+            WRITE ( Dailyunit(ivar) ) Nowyear, Nowmonth, Nowday, (Nhru_var_int(j,ivar), j=1,Nhru)
           ENDIF
         ELSE
           IF ( write_binary_cbh == OFF ) THEN
@@ -677,12 +676,11 @@ SUBROUTINE write_cbh_daily_values( ivar )
           ELSE
             WRITE ( Dailyunit(ivar) ) Nowyear, Nowmonth, Nowday
             first = 1
-            last = 0
             DO jj = 1, nrows
-              last = last + NhruOutNcol
+              last = first + NhruOutNcol
               IF ( last > Nhru ) last = Nhru
-              WRITE ( Dailyunit(ivar) ) Nhru_var_int(1,ivar), (',', Nhru_var_int(j,ivar), j=first+1,last)
-              first = first + NhruOutNcol
+              WRITE ( Dailyunit(ivar) ) (Nhru_var_int(j,ivar), j=first,last)
+              first = last + 1
             ENDDO
           ENDIF
         ENDIF