From 427dcca70c2db3c4bf41980c228055b9e826bbca Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Mon, 13 Mar 2023 18:25:28 -0400 Subject: [PATCH] cicecore: correct initial condition metadata (#818) * ice_history_write: fix initial condition metadata under 'hist_avg' When writing averaged history outputs (hist_avg=.true.), this setting also affects the initial condition. Even if the actual data variables written to the initial condition are not averaged (they are taken more or less directly from the restart or the hard-coded defaults, modulo aggregation over categories), their attributes ('cell_method' and 'time_rep') imply they are averaged, and the 'bound' attribute of the 'time' variable refers to the 'time_bounds' variable. Make the metadata of the initial condition more correct by: - not writing the 'time_bounds' variable (and the corresponding 'd2' dimension) - not writing the 'bounds' attribute of the 'time' variable - not writing the 'cell_method' attributes of each variable - writing the 'time_rep' attribute of each variable as 'instantaneous' instead of 'averaged'. Do this by checking 'write_ic' at all places where we check for the value of 'hist_avg' to write the above variables and attributes in each of the 3 IO backends (binary, netcdf, pio2). * drivers/{nemo_concepts,standalone}: write initial condition at initial time In CICE_InitMod::cice_init, we call ice_calendar::advance_timestep before writing the initial condition, such that the 'time' variable in the initial condition is not zero; it has a value of 1*dt (the model time step). The initial condition filename also reflects this, since 'msec' (model seconds) also has a value of 1*dt and is used in ice_history_shared::construct_filename. This leads to the initial condition filename not corresponding to the model initialization date/time but rather 1*dt later. Since we call 'accum_hist' after initializing the forcing, any forcing field written to the initial condition has values corresponding to msec=dt, whereas the ice state corresponds to msec=0, leading to an inconsistency. Fix that by calling 'accum_hist' to write the initial condition _before_ calling 'advance_timestep'. Since we now call 'accum_hist' before initializing the forcing, any forcing field written to the initial condition have its default, hard-coded value, instead of its value at time=dt. An improvement would be to read the forcing at time=dt, write the initial condition, advance the time step, and read the forcing again, but let's not complicate things too much for now. (cherry picked from commit 9a55ad9b71444e5600c517e687b5f5d641ee6dea) Cherry-pick notes: There were some conflicts in io_{netcdf,pio2}/ice_history_write.F90 due to 26d917a (Fix QC test, fix bug in history time axis, fix history output averaging for timestep output (#624), 2021-08-19), since that commit enabled averaging over multiple time steps and thus removed the assumption that histfreq(ns) = '1' means hist_avg = .false.. I also had to add additional changes to 'ice_history_write' in the conditions that are checked before writing the NetCDF attributes since we do not yet have the 'ice_write_hist_attrs' subroutine added in 078aab4 (Merge cgridDEV branch including C grid implementation and other fixes (#715), 2022-05-10). --- .../io/io_binary/ice_history_write.F90 | 17 ++++---- .../io/io_netcdf/ice_history_write.F90 | 27 ++++++------ .../io/io_pio2/ice_history_write.F90 | 43 ++++++++++--------- .../direct/nemo_concepts/CICE_InitMod.F90 | 4 +- .../drivers/standalone/cice/CICE_InitMod.F90 | 4 +- 5 files changed, 49 insertions(+), 46 deletions(-) diff --git a/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 index b98e09814..a8d6d467f 100644 --- a/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90 @@ -158,6 +158,7 @@ subroutine ice_write_hist(ns) trim(avail_hist_fields(n)%vcomment) if (histfreq(ns) == '1' .or. .not. hist_avg & + .or. write_ic & .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & @@ -186,7 +187,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 994) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -210,7 +211,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -234,7 +235,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -258,7 +259,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -282,7 +283,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -307,7 +308,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -333,7 +334,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -359,7 +360,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else diff --git a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 index 9c6b30ee1..22cb9a40f 100644 --- a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90 @@ -155,7 +155,7 @@ subroutine ice_write_hist (ns) ! define dimensions !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_def_dim(ncid,'d2',2,boundid) if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: defining dim d2') @@ -237,7 +237,7 @@ subroutine ice_write_hist (ns) call abort_ice(subname//'ERROR: invalid calendar settings') endif - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_put_att(ncid,varid,'bounds','time_bounds') if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: time bounds') @@ -247,7 +247,7 @@ subroutine ice_write_hist (ns) ! Define attributes for time bounds if hist_avg is true !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then dimid(1) = boundid dimid(2) = timid status = nf90_def_var(ncid,'time_bounds',lprecision,dimid(1:2),varid) @@ -563,7 +563,7 @@ subroutine ice_write_hist (ns) !----------------------------------------------------------------- ! Add cell_methods attribute to variables if averaged !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then if (TRIM(avail_hist_fields(n)%vname)/='sig1' & .or.TRIM(avail_hist_fields(n)%vname)/='sig2' & .or.TRIM(avail_hist_fields(n)%vname)/='sistreave' & @@ -576,6 +576,7 @@ subroutine ice_write_hist (ns) endif if (histfreq(ns) == '1' .or. .not. hist_avg & + .or. write_ic & .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & @@ -634,13 +635,13 @@ subroutine ice_write_hist (ns) !----------------------------------------------------------------- ! Add cell_methods attribute to variables if averaged !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_put_att(ncid,varid,'cell_methods','time: mean') if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: defining cell methods for '//avail_hist_fields(n)%vname) endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = nf90_put_att(ncid,varid,'time_rep','instantaneous') else status = nf90_put_att(ncid,varid,'time_rep','averaged') @@ -875,13 +876,13 @@ subroutine ice_write_hist (ns) !----------------------------------------------------------------- ! Add cell_methods attribute to variables if averaged !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_put_att(ncid,varid,'cell_methods','time: mean') if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: defining cell methods for '//avail_hist_fields(n)%vname) endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = nf90_put_att(ncid,varid,'time_rep','instantaneous') else status = nf90_put_att(ncid,varid,'time_rep','averaged') @@ -936,13 +937,13 @@ subroutine ice_write_hist (ns) !----------------------------------------------------------------- ! Add cell_methods attribute to variables if averaged !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_put_att(ncid,varid,'cell_methods','time: mean') if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: defining cell methods for '//avail_hist_fields(n)%vname) endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = nf90_put_att(ncid,varid,'time_rep','instantaneous') else status = nf90_put_att(ncid,varid,'time_rep','averaged') @@ -997,13 +998,13 @@ subroutine ice_write_hist (ns) !----------------------------------------------------------------- ! Add cell_methods attribute to variables if averaged !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_put_att(ncid,varid,'cell_methods','time: mean') if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: defining cell methods for '//avail_hist_fields(n)%vname) endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = nf90_put_att(ncid,varid,'time_rep','instantaneous') else status = nf90_put_att(ncid,varid,'time_rep','averaged') @@ -1098,7 +1099,7 @@ subroutine ice_write_hist (ns) ! write time_bounds info !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_inq_varid(ncid,'time_bounds',varid) if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: getting time_bounds id') diff --git a/cicecore/cicedynB/infrastructure/io/io_pio2/ice_history_write.F90 b/cicecore/cicedynB/infrastructure/io/io_pio2/ice_history_write.F90 index 72a1ed97f..c51b2d5c8 100644 --- a/cicecore/cicedynB/infrastructure/io/io_pio2/ice_history_write.F90 +++ b/cicecore/cicedynB/infrastructure/io/io_pio2/ice_history_write.F90 @@ -186,7 +186,7 @@ subroutine ice_write_hist (ns) ! define dimensions !----------------------------------------------------------------- - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_def_dim(File,'d2',2,boundid) endif @@ -224,12 +224,12 @@ subroutine ice_write_hist (ns) call abort_ice(subname//'ERROR: invalid calendar settings') endif - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'bounds','time_bounds') endif ! Define attributes for time_bounds if hist_avg is true - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then dimid2(1) = boundid dimid2(2) = timid !sgl status = pio_def_var(File,'time_bounds',pio_real,dimid2,varid) @@ -473,7 +473,7 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then if (TRIM(avail_hist_fields(n)%vname)/='sig1' & .or.TRIM(avail_hist_fields(n)%vname)/='sig2' & .or.TRIM(avail_hist_fields(n)%vname)/='sistreave' & @@ -484,6 +484,7 @@ subroutine ice_write_hist (ns) endif if (histfreq(ns) == '1' .or. .not. hist_avg & + .or. write_ic & .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & @@ -527,11 +528,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -569,11 +570,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -611,11 +612,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -653,11 +654,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -695,11 +696,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -743,11 +744,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -786,11 +787,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -830,11 +831,11 @@ subroutine ice_write_hist (ns) endif ! Add cell_methods attribute to variables if averaged - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_put_att(File,varid,'cell_methods','time: mean') endif - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then status = pio_put_att(File,varid,'time_rep','instantaneous') else status = pio_put_att(File,varid,'time_rep','averaged') @@ -908,7 +909,7 @@ subroutine ice_write_hist (ns) ! write time_bounds info !----------------------------------------------------------------- - if (hist_avg .and. histfreq(ns) /= '1') then + if (hist_avg .and. histfreq(ns) /= '1' .and. .not. write_ic) then status = pio_inq_varid(File,'time_bounds',varid) time_bounds=(/time_beg(ns),time_end(ns)/) bnd_start = (/1,1/) diff --git a/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 b/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 index b2a0e3cd1..4bce4055e 100644 --- a/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 +++ b/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 @@ -185,6 +185,8 @@ subroutine cice_init if (trim(runtype) == 'continue' .or. restart) & call init_shortwave ! initialize radiative transfer + if (write_ic) call accum_hist(dt) ! write initial conditions + ! determine the time and date at the end of the first timestep call advance_timestep() @@ -215,8 +217,6 @@ subroutine cice_init call init_flux_atm ! initialize atmosphere fluxes sent to coupler call init_flux_ocn ! initialize ocean fluxes sent to coupler - if (write_ic) call accum_hist(dt) ! write initial conditions - end subroutine cice_init !======================================================================= diff --git a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 b/cicecore/drivers/standalone/cice/CICE_InitMod.F90 index 60f71fa8a..c33a604d7 100644 --- a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 +++ b/cicecore/drivers/standalone/cice/CICE_InitMod.F90 @@ -190,6 +190,8 @@ subroutine cice_init if (trim(runtype) == 'continue' .or. restart) & call init_shortwave ! initialize radiative transfer + if (write_ic) call accum_hist(dt) ! write initial conditions + ! tcraig, use advance_timestep here ! istep = istep + 1 ! update time step counters ! istep1 = istep1 + 1 @@ -222,8 +224,6 @@ subroutine cice_init call init_flux_atm ! initialize atmosphere fluxes sent to coupler call init_flux_ocn ! initialize ocean fluxes sent to coupler - if (write_ic) call accum_hist(dt) ! write initial conditions - end subroutine cice_init !=======================================================================