From 39ba6895fc5936528fbbb420b0cba4de9b1f8d92 Mon Sep 17 00:00:00 2001 From: jinyuntang Date: Mon, 24 Aug 2026 12:17:29 -0700 Subject: [PATCH 1/3] fixed a bug in initInterp.F90 for initial field interpolation By avoiding some topo variables, it ensures the initial condition interpolation can move on. [BFB] --- components/elm/src/main/initInterp.F90 | 69 ++++++++++++++------------ 1 file changed, 37 insertions(+), 32 deletions(-) diff --git a/components/elm/src/main/initInterp.F90 b/components/elm/src/main/initInterp.F90 index a656f6338601..a62c6aba88fe 100644 --- a/components/elm/src/main/initInterp.F90 +++ b/components/elm/src/main/initInterp.F90 @@ -90,8 +90,8 @@ subroutine initInterp (filei, fileo, bounds) integer :: begc_o, endc_o ! output file column bounds integer :: begl_i, endl_i ! input file landunit bounds integer :: begl_o, endl_o ! output file landunit bounds - integer :: begg_i, endg_i ! input file gridcell bounds - integer :: begg_o, endg_o ! output file gridcell bounds + integer :: begg_i, endg_i ! input file gridcell bounds + integer :: begg_o, endg_o ! output file gridcell bounds integer :: nlevi,nlevo ! input/output number of levels type(file_desc_t) :: ncidi, ncido ! input/output pio fileids integer :: dimleni,dimleno ! input/output dimension length @@ -115,12 +115,12 @@ subroutine initInterp (filei, fileo, bounds) integer :: ivalue integer :: spinup_state_i, spinup_state_o integer :: decomp_cascade_state_i, decomp_cascade_state_o - integer :: npftsi, ncolsi, nlunsi, ngrcsi - integer :: npftso, ncolso, nlunso, ngrcso + integer :: npftsi, ncolsi, nlunsi, ngrcsi + integer :: npftso, ncolso, nlunso, ngrcso 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(:) @@ -148,13 +148,12 @@ subroutine initInterp (filei, fileo, bounds) call check_dim_subgrid(ncidi, ncido, dimname ='pft' , dimleni=npftsi, dimleno=npftso) 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) - + 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 write (iulog,*) 'input columns = ',ncolsi,' output columns = ',ncolso - write (iulog,*) 'input pfts = ',npftsi,' output pfts = ',npftso + write (iulog,*) 'input pfts = ',npftsi,' output pfts = ',npftso end if call check_dim_level(ncidi, ncido, dimname='levsno' ) @@ -195,28 +194,27 @@ subroutine initInterp (filei, fileo, bounds) begp_i = 1 ; endp_i = npftsi begc_i = 1 ; endc_i = ncolsi begl_i = 1 ; endl_i = nlunsi - begg_i = 1 ; endg_i = ngrcsi - + 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 - + 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(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(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)) - + allocate(grcindx(begg_o:endg_o)) ! For each output pft, find the input pft, pftindx, that is closest if (masterproc) then @@ -225,8 +223,7 @@ 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 + ! For each output column, find the input column, colindx, that is closest if (masterproc) then write(iulog,*)'finding minimum distance for columns' @@ -374,12 +371,19 @@ subroutine initInterp (filei, fileo, bounds) ! For 1D output variables !--------------------------------------------------- - else if ( ndimso == 1 ) then - - status = pio_inq_dimname(ncido, dimidso(1), vec_dimname) - if ( vec_dimname == 'pft' )then - begi = begp_i - endi = endp_i + else if ( ndimso == 1 ) then + + status = pio_inq_dimname(ncido, dimidso(1), vec_dimname) + if (trim(vec_dimname) == 'topounit') then + if (masterproc) then + write(iulog,*) 'Skipping topounit variable: ', trim(varname) + end if + CYCLE + end if + + if ( vec_dimname == 'pft' )then + begi = begp_i + endi = endp_i bego = begp_o endo = endp_o activei => pft_activei @@ -409,9 +413,9 @@ subroutine initInterp (filei, fileo, bounds) activei => grc_activei activeo => grc_activeo sgridindex => grcindx - else - call endrun(msg='ERROR interpinic: 1D variable '//trim(varname)//& - 'with unknown subgrid dimension: '//trim(vec_dimname)//& + else + call endrun(msg='ERROR interpinic: 1D variable '//trim(varname)//& + 'with unknown subgrid dimension: '//trim(vec_dimname)//& errMsg(__FILE__, __LINE__)) end if @@ -627,7 +631,7 @@ subroutine set_subgrid_dist(beg, end, dimname, ncid, active, subgrid) call ncd_io(ncid=ncid, varname='grid1d_lon' , flag='read', data=subgrid%lon , dim1name='gridunit') call ncd_io(ncid=ncid, varname='grid1d_lat' , flag='read', data=subgrid%lat , dim1name='gridunit') itemp(:) = 1 - end if + end if do n = beg,end if (itemp(n) > 0) then @@ -708,6 +712,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 @@ -863,8 +868,8 @@ logical function is_sametype (ni, no, subgridi, subgrido) end if else if (trim(subgridi%name) == 'gridcell' .and. trim(subgrido%name) == 'gridcell') then is_sametype = .true. - else - if (masterproc) then + else + if (masterproc) then write(iulog,*)'ERROR interpinic: is_sametype check on input and output type not supported' write(iulog,*)'typei = ',trim(subgridi%name) write(iulog,*)'typeo = ',trim(subgrido%name) From d7b8a34dcc08b36581800ed9eb7102f0884b7a58 Mon Sep 17 00:00:00 2001 From: jinyuntang Date: Mon, 31 Aug 2026 13:11:03 -0700 Subject: [PATCH 2/3] more restrict criterion for skip topounit variables in in-line interpolation Now for variables with topounit in the first dimension,interpinic_flag='skip' is enforced if they are reset for the beginning of a hybrid or startup run. --- .../elm/src/data_types/TopounitDataType.F90 | 4 ++-- components/elm/src/main/accumulMod.F90 | 18 ++++++++++++------ 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/components/elm/src/data_types/TopounitDataType.F90 b/components/elm/src/data_types/TopounitDataType.F90 index 7683b62f376d..9cbefa70868d 100644 --- a/components/elm/src/data_types/TopounitDataType.F90 +++ b/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 diff --git a/components/elm/src/main/accumulMod.F90 b/components/elm/src/main/accumulMod.F90 index a4f315435735..4d99e5306dd4 100644 --- a/components/elm/src/main/accumulMod.F90 +++ b/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 From d925f89090e94058708b45da45abff12d0acca15 Mon Sep 17 00:00:00 2001 From: jinyuntang Date: Mon, 31 Aug 2026 13:19:48 -0700 Subject: [PATCH 3/3] trigger endrun when topounit is requested in initInterp.F90 Instead of doing skip, now if variables are required to interpolate along topounit, the model will crash, assuming interpinic_flag is properly set when calling restartvar. This means if an interpolation is indeed needed, then the code needs proper revision. --- components/elm/src/main/initInterp.F90 | 67 ++++++++++++++------------ 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/components/elm/src/main/initInterp.F90 b/components/elm/src/main/initInterp.F90 index a62c6aba88fe..6ba470cf58b3 100644 --- a/components/elm/src/main/initInterp.F90 +++ b/components/elm/src/main/initInterp.F90 @@ -90,8 +90,8 @@ subroutine initInterp (filei, fileo, bounds) integer :: begc_o, endc_o ! output file column bounds integer :: begl_i, endl_i ! input file landunit bounds integer :: begl_o, endl_o ! output file landunit bounds - integer :: begg_i, endg_i ! input file gridcell bounds - integer :: begg_o, endg_o ! output file gridcell bounds + integer :: begg_i, endg_i ! input file gridcell bounds + integer :: begg_o, endg_o ! output file gridcell bounds integer :: nlevi,nlevo ! input/output number of levels type(file_desc_t) :: ncidi, ncido ! input/output pio fileids integer :: dimleni,dimleno ! input/output dimension length @@ -115,12 +115,12 @@ subroutine initInterp (filei, fileo, bounds) integer :: ivalue integer :: spinup_state_i, spinup_state_o integer :: decomp_cascade_state_i, decomp_cascade_state_o - integer :: npftsi, ncolsi, nlunsi, ngrcsi - integer :: npftso, ncolso, nlunso, ngrcso + integer :: npftsi, ncolsi, nlunsi, ngrcsi + integer :: npftso, ncolso, nlunso, ngrcso 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(:) @@ -148,12 +148,12 @@ subroutine initInterp (filei, fileo, bounds) call check_dim_subgrid(ncidi, ncido, dimname ='pft' , dimleni=npftsi, dimleno=npftso) 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) + 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 write (iulog,*) 'input columns = ',ncolsi,' output columns = ',ncolso - write (iulog,*) 'input pfts = ',npftsi,' output pfts = ',npftso + write (iulog,*) 'input pfts = ',npftsi,' output pfts = ',npftso end if call check_dim_level(ncidi, ncido, dimname='levsno' ) @@ -194,27 +194,27 @@ subroutine initInterp (filei, fileo, bounds) begp_i = 1 ; endp_i = npftsi begc_i = 1 ; endc_i = ncolsi begl_i = 1 ; endl_i = nlunsi - begg_i = 1 ; endg_i = ngrcsi + 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 + 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(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(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)) + allocate(grcindx(begg_o:endg_o)) ! For each output pft, find the input pft, pftindx, that is closest if (masterproc) then @@ -223,7 +223,7 @@ 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 + ! For each output column, find the input column, colindx, that is closest if (masterproc) then write(iulog,*)'finding minimum distance for columns' @@ -371,19 +371,22 @@ subroutine initInterp (filei, fileo, bounds) ! For 1D output variables !--------------------------------------------------- - else if ( ndimso == 1 ) then - - status = pio_inq_dimname(ncido, dimidso(1), vec_dimname) - if (trim(vec_dimname) == 'topounit') then - if (masterproc) then - write(iulog,*) 'Skipping topounit variable: ', trim(varname) - end if - CYCLE - end if - - if ( vec_dimname == 'pft' )then - begi = begp_i - endi = endp_i + 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 bego = begp_o endo = endp_o activei => pft_activei @@ -413,9 +416,9 @@ subroutine initInterp (filei, fileo, bounds) activei => grc_activei activeo => grc_activeo sgridindex => grcindx - else - call endrun(msg='ERROR interpinic: 1D variable '//trim(varname)//& - 'with unknown subgrid dimension: '//trim(vec_dimname)//& + else + call endrun(msg='ERROR interpinic: 1D variable '//trim(varname)//& + 'with unknown subgrid dimension: '//trim(vec_dimname)//& errMsg(__FILE__, __LINE__)) end if @@ -631,7 +634,7 @@ subroutine set_subgrid_dist(beg, end, dimname, ncid, active, subgrid) call ncd_io(ncid=ncid, varname='grid1d_lon' , flag='read', data=subgrid%lon , dim1name='gridunit') call ncd_io(ncid=ncid, varname='grid1d_lat' , flag='read', data=subgrid%lat , dim1name='gridunit') itemp(:) = 1 - end if + end if do n = beg,end if (itemp(n) > 0) then @@ -868,8 +871,8 @@ logical function is_sametype (ni, no, subgridi, subgrido) end if else if (trim(subgridi%name) == 'gridcell' .and. trim(subgrido%name) == 'gridcell') then is_sametype = .true. - else - if (masterproc) then + else + if (masterproc) then write(iulog,*)'ERROR interpinic: is_sametype check on input and output type not supported' write(iulog,*)'typei = ',trim(subgridi%name) write(iulog,*)'typeo = ',trim(subgrido%name)