Fix writing of annual lnd2glc aux history file

src/drivers/mct/cime_config/namelist_definition_drv.xml

@@ -1154,12 +1154,12 @@
     </values>
   </entry>
 
-  <entry id="histaux_l2x1yr">
+  <entry id="histaux_l2x1yrg">
     <type>logical</type>
     <category>history</category>
     <group>seq_infodata_inparm</group>
     <desc>
-      turns on coupler history stream for annual sno to coupler fields.
+      turns on coupler history stream for annual lnd to coupler glc forcing fields.
       default: false
     </desc>
     <values>

src/drivers/mct/main/cime_comp_mod.F90

@@ -429,7 +429,7 @@ module cime_comp_mod
   logical :: do_hist_r2x             ! create aux files: r2x
   logical :: do_hist_l2x             ! create aux files: l2x
   logical :: do_hist_a2x24hr         ! create aux files: a2x
-  logical :: do_hist_l2x1yr          ! create aux files: l2x
+  logical :: do_hist_l2x1yrg         ! create aux files: l2x 1yr glc forcings
   logical :: do_hist_a2x             ! create aux files: a2x
   logical :: do_hist_a2x3hrp         ! create aux files: a2x 3hr precip
   logical :: do_hist_a2x3hr          ! create aux files: a2x 3hr states
@@ -914,7 +914,7 @@ subroutine cime_pre_init2()
          histaux_a2x3hrp=do_hist_a2x3hrp           , &
          histaux_a2x24hr=do_hist_a2x24hr           , &
          histaux_l2x=do_hist_l2x                   , &
-         histaux_l2x1yr=do_hist_l2x1yr             , &
+         histaux_l2x1yrg=do_hist_l2x1yrg           , &
          histaux_r2x=do_hist_r2x                   , &
          run_barriers=run_barriers                 , &
          mct_usealltoall=mct_usealltoall           , &
@@ -2094,6 +2094,8 @@ subroutine cime_run()
     logical            :: drv_pause  ! Driver writes pause restart file
     character(len=CL)  :: drv_resume ! Driver resets state from restart file
 
+    logical            :: lnd2glc_averaged_now  ! Whether lnd2glc averages were taken this timestep
+
 101 format( A, i10.8, i8, 12A, A, F8.2, A, F8.2 )
 102 format( A, i10.8, i8, A, 8L3 )
 103 format( 5A )
@@ -2207,6 +2209,8 @@ subroutine cime_run()
        if (tod            == 0) t24hr_alarm = .true.
        if (month==1 .and. day==1 .and. tod==0) t1yr_alarm = .true.
 
+       lnd2glc_averaged_now = .false.
+
        if (seq_timemgr_alarmIsOn(EClock_d,seq_timemgr_alarm_datestop)) then
           if (iamroot_CPLID) then
              write(logunit,*) ' '
@@ -2935,6 +2939,7 @@ subroutine cime_run()
                 ! NOTE - only create appropriate input to glc if the avg_alarm is on
                 if (glcrun_avg_alarm) then
                    call prep_glc_accum_avg(timer='CPL:glcprep_avg')
+                   lnd2glc_averaged_now = .true.
 
                    ! Note that l2x_gx is obtained from mapping the module variable l2gacc_lx
                    call prep_glc_calc_l2x_gx(fractions_lx, timer='CPL:glcprep_lnd2glc')
@@ -3716,15 +3721,49 @@ subroutine cime_run()
              enddo
           endif
 
-          if (do_hist_l2x1yr .and. glcrun_alarm) then
-             ! Use yr_offset=-1 so the file with fields from year 1 has time stamp
-             ! 0001-01-01 rather than 0002-01-01, etc.
-             do eli = 1,num_inst_lnd
-                suffix = component_get_suffix(lnd(eli))
-                call seq_hist_writeaux(infodata, EClock_d, lnd(eli), flow='c2x', &
-                     aname='l2x'//trim(suffix), dname='doml', &
-                     nx=lnd_nx, ny=lnd_ny, nt=1, write_now=t1yr_alarm, yr_offset=-1)
-             enddo
+          if (do_hist_l2x1yrg) then
+             ! We use a different approach here than for other aux hist files: For other
+             ! files, we let seq_hist_writeaux accumulate fields in time. However, if we
+             ! stop in the middle of an accumulation period, these accumulated fields get
+             ! reset (because they aren't written to the cpl restart file); this is
+             ! potentially a problem for this year-long accumulation. Thus, here, we use
+             ! the existing accumulated fields from prep_glc_mod, because those *do*
+             ! continue properly through a restart.
+
+             ! The logic here assumes that we average the lnd2glc fields exactly at the
+             ! year boundary - no more and no less. If that's not the case, we're likely
+             ! to be writing the wrong thing to these aux files, so we check that
+             ! assumption here.
+             if (t1yr_alarm .and. .not. lnd2glc_averaged_now) then
+                write(logunit,*) 'ERROR: histaux_l2x1yrg requested;'
+                write(logunit,*) 'it is the year boundary, but lnd2glc fields were not averaged this time step.'
+                write(logunit,*) 'One possible reason is that you are running with a stub glc model.'
+                write(logunit,*) '(It only works to request histaux_l2x1yrg if running with a prognostic glc model.)'
+                call shr_sys_abort(subname// &
+                     ' do_hist_l2x1yrg and t1yr_alarm are true, but lnd2glc_averaged_now is false')
+             end if
+             if (lnd2glc_averaged_now .and. .not. t1yr_alarm) then
+                ! If we're averaging more frequently than yearly, then just writing the
+                ! current values of the averaged fields once per year won't give the true
+                ! annual averages.
+                write(logunit,*) 'ERROR: histaux_l2x1yrg requested;'
+                write(logunit,*) 'lnd2glc fields were averaged this time step, but it is not the year boundary.'
+                write(logunit,*) '(It only works to request histaux_l2x1yrg if GLC_AVG_PERIOD is yearly.)'
+                call shr_sys_abort(subname// &
+                     ' do_hist_l2x1yrg and lnd2glc_averaged_now are true, but t1yr_alarm is false')
+             end if
+
+             if (t1yr_alarm) then
+                do eli = 1,num_inst_lnd
+                   suffix = component_get_suffix(lnd(eli))
+                   ! Use yr_offset=-1 so the file with fields from year 1 has time stamp
+                   ! 0001-01-01 rather than 0002-01-01, etc.
+                   call seq_hist_writeaux(infodata, EClock_d, lnd(eli), flow='c2x', &
+                        aname='l2x1yr_glc'//trim(suffix), dname='doml', &
+                        nx=lnd_nx, ny=lnd_ny, nt=1, write_now=.true., yr_offset=-1, &
+                        av_to_write=prep_glc_get_l2gacc_lx_one_instance(eli))
+                enddo
+             endif
           endif
 
           if (do_hist_l2x) then

src/drivers/mct/main/prep_glc_mod.F90

@@ -40,6 +40,7 @@ module prep_glc_mod
 
   public :: prep_glc_get_l2x_gx
   public :: prep_glc_get_l2gacc_lx
+  public :: prep_glc_get_l2gacc_lx_one_instance
   public :: prep_glc_get_l2gacc_lx_cnt
   public :: prep_glc_get_mapper_Sl2g
   public :: prep_glc_get_mapper_Fl2g
@@ -1086,6 +1087,12 @@ function prep_glc_get_l2gacc_lx()
     prep_glc_get_l2gacc_lx => l2gacc_lx(:)
   end function prep_glc_get_l2gacc_lx
 
+  function prep_glc_get_l2gacc_lx_one_instance(lnd_inst)
+    integer, intent(in) :: lnd_inst
+    type(mct_aVect), pointer :: prep_glc_get_l2gacc_lx_one_instance
+    prep_glc_get_l2gacc_lx_one_instance => l2gacc_lx(lnd_inst)
+  end function prep_glc_get_l2gacc_lx_one_instance
+
   function prep_glc_get_l2gacc_lx_cnt()
     integer, pointer :: prep_glc_get_l2gacc_lx_cnt
     prep_glc_get_l2gacc_lx_cnt => l2gacc_lx_cnt

src/drivers/mct/main/seq_hist_mod.F90

@@ -971,7 +971,7 @@ end subroutine seq_hist_writeavg
   !===============================================================================
 
   subroutine seq_hist_writeaux(infodata, EClock_d, comp, flow, aname, dname, &
-       nx, ny, nt, write_now, flds, yr_offset)
+       nx, ny, nt, write_now, flds, yr_offset, av_to_write)
 
     implicit none
 
@@ -989,6 +989,10 @@ subroutine seq_hist_writeaux(infodata, EClock_d, comp, flow, aname, dname, &
     character(*)             , optional, intent(in) :: flds      ! list of fields to write
     integer                  , optional, intent(in) :: yr_offset ! offset to apply to current year when generating file name
 
+    ! If av_to_write is provided, then write fields from this attribute vector.
+    ! Otherwise, get the attribute vector from 'comp', based on 'flow'.
+    type(mct_avect), target  , optional, intent(in) :: av_to_write
+
     !--- local ---
     type(mct_gGrid), pointer :: dom
     type(mct_avect), pointer :: av
@@ -1079,10 +1083,14 @@ subroutine seq_hist_writeaux(infodata, EClock_d, comp, flow, aname, dname, &
     enddo
 
     if (iamin_CPLID) then
-       if (flow == 'c2x') then
-          av => component_get_c2x_cx(comp)
-       else if (flow == 'x2c') then
-          av => component_get_x2c_cx(comp)
+       if (present(av_to_write)) then
+          av => av_to_write
+       else
+          if (flow == 'c2x') then
+             av => component_get_c2x_cx(comp)
+          else if (flow == 'x2c') then
+             av => component_get_x2c_cx(comp)
+          end if
        end if
        dom   => component_get_dom_cx(comp)
        gsmap => component_get_gsmap_cx(comp)

src/drivers/mct/shr/seq_infodata_mod.F90

@@ -153,7 +153,7 @@ MODULE seq_infodata_mod
      logical                 :: histaux_a2x3hr  ! cpl writes aux hist files: a2x 3hr states
      logical                 :: histaux_a2x3hrp ! cpl writes aux hist files: a2x 3hr precip
      logical                 :: histaux_a2x24hr ! cpl writes aux hist files: a2x daily all
-     logical                 :: histaux_l2x1yr  ! cpl writes aux hist files: l2x annual all
+     logical                 :: histaux_l2x1yrg ! cpl writes aux hist files: l2x annual glc forcings
      logical                 :: histaux_l2x     ! cpl writes aux hist files: l2x every c2l comm
      logical                 :: histaux_r2x     ! cpl writes aux hist files: r2x every c2o comm
      logical                 :: histavg_atm     ! cpl writes atm fields in average history file
@@ -389,7 +389,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid, cpl_tag)
     logical                :: histaux_a2x3hr     ! cpl writes aux hist files: a2x 3hr states
     logical                :: histaux_a2x3hrp    ! cpl writes aux hist files: a2x 2hr precip
     logical                :: histaux_a2x24hr    ! cpl writes aux hist files: a2x daily all
-    logical                :: histaux_l2x1yr     ! cpl writes aux hist files: l2x annual all
+    logical                :: histaux_l2x1yrg    ! cpl writes aux hist files: l2x annual glc forcings
     logical                :: histaux_l2x        ! cpl writes aux hist files: l2x every c2l comm
     logical                :: histaux_r2x        ! cpl writes aux hist files: r2x every c2o comm
     logical                :: histavg_atm        ! cpl writes atm fields in average history file
@@ -443,7 +443,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid, cpl_tag)
          histaux_a2x24hr,histaux_l2x   ,histaux_r2x,       &
          histavg_atm, histavg_lnd, histavg_ocn, histavg_ice, &
          histavg_rof, histavg_glc, histavg_wav, histavg_xao, &
-         histaux_l2x1yr, cpl_seq_option,                   &
+         histaux_l2x1yrg, cpl_seq_option,                   &
          eps_frac, eps_amask,                   &
          eps_agrid, eps_aarea, eps_omask, eps_ogrid,       &
          eps_oarea, esmf_map_flag,                         &
@@ -534,7 +534,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid, cpl_tag)
        histaux_a2x3hr        = .false.
        histaux_a2x3hrp       = .false.
        histaux_a2x24hr       = .false.
-       histaux_l2x1yr        = .false.
+       histaux_l2x1yrg       = .false.
        histaux_l2x           = .false.
        histaux_r2x           = .false.
        histavg_atm           = .true.
@@ -657,7 +657,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid, cpl_tag)
        infodata%histaux_a2x3hr        = histaux_a2x3hr
        infodata%histaux_a2x3hrp       = histaux_a2x3hrp
        infodata%histaux_a2x24hr       = histaux_a2x24hr
-       infodata%histaux_l2x1yr        = histaux_l2x1yr
+       infodata%histaux_l2x1yrg       = histaux_l2x1yrg
        infodata%histaux_l2x           = histaux_l2x
        infodata%histaux_r2x           = histaux_r2x
        infodata%histavg_atm           = histavg_atm
@@ -955,7 +955,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_
        budget_inst, budget_daily, budget_month, wall_time_limit,          &
        budget_ann, budget_ltann, budget_ltend , force_stop_at,            &
        histaux_a2x    , histaux_a2x1hri, histaux_a2x1hr,                  &
-       histaux_a2x3hr, histaux_a2x3hrp , histaux_l2x1yr,                  &
+       histaux_a2x3hr, histaux_a2x3hrp , histaux_l2x1yrg,                 &
        histaux_a2x24hr, histaux_l2x   , histaux_r2x     , orb_obliq,      &
        histavg_atm, histavg_lnd, histavg_ocn, histavg_ice,                &
        histavg_rof, histavg_glc, histavg_wav, histavg_xao,                &
@@ -1054,7 +1054,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_
     logical,                optional, intent(OUT) :: histaux_a2x3hr
     logical,                optional, intent(OUT) :: histaux_a2x3hrp
     logical,                optional, intent(OUT) :: histaux_a2x24hr
-    logical,                optional, intent(OUT) :: histaux_l2x1yr
+    logical,                optional, intent(OUT) :: histaux_l2x1yrg
     logical,                optional, intent(OUT) :: histaux_l2x
     logical,                optional, intent(OUT) :: histaux_r2x
     logical,                optional, intent(OUT) :: histavg_atm
@@ -1228,7 +1228,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_
     if ( present(histaux_a2x3hr) ) histaux_a2x3hr = infodata%histaux_a2x3hr
     if ( present(histaux_a2x3hrp)) histaux_a2x3hrp= infodata%histaux_a2x3hrp
     if ( present(histaux_a2x24hr)) histaux_a2x24hr= infodata%histaux_a2x24hr
-    if ( present(histaux_l2x1yr) ) histaux_l2x1yr = infodata%histaux_l2x1yr
+    if ( present(histaux_l2x1yrg)) histaux_l2x1yrg= infodata%histaux_l2x1yrg
     if ( present(histaux_l2x)    ) histaux_l2x    = infodata%histaux_l2x
     if ( present(histaux_r2x)    ) histaux_r2x    = infodata%histaux_r2x
     if ( present(histavg_atm)    ) histavg_atm    = infodata%histavg_atm
@@ -1527,7 +1527,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_
        budget_inst, budget_daily, budget_month, force_stop_at,            &
        budget_ann, budget_ltann, budget_ltend ,                           &
        histaux_a2x    , histaux_a2x1hri, histaux_a2x1hr,                  &
-       histaux_a2x3hr, histaux_a2x3hrp , histaux_l2x1yr,                  &
+       histaux_a2x3hr, histaux_a2x3hrp , histaux_l2x1yrg,                 &
        histaux_a2x24hr, histaux_l2x   , histaux_r2x     , orb_obliq,      &
        histavg_atm, histavg_lnd, histavg_ocn, histavg_ice,                &
        histavg_rof, histavg_glc, histavg_wav, histavg_xao,                &
@@ -1625,7 +1625,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_
     logical,                optional, intent(IN)    :: histaux_a2x3hr
     logical,                optional, intent(IN)    :: histaux_a2x3hrp
     logical,                optional, intent(IN)    :: histaux_a2x24hr
-    logical,                optional, intent(IN)    :: histaux_l2x1yr
+    logical,                optional, intent(IN)    :: histaux_l2x1yrg
     logical,                optional, intent(IN)    :: histaux_l2x
     logical,                optional, intent(IN)    :: histaux_r2x
     logical,                optional, intent(IN)    :: histavg_atm
@@ -1798,7 +1798,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_
     if ( present(histaux_a2x3hr) ) infodata%histaux_a2x3hr = histaux_a2x3hr
     if ( present(histaux_a2x3hrp)) infodata%histaux_a2x3hrp= histaux_a2x3hrp
     if ( present(histaux_a2x24hr)) infodata%histaux_a2x24hr= histaux_a2x24hr
-    if ( present(histaux_l2x1yr) ) infodata%histaux_l2x1yr = histaux_l2x1yr
+    if ( present(histaux_l2x1yrg)) infodata%histaux_l2x1yrg= histaux_l2x1yrg
     if ( present(histaux_l2x)    ) infodata%histaux_l2x    = histaux_l2x
     if ( present(histaux_r2x)    ) infodata%histaux_r2x    = histaux_r2x
     if ( present(histavg_atm)    ) infodata%histavg_atm    = histavg_atm
@@ -2219,7 +2219,7 @@ subroutine seq_infodata_bcast(infodata,mpicom)
     call shr_mpi_bcast(infodata%histaux_a2x3hr        ,  mpicom)
     call shr_mpi_bcast(infodata%histaux_a2x3hrp       ,  mpicom)
     call shr_mpi_bcast(infodata%histaux_a2x24hr       ,  mpicom)
-    call shr_mpi_bcast(infodata%histaux_l2x1yr        ,  mpicom)
+    call shr_mpi_bcast(infodata%histaux_l2x1yrg       ,  mpicom)
     call shr_mpi_bcast(infodata%histaux_l2x           ,  mpicom)
     call shr_mpi_bcast(infodata%histaux_r2x           ,  mpicom)
     call shr_mpi_bcast(infodata%histavg_atm           ,  mpicom)
@@ -2889,7 +2889,7 @@ SUBROUTINE seq_infodata_print( infodata )
     write(logunit,F0L) subname,'histaux_a2x3hr           = ', infodata%histaux_a2x3hr
     write(logunit,F0L) subname,'histaux_a2x3hrp          = ', infodata%histaux_a2x3hrp
     write(logunit,F0L) subname,'histaux_a2x24hr          = ', infodata%histaux_a2x24hr
-    write(logunit,F0L) subname,'histaux_l2x1yr           = ', infodata%histaux_l2x1yr
+    write(logunit,F0L) subname,'histaux_l2x1yrg          = ', infodata%histaux_l2x1yrg
     write(logunit,F0L) subname,'histaux_l2x              = ', infodata%histaux_l2x
     write(logunit,F0L) subname,'histaux_r2x              = ', infodata%histaux_r2x
     write(logunit,F0L) subname,'histavg_atm              = ', infodata%histavg_atm