Stop running 0th time step

python/ctsm/crop_calendars/generate_gdds_functions.py

@@ -60,6 +60,7 @@ def error(logger, string):
     """
     Simultaneously print ERROR messages to console and to log file
     """
+    print(string)
     logger.error(string)
     raise RuntimeError(string)
 
@@ -271,7 +272,7 @@ def import_and_process_1yr(
     else:
         chunks = None
 
-    # Get h2 file (list)
+    # Get h1 file (list)
     h1_pattern = os.path.join(indir, "*h1.*.nc")
     h1_filelist = glob.glob(h1_pattern)
     if not h1_filelist:
@@ -550,15 +551,16 @@ def import_and_process_1yr(
 
     log(logger, "   Importing accumulated GDDs...")
     clm_gdd_var = "GDDACCUM"
-    my_vars = [clm_gdd_var, "GDDHARV"]
-    pattern = os.path.join(indir, f"*h2.{this_year-1}-01-01*.nc")
-    h2_files = glob.glob(pattern)
-    if not h2_files:
-        pattern = os.path.join(indir, f"*h2.{this_year-1}-01-01*.nc.base")
+    myVars = [clm_gdd_var, "GDDHARV"]
+    patterns = [f"*h2.{this_year-1}-01*.nc", f"*h2.{this_year-1}-01*.nc.base"]
+    for p in patterns:
+        pattern = os.path.join(indir, p)
         h2_files = glob.glob(pattern)
-        if not h2_files:
-            error(logger, f"No files found matching pattern '*h2.{this_year-1}-01-01*.nc(.base)'")
-    h2_ds = import_ds(
+        if h2_files:
+            break
+    if not h2_files:
+        error(logger, f"No files found matching patterns: {patterns}")
+    h2_ds = utils.import_ds(
         h2_files,
         my_vars=my_vars,
         my_vegtypes=crops_to_read,

src/biogeochem/CNPhenologyMod.F90

@@ -392,19 +392,7 @@ subroutine CNPhenology (bounds, num_soilc, filter_soilc, num_soilp, &
             soilstate_inst, temperature_inst, atm2lnd_inst, wateratm2lndbulk_inst, cnveg_state_inst, &
             cnveg_carbonstate_inst, cnveg_nitrogenstate_inst, cnveg_carbonflux_inst, cnveg_nitrogenflux_inst)
 
-       ! BACKWARDS_COMPATIBILITY(wjs, 2022-02-03) Old restart files generated at the end
-       ! of the year can indicate that a crop was panted on Jan 1, because that used to be
-       ! the time given to the last time step of the year. This would cause problems if we
-       ! ran CropPhenology in time step 0, because now time step 0 is labeled as Dec 31,
-       ! so CropPhenology would see the crop as having been planted 364 days ago, and so
-       ! would want to harvest this newly-planted crop. To avoid this situation, we avoid
-       ! calling CropPhenology on time step 0.
-       !
-       ! This .not. is_first_step() condition can be removed either when we can rely on
-       ! all restart files having been generated with
-       ! https://github.com/ESCOMP/CTSM/issues/1623 resolved, or we stop having a time
-       ! step 0 (https://github.com/ESCOMP/CTSM/issues/925).
-       if (num_pcropp > 0 .and. .not. is_first_step()) then
+       if (num_pcropp > 0) then
           call CropPhenology(num_pcropp, filter_pcropp, &
                waterdiagnosticbulk_inst, temperature_inst, crop_inst, canopystate_inst, cnveg_state_inst, &
                cnveg_carbonstate_inst, cnveg_nitrogenstate_inst, cnveg_carbonflux_inst, cnveg_nitrogenflux_inst, &

src/biogeochem/CNVegetationFacade.F90

@@ -1285,7 +1285,7 @@ subroutine EndOfTimeStepVegDynamics(this, bounds, num_natvegp, filter_natvegp, &
        ! Call dv (dynamic vegetation) at last time step of year
 
        call t_startf('d2dgvm')
-       if (is_end_curr_year() .and. .not. is_first_step())  then
+       if (is_end_curr_year())  then
 
           ! Get date info.  kyr is used in lpj().  At end of first year, kyr = 2.
           call get_curr_date(yr, mon, day, sec)
@@ -1337,7 +1337,7 @@ subroutine WriteHistory(this, bounds)
 
     ! Write to CNDV history buffer if appropriate
     if (use_cndv) then
-       if (is_end_curr_year() .and. .not. is_first_step())  then
+       if (is_end_curr_year())  then
           call t_startf('clm_drv_io_hdgvm')
           call CNDVHist( bounds, this%dgvs_inst )
           if (masterproc) write(iulog,*) 'Annual CNDV calculations are complete'

src/biogeochem/CropType.F90

@@ -930,7 +930,7 @@ subroutine CropIncrementYear (this, num_pcropp, filter_pcropp)
     ! Update nyrs when it's the end of the year (unless it's the very start of the
     ! run). This assumes that, if this patch is active at the end of the year, then it was
     ! active for the whole year.
-    if ((kmo == 1 .and. kda == 1 .and. mcsec == 0) .and. .not. is_first_step()) then
+    if ((kmo == 1 .and. kda == 1 .and. mcsec == 0)) then
        do fp = 1, num_pcropp
           p = filter_pcropp(fp)
 

src/biogeochem/dynCNDVMod.F90

@@ -99,7 +99,7 @@ subroutine dynCNDV_interp( bounds, dgvs_inst)
           patch%wtcol(p)   = dgvs_inst%fpcgrid_patch(p) + &
                     wt1 * (dgvs_inst%fpcgridold_patch(p) - dgvs_inst%fpcgrid_patch(p))
 
-          if (mon==1 .and. day==1 .and. sec==dtime .and. nstep>0) then
+          if (mon==1 .and. day==1 .and. sec==dtime) then
              dgvs_inst%fpcgridold_patch(p) = dgvs_inst%fpcgrid_patch(p)
           end if
        end if

src/cpl/lilac/lnd_comp_esmf.F90

@@ -559,7 +559,6 @@ subroutine lnd_run(gcomp, import_state, export_state, clock, rc)
     integer                :: nstep          ! time step index
     logical                :: rstwr          ! .true. ==> write restart file before returning
     logical                :: nlend          ! .true. ==> last time-step
-    logical                :: dosend         ! true => send data back to driver
     logical                :: doalb          ! .true. ==> do albedo calculation on this time step
     real(r8)               :: nextsw_cday    ! calday from clock of next radiation computation
     real(r8)               :: caldayp1       ! ctsm calday plus dtime offset
@@ -626,161 +625,146 @@ subroutine lnd_run(gcomp, import_state, export_state, clock, rc)
     !--------------------------------
 
     dtime = get_step_size()
-    dosend = .false.
-    do while(.not. dosend)
-
-       ! We assume that the land model time step matches the coupling interval. However,
-       ! we still need this while loop to handle the initial time step (time 0). We may
-       ! want to get rid of this time step 0 in the lilac coupling, at which point we
-       ! should be able to remove this while loop and dosend variable.
-       !
-       ! See also https://github.com/ESCOMP/CTSM/issues/925
-       nstep = get_nstep()
-       if (nstep > 0) then
-          dosend = .true.
-       end if
 
-       !--------------------------------
-       ! Determine calendar day info
-       !--------------------------------
+    ! We assume that the land model time step matches the coupling interval.
+    nstep = get_nstep()
 
-       calday = get_curr_calday(reuse_day_365_for_day_366=.true.)
-       caldayp1 = get_curr_calday(offset=dtime, reuse_day_365_for_day_366=.true.)
+    !--------------------------------
+    ! Determine calendar day info
+    !--------------------------------
 
-       !--------------------------------
-       ! Get time of next atmospheric shortwave calculation
-       !--------------------------------
+    calday = get_curr_calday(reuse_day_365_for_day_366=.true.)
+    caldayp1 = get_curr_calday(offset=dtime, reuse_day_365_for_day_366=.true.)
 
-       ! TODO(NS): nextsw_cday should come directly from atmosphere!
-       ! For now I am setting nextsw_cday to be the same caldayp1
-       !
-       ! See also https://github.com/ESCOMP/CTSM/issues/860
+    !--------------------------------
+    ! Get time of next atmospheric shortwave calculation
+    !--------------------------------
 
-       nextsw_cday = calday
-       if (masterproc) then
-          write(iulog,*) trim(subname) // '... nextsw_cday is : ', nextsw_cday
-       end if
+    ! TODO(NS): nextsw_cday should come directly from atmosphere!
+    ! For now I am setting nextsw_cday to be the same caldayp1
+    !
+    ! See also https://github.com/ESCOMP/CTSM/issues/860
 
-       !--------------------------------
-       ! Obtain orbital values
-       !--------------------------------
-
-       call shr_orb_decl( calday     , eccen, mvelpp, lambm0, obliqr, declin  , eccf )
-       call shr_orb_decl( nextsw_cday, eccen, mvelpp, lambm0, obliqr, declinp1, eccf )
-
-       if (masterproc) then
-          write(iulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-          write(iulog,F02) 'nextsw_cday is :  ', nextsw_cday
-          write(iulog,F02) 'calday is      :  ', calday
-          write(iulog,F02) 'eccen is       :  ', eccen
-          write(iulog,F02) 'mvelpp is      :  ', mvelpp
-          write(iulog,F02) 'lambm0 is      :  ', lambm0
-          write(iulog,F02) 'obliqr is      :  ', obliqr
-          write(iulog,F02) 'declin is      :  ', declin
-          write(iulog,F02) 'declinp1 is    :  ', declinp1
-          write(iulog,*  ) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-       end if
+    nextsw_cday = calday
+    if (masterproc) then
+       write(iulog,*) trim(subname) // '... nextsw_cday is : ', nextsw_cday
+    end if
 
-       !--------------------------------
-       ! Determine doalb based on nextsw_cday sent from atm model
-       !--------------------------------
-
-       if (nstep == 0) then
-          doalb = .false.
-          nextsw_cday = caldayp1
-       else if (nstep == 1) then
-          !doalb = (abs(nextsw_cday- caldayp1) < 1.e-10_r8)
-          doalb = .false.
-       else
-          doalb = (nextsw_cday >= -0.5_r8)
-       end if
+    !--------------------------------
+    ! Obtain orbital values
+    !--------------------------------
 
-       if (masterproc) then
-          write(iulog,*) '------------  LILAC  ----------------'
-          write(iulog,*) 'nstep       : ', nstep
-          write(iulog,*) 'calday      : ', calday
-          write(iulog,*) 'caldayp1    : ', caldayp1
-          write(iulog,*) 'nextsw_cday : ', nextsw_cday
-          write(iulog,*) 'doalb       : ', doalb
-          write(iulog,*) '-------------------------------------'
-       end if
+    call shr_orb_decl( calday     , eccen, mvelpp, lambm0, obliqr, declin  , eccf )
+    call shr_orb_decl( nextsw_cday, eccen, mvelpp, lambm0, obliqr, declinp1, eccf )
 
-       call update_rad_dtime(doalb)
+    if (masterproc) then
+       write(iulog,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+       write(iulog,F02) 'nextsw_cday is :  ', nextsw_cday
+       write(iulog,F02) 'calday is      :  ', calday
+       write(iulog,F02) 'eccen is       :  ', eccen
+       write(iulog,F02) 'mvelpp is      :  ', mvelpp
+       write(iulog,F02) 'lambm0 is      :  ', lambm0
+       write(iulog,F02) 'obliqr is      :  ', obliqr
+       write(iulog,F02) 'declin is      :  ', declin
+       write(iulog,F02) 'declinp1 is    :  ', declinp1
+       write(iulog,*  ) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+    end if
 
-       !--------------------------------
-       ! Determine if time to write restart
-       !--------------------------------
+    !--------------------------------
+    ! Determine doalb based on nextsw_cday sent from atm model
+    !--------------------------------
 
-       call ESMF_ClockGetAlarm(clock, alarmname='lilac_restart_alarm', alarm=alarm, rc=rc)
-       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    if (nstep == 1) then
+       !doalb = (abs(nextsw_cday- caldayp1) < 1.e-10_r8)
+       doalb = .false.
+    else
+       doalb = (nextsw_cday >= -0.5_r8)
+    end if
 
-       if (ESMF_AlarmIsRinging(alarm, rc=rc)) then
-          if (ChkErr(rc,__LINE__,u_FILE_u)) return
-          rstwr = .true.
-          call ESMF_AlarmRingerOff( alarm, rc=rc )
-          if (ChkErr(rc,__LINE__,u_FILE_u)) return
-       else
-          rstwr = .false.
-       endif
-       if (masterproc) then
-          write(iulog,*)' restart alarm is ',rstwr
-       end if
+    if (masterproc) then
+       write(iulog,*) '------------  LILAC  ----------------'
+       write(iulog,*) 'nstep       : ', nstep
+       write(iulog,*) 'calday      : ', calday
+       write(iulog,*) 'caldayp1    : ', caldayp1
+       write(iulog,*) 'nextsw_cday : ', nextsw_cday
+       write(iulog,*) 'doalb       : ', doalb
+       write(iulog,*) '-------------------------------------'
+    end if
 
-       !--------------------------------
-       ! Determine if time to stop
-       !--------------------------------
+    call update_rad_dtime(doalb)
 
-       call ESMF_ClockGetAlarm(clock, alarmname='lilac_stop_alarm', alarm=alarm, rc=rc)
-       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    !--------------------------------
+    ! Determine if time to write restart
+    !--------------------------------
 
-       if (ESMF_AlarmIsRinging(alarm, rc=rc)) then
-          if (ChkErr(rc,__LINE__,u_FILE_u)) return
-          nlend = .true.
-          call ESMF_AlarmRingerOff( alarm, rc=rc )
-          if (ChkErr(rc,__LINE__,u_FILE_u)) return
-       else
-          nlend = .false.
-       endif
-       if (masterproc) then
-          write(iulog,*)' stop alarm is ',nlend
-       end if
+    call ESMF_ClockGetAlarm(clock, alarmname='lilac_restart_alarm', alarm=alarm, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-       !--------------------------------
-       ! Run CTSM
-       !--------------------------------
+    if (ESMF_AlarmIsRinging(alarm, rc=rc)) then
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       rstwr = .true.
+       call ESMF_AlarmRingerOff( alarm, rc=rc )
+       if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    else
+       rstwr = .false.
+    endif
+    if (masterproc) then
+       write(iulog,*)' restart alarm is ',rstwr
+    end if
 
-       call t_barrierf('sync_ctsm_run1', mpicom)
+    !--------------------------------
+    ! Determine if time to stop
+    !--------------------------------
 
-       ! Restart File - use nexttimestr rather than currtimestr here since that is the time at the end of
-       ! the timestep and is preferred for restart file names
-       ! TODO: is this correct for lilac?
+    call ESMF_ClockGetAlarm(clock, alarmname='lilac_stop_alarm', alarm=alarm, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-       call ESMF_ClockGetNextTime(clock, nextTime=nextTime, rc=rc)
+    if (ESMF_AlarmIsRinging(alarm, rc=rc)) then
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
-       call ESMF_TimeGet(nexttime, yy=yr_lilac, mm=mon_lilac, dd=day_lilac, s=tod_lilac, rc=rc)
+       nlend = .true.
+       call ESMF_AlarmRingerOff( alarm, rc=rc )
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
-       write(rdate,'(i4.4,"-",i2.2,"-",i2.2,"-",i5.5)') yr_lilac, mon_lilac, day_lilac, tod_lilac
+    else
+       nlend = .false.
+    endif
+    if (masterproc) then
+       write(iulog,*)' stop alarm is ',nlend
+    end if
 
-       call t_startf ('ctsm_run')
-       call clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, rof_prognostic=.false.)
-       call t_stopf ('ctsm_run')
+    !--------------------------------
+    ! Run CTSM
+    !--------------------------------
 
-       !--------------------------------
-       ! Pack export state
-       !--------------------------------
+    call t_barrierf('sync_ctsm_run1', mpicom)
 
-       call t_startf ('lc_lnd_export')
-       call export_fields(export_state, bounds,  rc)
-       if (ChkErr(rc,__LINE__,u_FILE_u)) return
-       call t_stopf ('lc_lnd_export')
+    ! Restart File - use nexttimestr rather than currtimestr here since that is the time at the end of
+    ! the timestep and is preferred for restart file names
+    ! TODO: is this correct for lilac?
 
-       !--------------------------------
-       ! Advance ctsm time step
-       !--------------------------------
+    call ESMF_ClockGetNextTime(clock, nextTime=nextTime, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call ESMF_TimeGet(nexttime, yy=yr_lilac, mm=mon_lilac, dd=day_lilac, s=tod_lilac, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    write(rdate,'(i4.4,"-",i2.2,"-",i2.2,"-",i5.5)') yr_lilac, mon_lilac, day_lilac, tod_lilac
 
-       call advance_timestep()
+    call t_startf ('ctsm_run')
+    call clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, rof_prognostic=.false.)
+    call t_stopf ('ctsm_run')
 
-    end do
+    !--------------------------------
+    ! Pack export state
+    !--------------------------------
+
+    call t_startf ('lc_lnd_export')
+    call export_fields(export_state, bounds,  rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    call t_stopf ('lc_lnd_export')
+
+    !--------------------------------
+    ! Advance ctsm time step
+    !--------------------------------
+
+    call advance_timestep()
 
     !--------------------------------
     ! Check that internal clock is in sync with lilac driver clock