From ced0dd23b8de39643930da943739a9512b3f0aab Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Mon, 10 Oct 2022 12:30:18 -0600 Subject: [PATCH 1/4] Correct units in history output --- cicecore/cicedynB/analysis/ice_history_fsd.F90 | 8 ++++---- doc/source/science_guide/sg_fundvars.rst | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_history_fsd.F90 b/cicecore/cicedynB/analysis/ice_history_fsd.F90 index 50fee99e7..18e936e13 100644 --- a/cicecore/cicedynB/analysis/ice_history_fsd.F90 +++ b/cicecore/cicedynB/analysis/ice_history_fsd.F90 @@ -137,7 +137,7 @@ subroutine init_hist_fsd_2D do ns = 1, nstreams if (f_wave_sig_ht(1:1) /= 'x') & - call define_hist_field(n_wave_sig_ht,"wave_sig_ht","1",tstr2D, tcstr, & + call define_hist_field(n_wave_sig_ht,"wave_sig_ht","m",tstr2D, tcstr, & "significant height of wind and swell waves", & "from attenuated spectrum in ice", c1, c0, & ns, f_wave_sig_ht) @@ -147,7 +147,7 @@ subroutine init_hist_fsd_2D "for waves", c1, c0, & ns, f_aice_ww) if (f_diam_ww(1:1) /= 'x') & - call define_hist_field(n_diam_ww,"diam_ww","1",tstr2D, tcstr, & + call define_hist_field(n_diam_ww,"diam_ww","m",tstr2D, tcstr, & "Average (number) diameter of floes > Dmin", & "for waves", c1, c0, & ns, f_diam_ww) @@ -216,7 +216,7 @@ subroutine init_hist_fsd_3Df if (histfreq(ns) /= 'x') then if (f_afsd(1:1) /= 'x') & - call define_hist_field(n_afsd,"afsd", "1", tstr3Df, tcstr, & + call define_hist_field(n_afsd,"afsd", "1/m", tstr3Df, tcstr, & "areal floe size distribution", & "per unit bin width ", c1, c0, ns, f_afsd) if (f_dafsd_newi(1:1) /= 'x') & @@ -272,7 +272,7 @@ subroutine init_hist_fsd_4Df if (histfreq(ns) /= 'x') then if (f_afsdn(1:1) /= 'x') & - call define_hist_field(n_afsdn,"afsdn","1",tstr4Df, tcstr, & + call define_hist_field(n_afsdn,"afsdn","1/m",tstr4Df, tcstr, & "areal floe size and thickness distribution", & "per unit bin width", c1, c0, ns, f_afsdn) diff --git a/doc/source/science_guide/sg_fundvars.rst b/doc/source/science_guide/sg_fundvars.rst index e80b68f4a..5b5703266 100644 --- a/doc/source/science_guide/sg_fundvars.rst +++ b/doc/source/science_guide/sg_fundvars.rst @@ -14,7 +14,7 @@ modeling is to describe the evolution of the ice thickness distribution In addition to an ice thickness distribution, CICE includes an optional capability for a floe size distribution. -Ice floe horizontal size may change through vertical and lateral growth and melting of existing floes, freezing of new ice, wave breaking, and welding of floes in freezing conditions. The floe size distribution (FSD) is a probability function that characterizes this variability. The scheme is based on the theoretical framework described in :cite:`Horvat15` for a joint floe size and thickness distribution (FSTD), and was implemented by :cite:`Roach18` and :cite:`Roach19`. The joint floe size distribution is carried as an area-weighted tracer, defined as the fraction of ice belonging to a given thickness category with lateral floe size belong to a given floe size class. This development includes interactions between sea ice and ocean surface waves. Input data on ocean surface wave spectra at a single time is provided for testing, but as with the other CICE datasets, it should not be used for production runs or publications. +Ice floe horizontal size may change through vertical and lateral growth and melting of existing floes, freezing of new ice, wave breaking, and welding of floes in freezing conditions. The floe size distribution (FSD) is a probability function that characterizes this variability. The scheme is based on the theoretical framework described in :cite:`Horvat15` for a joint floe size and thickness distribution (FSTD), and was implemented by :cite:`Roach18` and :cite:`Roach19`. The joint floe size distribution is carried as an area-weighted tracer, defined as the fraction of ice belonging to a given thickness category with lateral floe size belong to a given floe size class. This development includes interactions between sea ice and ocean surface waves. Input data on ocean surface wave spectra at a single time is provided for testing, but as with the other CICE datasets, it should not be used for production runs or publications. It is not recommended to use the FSD without ocean surface waves. Additional information about the ITD and joint FSTD for CICE can be found in the `Icepack documentation `_. From 83d5256ff8ee8a8d0e7f15e945dfb16521a2897f Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Fri, 14 Oct 2022 13:07:41 -0600 Subject: [PATCH 2/4] Add warning for no waves with FSD, make FSD tendencies per s --- cicecore/cicedynB/analysis/ice_history.F90 | 2 +- .../cicedynB/analysis/ice_history_fsd.F90 | 25 +++++++++++-------- cicecore/cicedynB/general/ice_init.F90 | 11 +++++--- 3 files changed, 23 insertions(+), 15 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_history.F90 b/cicecore/cicedynB/analysis/ice_history.F90 index caaa56295..f5e7d0d16 100644 --- a/cicecore/cicedynB/analysis/ice_history.F90 +++ b/cicecore/cicedynB/analysis/ice_history.F90 @@ -3546,7 +3546,7 @@ subroutine accum_hist (dt) call accum_hist_drag (iblk) ! floe size distribution - call accum_hist_fsd (iblk) + call accum_hist_fsd (dt, iblk) ! advanced snow physics call accum_hist_snow (iblk) diff --git a/cicecore/cicedynB/analysis/ice_history_fsd.F90 b/cicecore/cicedynB/analysis/ice_history_fsd.F90 index 18e936e13..bf7fa40b4 100644 --- a/cicecore/cicedynB/analysis/ice_history_fsd.F90 +++ b/cicecore/cicedynB/analysis/ice_history_fsd.F90 @@ -220,23 +220,23 @@ subroutine init_hist_fsd_3Df "areal floe size distribution", & "per unit bin width ", c1, c0, ns, f_afsd) if (f_dafsd_newi(1:1) /= 'x') & - call define_hist_field(n_dafsd_newi,"dafsd_newi","1",tstr3Df, tcstr, & + call define_hist_field(n_dafsd_newi,"dafsd_newi","1/s",tstr3Df, tcstr, & "Change in fsd: new ice", & "Avg over freq period", c1, c0, ns, f_dafsd_newi) if (f_dafsd_latg(1:1) /= 'x') & - call define_hist_field(n_dafsd_latg,"dafsd_latg","1",tstr3Df, tcstr, & + call define_hist_field(n_dafsd_latg,"dafsd_latg","1/s",tstr3Df, tcstr, & "Change in fsd: lateral growth", & "Avg over freq period", c1, c0, ns, f_dafsd_latg) if (f_dafsd_latm(1:1) /= 'x') & - call define_hist_field(n_dafsd_latm,"dafsd_latm","1",tstr3Df, tcstr, & + call define_hist_field(n_dafsd_latm,"dafsd_latm","1/s",tstr3Df, tcstr, & "Change in fsd: lateral melt", & "Avg over freq period", c1, c0, ns, f_dafsd_latm) if (f_dafsd_wave(1:1) /= 'x') & - call define_hist_field(n_dafsd_wave,"dafsd_wave","1",tstr3Df, tcstr, & + call define_hist_field(n_dafsd_wave,"dafsd_wave","1/s",tstr3Df, tcstr, & "Change in fsd: waves", & "Avg over freq period", c1, c0, ns, f_dafsd_wave) if (f_dafsd_weld(1:1) /= 'x') & - call define_hist_field(n_dafsd_weld,"dafsd_weld","1",tstr3Df, tcstr, & + call define_hist_field(n_dafsd_weld,"dafsd_weld","1/s",tstr3Df, tcstr, & "Change in fsd: welding", & "Avg over freq period", c1, c0, ns, f_dafsd_weld) endif ! if (histfreq(ns) /= 'x') @@ -288,7 +288,7 @@ end subroutine init_hist_fsd_4Df ! accumulate average ice quantities or snapshots ! author: Elizabeth C. Hunke, LANL - subroutine accum_hist_fsd (iblk) + subroutine accum_hist_fsd (dt, iblk) use ice_blocks, only: nx_block, ny_block use ice_constants, only: c0, c1, c2, c4 @@ -298,6 +298,9 @@ subroutine accum_hist_fsd (iblk) use ice_arrays_column, only: wave_sig_ht, floe_rad_c, floe_binwidth, & d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld + real (kind=dbl_kind), intent(in) :: & + dt ! time step + integer (kind=int_kind), intent(in) :: & iblk ! block index @@ -452,19 +455,19 @@ subroutine accum_hist_fsd (iblk) if (f_dafsd_newi(1:1)/= 'x') & call accum_hist_field(n_dafsd_newi-n3Dacum, iblk, nfsd_hist, & - d_afsd_newi(:,:,1:nfsd_hist,iblk), a3Df) + d_afsd_newi(:,:,1:nfsd_hist,iblk)/dt, a3Df) if (f_dafsd_latg(1:1)/= 'x') & call accum_hist_field(n_dafsd_latg-n3Dacum, iblk, nfsd_hist, & - d_afsd_latg(:,:,1:nfsd_hist,iblk), a3Df) + d_afsd_latg(:,:,1:nfsd_hist,iblk)/dt, a3Df) if (f_dafsd_latm(1:1)/= 'x') & call accum_hist_field(n_dafsd_latm-n3Dacum, iblk, nfsd_hist, & - d_afsd_latm(:,:,1:nfsd_hist,iblk), a3Df) + d_afsd_latm(:,:,1:nfsd_hist,iblk)/dt, a3Df) if (f_dafsd_wave(1:1)/= 'x') & call accum_hist_field(n_dafsd_wave-n3Dacum, iblk, nfsd_hist, & - d_afsd_wave(:,:,1:nfsd_hist,iblk), a3Df) + d_afsd_wave(:,:,1:nfsd_hist,iblk)/dt, a3Df) if (f_dafsd_weld(1:1)/= 'x') & call accum_hist_field(n_dafsd_weld-n3Dacum, iblk, nfsd_hist, & - d_afsd_weld(:,:,1:nfsd_hist,iblk), a3Df) + d_afsd_weld(:,:,1:nfsd_hist,iblk)/dt, a3Df) endif ! a3Df allocated ! 4D floe size, thickness category fields diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index e0ebdfbed..39f250f98 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -1582,7 +1582,12 @@ subroutine input_data wave_spec = .false. if (tr_fsd .and. (trim(wave_spec_type) /= 'none')) wave_spec = .true. - + if (tr_fsd .and. (trim(wave_spec_type) == 'none')) then + if (my_task == master_task) then + write(nu_diag,*) subname//' WARNING: tr_fsd=T but wave_spec=F - not recommended' + endif + end if + ! compute grid locations for thermo, u and v fields grid_ice_thrm = 'T' @@ -2074,7 +2079,7 @@ subroutine input_data if (wave_spec) then tmpstr2 = ' : use wave spectrum for floe size distribution' else - tmpstr2 = ' : floe size distribution does not use wave spectrum' + tmpstr2 = 'WARNING : floe size distribution does not use wave spectrum' endif write(nu_diag,1010) ' wave_spec = ', wave_spec,trim(tmpstr2) if (wave_spec) then @@ -2085,7 +2090,7 @@ subroutine input_data elseif (trim(wave_spec_type) == 'constant') then tmpstr2 = ' : constant wave spectrum data file provided for testing' elseif (trim(wave_spec_type) == 'random') then - tmpstr2 = ' : wave data file provided, spectrum generated using random number' + tmpstr2 = ' : wave data file provided, sea surface height generated using random number' else tmpstr2 = ' : unknown value' endif From dcb8792c8731f4a21b9534132b92b13cf3b523df Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Fri, 14 Oct 2022 14:54:28 -0600 Subject: [PATCH 3/4] Add capability for time-varying wave forcing --- cicecore/cicedynB/general/ice_forcing.F90 | 191 +++++++++++++++++++++- cicecore/cicedynB/general/ice_init.F90 | 6 +- 2 files changed, 187 insertions(+), 10 deletions(-) diff --git a/cicecore/cicedynB/general/ice_forcing.F90 b/cicecore/cicedynB/general/ice_forcing.F90 index 381686c9b..ff79778c5 100755 --- a/cicecore/cicedynB/general/ice_forcing.F90 +++ b/cicecore/cicedynB/general/ice_forcing.F90 @@ -22,7 +22,7 @@ module ice_forcing use ice_boundary, only: ice_HaloUpdate use ice_blocks, only: nx_block, ny_block use ice_domain, only: halo_info - use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global + use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global, nfreq use ice_communicate, only: my_task, master_task use ice_calendar, only: istep, istep1, & msec, mday, mmonth, myear, yday, daycal, & @@ -116,6 +116,9 @@ module ice_forcing topmelt_data, & botmelt_data + real (kind=dbl_kind), dimension(:,:,:,:,:), allocatable :: & + wave_spectrum_data ! field values at 2 temporal data points + character(char_len), public :: & atm_data_format, & ! 'bin'=binary or 'nc'=netcdf ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf @@ -229,6 +232,7 @@ subroutine alloc_forcing ocn_frc_m(nx_block,ny_block, max_blocks,nfld,12), & ! ocn data for 12 months topmelt_file(ncat), & botmelt_file(ncat), & + wave_spectrum_data(nx_block,ny_block,nfreq,2,max_blocks), & stat=ierr) if (ierr/=0) call abort_ice('(alloc_forcing): Out of Memory') @@ -1527,6 +1531,44 @@ subroutine interpolate_data (field_data, field) end subroutine interpolate_data +!======================================================================= + + subroutine interpolate_wavespec_data (field_data, field) + +! Linear interpolation + +! author: Elizabeth C. Hunke, LANL + + use ice_domain, only: nblocks + + real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,2,max_blocks), intent(in) :: & + field_data ! 2 values used for interpolation + + real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,max_blocks), intent(out) :: & + field ! interpolated field + + ! local variables + + integer (kind=int_kind) :: i,j, iblk, freq + + character(len=*), parameter :: subname = '(interpolate data)' + + !$OMP PARALLEL DO PRIVATE(iblk,i,j) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + do freq = 1, nfreq + field(i,j,freq,iblk) = c1intp * field_data(i,j,freq,1,iblk) & + + c2intp * field_data(i,j,freq,2,iblk) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + + end subroutine interpolate_wavespec_data + + !======================================================================= subroutine file_year (data_file, yr) @@ -5566,7 +5608,6 @@ subroutine get_wave_spec use ice_arrays_column, only: wave_spectrum, & dwavefreq, wavefreq use ice_constants, only: c0 - use ice_domain_size, only: nfreq use ice_timers, only: ice_timer_start, ice_timer_stop, timer_fsd ! local variables @@ -5592,7 +5633,6 @@ subroutine get_wave_spec ! if no wave data is provided, wave_spectrum is zero everywhere wave_spectrum(:,:,:,:) = c0 - wave_spec_dir = ocn_data_dir debug_forcing = .false. ! wave spectrum and frequencies @@ -5610,10 +5650,7 @@ subroutine get_wave_spec file=__FILE__, line=__LINE__) else #ifdef USE_NETCDF - call ice_open_nc(wave_spec_file,fid) - call ice_read_nc_xyf (fid, 1, 'efreq', wave_spectrum(:,:,:,:), debug_forcing, & - field_loc_center, field_type_scalar) - call ice_close_nc(fid) + call wave_spec_data #else write (nu_diag,*) "wave spectrum file not available, requires cpp USE_NETCDF" write (nu_diag,*) "wave spectrum file not available, using default profile" @@ -5628,6 +5665,146 @@ subroutine get_wave_spec end subroutine get_wave_spec +!======================================================================= +! +! Read in wave spectrum forcing as a function of time. 6 hourly +! LR started working from JRA55_data routine +! Changed fields, and changed 3 hourly to 6 hourly +! + subroutine wave_spec_data + + use ice_blocks, only: block, get_block + use ice_global_reductions, only: global_minval, global_maxval + use ice_domain, only: nblocks, distrb_info, blocks_ice + use ice_arrays_column, only: wave_spectrum, & + dwavefreq, wavefreq + use ice_read_write, only: ice_read_nc_xyf + use ice_grid, only: hm, tlon, tlat, tmask, umask + use ice_calendar, only: days_per_year, use_leap_years + + integer (kind=int_kind) :: & + ncid , & ! netcdf file id + i, j, freq , & + ixm,ixx,ixp , & ! record numbers for neighboring months + recnum , & ! record number + maxrec , & ! maximum record number + recslot , & ! spline slot for current record + midmonth , & ! middle day of month + dataloc , & ! = 1 for data located in middle of time interval + ! = 2 for date located at end of time interval + iblk , & ! block index + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + yr ! current forcing year + + real (kind=dbl_kind) :: & + sec6hr , & ! number of seconds in 3 hours + secday , & ! number of seconds in day + vmin, vmax + + logical (kind=log_kind) :: readm, read6,debug_n_d + + type (block) :: & + this_block ! block information for current block + + real(kind=dbl_kind), dimension(nfreq) :: & + wave_spectrum_profile ! wave spectrum + + character(len=64) :: fieldname !netcdf field name + character(char_len_long) :: spec_file + character(char_len) :: wave_spec_type + logical (kind=log_kind) :: wave_spec + character(len=*), parameter :: subname = '(wave_spec_data)' + + + + debug_n_d = .false. !usually false + + call icepack_query_parameters(secday_out=secday) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + call icepack_init_wave(nfreq, & + wave_spectrum_profile, & + wavefreq, dwavefreq) + + + !spec_file = trim(ocn_data_dir)//'/'//trim(wave_spec_file) + spec_file = trim(wave_spec_file) + wave_spectrum_data = c0 + wave_spectrum = c0 + yr = fyear ! current year + !------------------------------------------------------------------- + ! 6-hourly data + ! + ! Assume that the 6-hourly value is located at the end of the + ! 6-hour period. This is the convention for NCEP reanalysis data. + ! E.g. record 1 gives conditions at 6 am GMT on 1 January. + !------------------------------------------------------------------- + + dataloc = 2 ! data located at end of interval + sec6hr = secday/c4 ! seconds in 6 hours + !maxrec = 2920 ! 365*8; for leap years = 366*8 + + if (use_leap_years) days_per_year = 366 !overrides setting of 365 in ice_calendar + maxrec = days_per_year*4 + + if(days_per_year == 365 .and. (mod(yr, 4) == 0)) then + call abort_ice('days_per_year should be set to 366 for leap years') + end if + + ! current record number + recnum = 4*int(yday) - 3 + int(real(msec,kind=dbl_kind)/sec6hr) + + ! Compute record numbers for surrounding data (2 on each side) + + ixm = mod(recnum+maxrec-2,maxrec) + 1 + ixx = mod(recnum-1, maxrec) + 1 + + ! Compute interpolation coefficients + ! If data is located at the end of the time interval, then the + ! data value for the current record goes in slot 2 + + recslot = 2 + ixp = -99 + call interp_coeff (recnum, recslot, sec6hr, dataloc) + + ! Read + read6 = .false. + if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true. + !------------------------------------------------------------------- + ! File is NETCDF + ! file variable names are: + ! efreq (wave spectrum, energy as a function of wave frequency UNITS) + !------------------------------------------------------------------- + call ice_open_nc(spec_file,ncid) + + call ice_read_nc_xyf(ncid,recnum,'efreq',wave_spectrum_data(:,:,:,1,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_read_nc_xyf(ncid,recnum,'efreq',wave_spectrum_data(:,:,:,2,:),debug_n_d, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + call ice_close_nc(ncid) + + + ! Interpolate + call interpolate_wavespec_data (wave_spectrum_data, wave_spectrum) + + ! Save record number + oldrecnum = recnum + + if (local_debug) then + if (my_task == master_task) write (nu_diag,*) & + 'wave_spec_data ',spec_file + if (my_task.eq.master_task) & + write (nu_diag,*) 'maxrec',maxrec + write (nu_diag,*) 'days_per_year', days_per_year + + endif ! local debug + + end subroutine wave_spec_data + !======================================================================= ! initial snow aging lookup table diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 39f250f98..a2fa268ce 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -2086,11 +2086,11 @@ subroutine input_data if (trim(wave_spec_type) == 'none') then tmpstr2 = ' : no wave data provided, no wave-ice interactions' elseif (trim(wave_spec_type) == 'profile') then - tmpstr2 = ' : use fixed dummy wave spectrum for testing' + tmpstr2 = ' : use fixed dummy wave spectrum for testing, sea surface height generated using constant phase (1 iteration of wave fracture)' elseif (trim(wave_spec_type) == 'constant') then - tmpstr2 = ' : constant wave spectrum data file provided for testing' + tmpstr2 = ' : wave spectrum data file provided, sea surface height generated using constant phase (1 iteration of wave fracture)' elseif (trim(wave_spec_type) == 'random') then - tmpstr2 = ' : wave data file provided, sea surface height generated using random number' + tmpstr2 = ' : wave spectrum data file provided, sea surface height generated using random number (multiple iterations of wave fracture to convergence)' else tmpstr2 = ' : unknown value' endif From f88a2fef3b02f654bdc466a19829278114926b4f Mon Sep 17 00:00:00 2001 From: lettie-roach Date: Tue, 18 Oct 2022 09:39:48 -0400 Subject: [PATCH 4/4] Make FSD history output use aice, not aice_init --- .../cicedynB/analysis/ice_history_fsd.F90 | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_history_fsd.F90 b/cicecore/cicedynB/analysis/ice_history_fsd.F90 index bf7fa40b4..b52db4e05 100644 --- a/cicecore/cicedynB/analysis/ice_history_fsd.F90 +++ b/cicecore/cicedynB/analysis/ice_history_fsd.F90 @@ -294,7 +294,7 @@ subroutine accum_hist_fsd (dt, iblk) use ice_constants, only: c0, c1, c2, c4 use ice_history_shared, only: a2D, a3Df, a4Df, nfsd_hist, & ncat_hist, accum_hist_field, n3Dacum, n4Dscum - use ice_state, only: trcrn, aicen_init, vicen, aice_init + use ice_state, only: trcrn, aicen, vicen, aice use ice_arrays_column, only: wave_sig_ht, floe_rad_c, floe_binwidth, & d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld @@ -345,7 +345,7 @@ subroutine accum_hist_fsd (dt, iblk) worka(i,j) = c0 do n = 1, ncat_hist do k = 1, nfsd_hist - worka(i,j) = worka(i,j) + aicen_init(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) + worka(i,j) = worka(i,j) + aicen(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) end do end do end do @@ -360,7 +360,7 @@ subroutine accum_hist_fsd (dt, iblk) workb = c0 do n = 1, ncat_hist do k = 1, nfsd_hist - workc = aicen_init(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) & + workc = aicen(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) & / (c4*floeshape*floe_rad_c(k)**2) ! number-mean radius worka(i,j) = worka(i,j) + workc * floe_rad_c(k) @@ -383,7 +383,7 @@ subroutine accum_hist_fsd (dt, iblk) workb = c0 do n = 1, ncat_hist do k = 1, nfsd_hist - workb = workb + aicen_init(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) + workb = workb + aicen(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) worka(i,j) = worka(i,j) + vicen(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) end do end do @@ -401,12 +401,12 @@ subroutine accum_hist_fsd (dt, iblk) do j = 1, ny_block do i = 1, nx_block worka(i,j) = c0 - if (aice_init(i,j,iblk) > puny) then + if (aice(i,j,iblk) > puny) then do k = 1, nfsd_hist do n = 1, ncat_hist worka(i,j) = worka(i,j) & + (trcrn(i,j,nt_fsd+k-1,n,iblk) * floe_rad_c(k) & - * aicen_init(i,j,n,iblk)/aice_init(i,j,iblk)) + * aicen(i,j,n,iblk)/aice(i,j,iblk)) end do end do endif @@ -419,12 +419,12 @@ subroutine accum_hist_fsd (dt, iblk) do j = 1, ny_block do i = 1, nx_block worka(i,j) = c0 - if (aice_init(i,j,iblk) > puny) then + if (aice(i,j,iblk) > puny) then do k = 1, nfsd_hist do n = 1, ncat_hist worka(i,j) = worka(i,j) & + (c8*floeshape*trcrn(i,j,nt_fsd+k-1,n,iblk)*floe_rad_c(k) & - *aicen_init(i,j,n,iblk)/(c4*floeshape*floe_rad_c(k)**2 *aice_init(i,j,iblk))) + *aicen(i,j,n,iblk)/(c4*floeshape*floe_rad_c(k)**2 *aice(i,j,iblk))) end do end do endif @@ -445,7 +445,7 @@ subroutine accum_hist_fsd (dt, iblk) worke(i,j,k)=c0 do n = 1, ncat_hist worke(i,j,k) = worke(i,j,k) + (trcrn(i,j,nt_fsd+k-1,n,iblk) & - * aicen_init(i,j,n,iblk)/floe_binwidth(k)) + * aicen(i,j,n,iblk)/floe_binwidth(k)) end do end do end do @@ -479,7 +479,7 @@ subroutine accum_hist_fsd (dt, iblk) do j = 1, ny_block do i = 1, nx_block workd(i,j,k,n) = trcrn(i,j,nt_fsd+k-1,n,iblk) & - * aicen_init(i,j,n,iblk)/floe_binwidth(k) + * aicen(i,j,n,iblk)/floe_binwidth(k) end do end do end do