Fix inline initial condition interpolation

components/elm/src/data_types/TopounitDataType.F90

@@ -648,7 +648,7 @@ subroutine restart_top_es(this, bounds, ncid, flag)
 
     call restartvar(ncid=ncid, flag=flag, varname='TS_TOPO', xtype=ncd_double, &
         dim1name='topounit', long_name='surface radiative temperature', &
-        units='K', interpinic_flag='copy', readvar=readvar, data=this%t_rad)
+        units='K', interpinic_flag='skip', readvar=readvar, data=this%t_rad)
 
   end subroutine restart_top_es
 
@@ -696,7 +696,7 @@ subroutine restart_top_ws(this, bounds, ncid, flag)
 
     call restartvar(ncid=ncid, flag=flag, varname='FROM_UPHILL', xtype=ncd_double, &
         dim1name='topounit', long_name='water received from uphill topounit(s)', &
-        units='kg/m2', interpinic_flag='interp', readvar=readvar, data=this%from_uphill)
+        units='kg/m2', interpinic_flag='skip', readvar=readvar, data=this%from_uphill)
     if (flag=='read' .and. .not. readvar) then
        this%from_uphill(bounds%begt:bounds%endt) = 0.0_r8
     end if

components/elm/src/main/accumulMod.F90

@@ -278,7 +278,7 @@ subroutine extract_accum_field_sl (name, field, nstep)
 
     beg = accum(nf)%beg1d
     end = accum(nf)%end1d
-    if (size(field,dim=1) < end-beg+1) then
+    if (size(field,dim=1) /= end-beg+1) then
        write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name
        write(iulog,*)'size of first dimension of field is ',&
             size(field,dim=1),' and should be ',end-beg+1
@@ -341,7 +341,7 @@ subroutine extract_accum_field_ml (name, field, nstep)
     numlev = accum(nf)%numlev
     beg = accum(nf)%beg1d
     end = accum(nf)%end1d
-    if (size(field,dim=1) < end-beg+1) then
+    if (size(field,dim=1) /= end-beg+1) then
        write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name
        write(iulog,*)'size of first dimension of field is ',&
             size(field,dim=1),' and should be ',end-beg+1
@@ -406,7 +406,7 @@ subroutine update_accum_field_sl (name, field, nstep)
 
     beg = accum(nf)%beg1d
     end = accum(nf)%end1d
-    if (size(field,dim=1) < end-beg+1) then
+    if (size(field,dim=1) /= end-beg+1) then
        write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_SL for field ',accum(nf)%name
        write(iulog,*)'size of first dimension of field is ',size(field,dim=1),&
             ' and should be ',end-beg+1
@@ -500,7 +500,7 @@ subroutine update_accum_field_ml (name, field, nstep)
     numlev = accum(nf)%numlev
     beg = accum(nf)%beg1d
     end = accum(nf)%end1d
-    if (size(field,dim=1) < end-beg+1) then
+    if (size(field,dim=1) /= end-beg+1) then
        write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_ML for field ',accum(nf)%name
        write(iulog,*)'size of first dimension of field is ',size(field,dim=1),&
             ' and should be ',end-beg+1
@@ -588,6 +588,7 @@ subroutine accumulRest( ncid, flag )
     ! !LOCAL VARIABLES:
     integer :: nf                             ! indices
     logical :: readvar                        ! determine if variable is on initial file
+    character(len=  8) :: value_interpinic_flag ! initial condition interpolation policy    
     character(len=128) :: varname             ! temporary
     character(len= 32) :: subname='AccumRest' ! subroutine name
     !------------------------------------------------------------------------
@@ -598,17 +599,22 @@ subroutine accumulRest( ncid, flag )
        ! accum(nf)%val is always 2d
 
        varname = trim(accum(nf)%name) // '_VALUE'
+       if (trim(accum(nf)%type1d) == 'topounit') then
+          value_interpinic_flag = 'skip'
+       else
+          value_interpinic_flag = 'interp'
+       end if       
        if (accum(nf)%numlev == 1) then
           call restartvar(ncid=ncid, flag=flag, varname=varname, xtype=ncd_double, &
                dim1name=accum(nf)%type1d, &
                long_name=accum(nf)%desc, units=accum(nf)%units, &
-               interpinic_flag='interp', &
+               interpinic_flag=value_interpinic_flag, &
                data=accum(nf)%val, readvar=readvar)
        else
           call restartvar(ncid=ncid, flag=flag, varname=varname, xtype=ncd_double, &
                dim1name=accum(nf)%type1d, dim2name=accum(nf)%type2d, &
                long_name=accum(nf)%desc, units=accum(nf)%units, &
-               interpinic_flag='interp', &
+               interpinic_flag=value_interpinic_flag, &
                data=accum(nf)%val, readvar=readvar)
        end if
 

components/elm/src/main/initInterp.F90

@@ -120,7 +120,7 @@ subroutine initInterp (filei, fileo, bounds)
     integer , pointer  :: pftindx(:)     
     integer , pointer  :: colindx(:)     
     integer , pointer  :: lunindx(:)     
-    integer , pointer  :: grcindx(:) 
+    integer , pointer  :: grcindx(:)
     logical , pointer  :: pft_activei(:), pft_activeo(:) 
     logical , pointer  :: col_activei(:), col_activeo(:) 
     logical , pointer  :: lun_activei(:), lun_activeo(:) 
@@ -149,7 +149,6 @@ subroutine initInterp (filei, fileo, bounds)
     call check_dim_subgrid(ncidi, ncido, dimname ='column'  , dimleni=ncolsi, dimleno=ncolso)
     call check_dim_subgrid(ncidi, ncido, dimname ='landunit', dimleni=nlunsi, dimleno=nlunso)
     call check_dim_subgrid(ncidi, ncido, dimname ='gridcell', dimleni=ngrcsi, dimleno=ngrcso)
-
     if (masterproc) then
        write (iulog,*) 'input gridcells = ',ngrcsi,' output gridcells = ',ngrcso
        write (iulog,*) 'input landuntis = ',nlunsi,' output landunits = ',nlunso
@@ -196,27 +195,26 @@ subroutine initInterp (filei, fileo, bounds)
     begc_i = 1 ;  endc_i = ncolsi
     begl_i = 1 ;  endl_i = nlunsi
     begg_i = 1 ;  endg_i = ngrcsi
-
+    
     begp_o = bounds%begp ;  endp_o = bounds%endp
     begc_o = bounds%begc ;  endc_o = bounds%endc
     begl_o = bounds%begl ;  endl_o = bounds%endl
     begg_o = bounds%begg ;  endg_o = bounds%endg
-
+    
     allocate(pft_activei(begp_i:endp_i))
     allocate(col_activei(begc_i:endc_i))
     allocate(lun_activei(begl_i:endl_i))
     allocate(grc_activei(begg_i:endg_i))
-
+    
     allocate(pft_activeo(begp_o:endp_o))
     allocate(col_activeo(begc_o:endc_o))
     allocate(lun_activeo(begl_o:endl_o))
     allocate(grc_activeo(begg_o:endg_o))
-
+    
     allocate(pftindx(begp_o:endp_o))
     allocate(colindx(begc_o:endc_o))
     allocate(lunindx(begl_o:endl_o))
     allocate(grcindx(begg_o:endg_o))
-
     ! For each output pft, find the input pft, pftindx, that is closest
 
     if (masterproc) then
@@ -225,7 +223,6 @@ subroutine initInterp (filei, fileo, bounds)
     vec_dimname = 'pft'
     call findMinDist(vec_dimname, begp_i, endp_i, begp_o, endp_o, ncidi, ncido, &
          pft_activei, pft_activeo, pftindx )
-
     ! For each output column, find the input column, colindx, that is closest
 
     if (masterproc) then
@@ -377,6 +374,16 @@ subroutine initInterp (filei, fileo, bounds)
        else if ( ndimso == 1 ) then
 
           status = pio_inq_dimname(ncido, dimidso(1), vec_dimname)
+          if (trim(vec_dimname) == 'topounit') then
+             if (masterproc) then
+               call endrun(msg='ERROR interpinic: 1D topounit variable '//trim(varname)// &
+                  ' is not marked skip, but topounit interpolation is not supported. '// &
+                  'Set interpinic_flag=skip or implement a valid topounit mapping.'// &
+                  errMsg(__FILE__, __LINE__))
+             end if
+             CYCLE
+          end if
+
           if ( vec_dimname == 'pft' )then
              begi = begp_i
              endi = endp_i
@@ -708,6 +715,7 @@ subroutine set_subgrid_glob(beg, end, dimname, ncid, active, subgrid)
     else if (dimname == 'gridcell') then
        call ncd_io(ncid=ncid, varname='grid1d_lon'    , flag='read', data=subgrid%lon  )
        call ncd_io(ncid=ncid, varname='grid1d_lat'    , flag='read', data=subgrid%lat  )
+       itemp(:)=1
     end if
 
     do n = beg,end