From 937c56d90f04bc81e1a39b8a969d17bc642c54a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Dupont?= Date: Fri, 14 Apr 2023 23:16:30 +0000 Subject: [PATCH] ice_init: add NetCDF coldstart capability This new feature consists of reading the following fields from a NetCDF input file (currently hardcoded to 'init_cice.nc') on the same grid as the model: * aice (ice concentration) * hice (ice thickness) * hsnow (snow thicknes) * Tair (surface air temperature) This is done in ice_init::init_state, just prior to calling 'set_state_var'. Then, in 'set_state_var', when computing the list of cells to initialize (arrays indx[ij]), we add a condition on 'ice_data_type' and choose all ocean points with non-zero concentration (since at that point 'aice' is already read from the NetCDF file). Arrays 'aice' and 'hice' are then used to distribute the volume among the different existing thickness category, using a procedure developed by Mathieu Chevalier in 2015 (ice_data_type='initnc_sprd'), or to set the volume of the closest category (ice_data_type='initnc'). Note that we need to make sure not to overwrite 'aicen', 'vicen' and 'vsnon' in the following loop on 'ncat'. Reading 'Tair' allows us to initialize the vertical temperature profile in the ice (in icepack_init_trcr) to more realistic values. In the future, we might want to also read sss and sst, and adjust the latter to ensure consistency between sss and the freezing temperature Tf (taking the maximum of Tf and sst). This is left for future work, but let's already add 'sst' as in 'inout' argument in 'readnc'. This way the existing values are left unmodified for now. Co-authored-by: Philippe Blain --- cicecore/cicedynB/general/ice_init.F90 | 313 ++++++++++++++++++++++++- 1 file changed, 310 insertions(+), 3 deletions(-) diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 8dab6e8fa..a0bea838c 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -2267,6 +2267,7 @@ subroutine init_state use ice_state, only: trcr_depend, aicen, trcrn, vicen, vsnon, & aice0, aice, vice, vsno, trcr, aice_init, bound_state, & n_trcr_strata, nt_strata, trcr_base, uvel, vvel + use ice_forcing, only: ice_data_type integer (kind=int_kind) :: & ilo, ihi , & ! physical domain indices @@ -2457,6 +2458,11 @@ subroutine init_state nt_strata (nt_ipnd,1) = nt_apnd ! on melt pond area endif + if (trim(ice_data_type(1:6)) == 'initnc') then + ! note that the volume arrays for ice and snow + ! are actually read as thickness, we are just using the arrays as placeholders + call readnc(aice, vice, vsno, Tair, sst) + endif !----------------------------------------------------------------- ! Set state variables !----------------------------------------------------------------- @@ -2484,7 +2490,9 @@ subroutine init_state salinz(:,:,:, iblk), Tmltz(:,:,:, iblk), & aicen(:,:, :,iblk), trcrn(:,:,:,:,iblk), & vicen(:,:, :,iblk), vsnon(:,:, :,iblk), & - uvel (:,:, iblk), vvel (:,:, iblk)) + uvel (:,:, iblk), vvel (:,:, iblk), & + aice (:,:, iblk), vice (:,:, iblk), & + vsno (:,:, iblk)) enddo ! iblk !$OMP END PARALLEL DO @@ -2562,7 +2570,9 @@ subroutine set_state_var (nx_block, ny_block, & salinz, Tmltz, & aicen, trcrn, & vicen, vsnon, & - uvel, vvel) + uvel, vvel, & + aice, hice, & + hsno) use ice_arrays_column, only: hin_max use ice_domain_size, only: nilyr, nslyr, nx_global, ny_global, ncat @@ -2589,7 +2599,10 @@ subroutine set_state_var (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & Tair , & ! air temperature (K) Tf , & ! freezing temperature (C) - sst ! sea surface temperature (C) + sst , & ! sea surface temperature (C) + aice , & ! ice total concentration + hice , & ! ice thickness (m) + hsno ! snow thickness (m) real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr), intent(in) :: & salinz , & ! initial salinity profile @@ -2799,6 +2812,34 @@ subroutine set_state_var (nx_block, ny_block, & endif ! rectgrid + ! ---------------------------------------------------------- + ! interface for choosing aicen,vicen,vsnon based on averaged + ! fields of aice,hice,hsno + ! ---------------------------------------------------------- + if (ice_data_type(1:6) == 'initnc') then + ! find valid ice points + icells = 0 + do j = jlo, jhi + do i = ilo, ihi + if ( tmask(i,j) .and. aice(i,j) > c0 ) then + icells = icells + 1 + indxi(icells) = i + indxj(icells) = j + endif ! tmask + enddo ! i + enddo ! j + + if (trim(ice_data_type) == 'initnc_sprd') then + ! spread mean values to ITD + call spread_avg_itd(icells, indxi, indxj, & + aice, hice, hsno, aicen, vicen, vsnon) + else + ! put the mean in closest category ITD + call put_avg_itd(icells, indxi, indxj, & + aice, hice, hsno, aicen, vicen, vsnon) + endif + endif + do n = 1, ncat ! ice volume, snow volume @@ -2806,6 +2847,8 @@ subroutine set_state_var (nx_block, ny_block, & i = indxi(ij) j = indxj(ij) + ! avoid overwriting aicen,vicen,vsnon if already defined + if (ice_data_type(1:6) /= 'initnc') then aicen(i,j,n) = ainit(n) if (trim(ice_data_type) == 'box2001') then @@ -2842,6 +2885,7 @@ subroutine set_state_var (nx_block, ny_block, & vicen(i,j,n) = hinit(n) * ainit(n) ! m endif vsnon(i,j,n) = min(aicen(i,j,n)*hsno_init,p2*vicen(i,j,n)) + endif ! on ice_data_type /= initnc call icepack_init_trcr(Tair = Tair(i,j), Tf = Tf(i,j), & Sprofile = salinz(i,j,:), & @@ -3018,6 +3062,269 @@ end subroutine boxslotcyl_data_vel !======================================================================= +! Read single-category state variables from netcdf +! +! author: FD 2012-11 (from initnc) + + subroutine readnc(aice, hice, hsno, Tair, sst) + + use ice_fileunits, only: nu_diag + use ice_read_write + use ice_constants, only: field_loc_center, field_type_scalar + use ice_domain_size, only: max_blocks + use ice_blocks, only: nx_block, ny_block + + ! arguments + real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks), intent(inout) :: & + Tair , & ! air temperature (K) + sst , & ! sea surface temperature (C) + aice , & ! ice total concentration + hice , & ! ice thickness (m) + hsno ! snow thickness (m) + + ! locals + integer (kind=int_kind) :: & + i, j, & ! index inside block + iblk, & ! block index + n, & ! ice category index + it, & ! tracer index + fid_init ! file id for netCDF init file + + character (char_len) :: & + fieldname ! field name in netcdf file + character (char_len_long) :: & ! input data file names + init_file = 'init_cice.nc' + logical (kind=log_kind) :: diag=.true. + + if (my_task == master_task) then + + write (nu_diag,*) ' ' + write (nu_diag,*) 'Initial ice file: ', trim(init_file) + call flush(nu_diag) + + call ice_open_nc(init_file,fid_init) + + endif + + fieldname='aice' + + call ice_read_nc(fid_init,1,fieldname,aice,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname='hice' + + call ice_read_nc(fid_init,1,fieldname,hice,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname='hsnow' + + call ice_read_nc(fid_init,1,fieldname,hsno,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + fieldname='tair' + + call ice_read_nc(fid_init,1,fieldname,Tair,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + +! for later +! fieldname='sst' +! +! call ice_read_nc(fid_init,1,fieldname,sst ,diag, & +! field_loc=field_loc_center, & +! field_type=field_type_scalar) +! + if (my_task == master_task) then + call ice_close_nc(fid_init) + write(nu_diag,*) 'closing file ',TRIM(init_file) + call flush(nu_diag) + endif + end subroutine readnc + +!======================================================================= + +! Initialize from single-category state variables +! +! author: FD 2012-11 (from initnc) + + subroutine put_avg_itd(icells, indxi, indxj, aice, hice, hsno, aicen, vicen, vsnon) + + use ice_constants, only: c0 + use ice_domain_size, only: max_blocks, ncat + use ice_blocks, only: nx_block, ny_block + use ice_arrays_column, only: hin_max + + ! arguments + integer, intent(in) :: icells + + integer (kind=int_kind), dimension(nx_block*ny_block), intent(in) :: & + indxi, indxj ! compressed indices for cells with aicen > puny + + real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & + aice, hice, hsno ! input array (real, 8-byte) + + real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), intent(out) :: & + aicen, vicen, vsnon ! output array (real, 8-byte) + + ! locals + integer (kind=int_kind) :: & + ij, & ! compact index + i, j, & ! index inside block + n ! ice category index + + aicen(:,:,:) = c0 + vicen(:,:,:) = c0 + vsnon(:,:,:) = c0 + + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + + !----------------------------------------------------------------- + ! Place ice at given thickness category + do n = 1, ncat + if ( hin_max(n-1) <= hice(i,j) .and. & + & hice(i,j) < hin_max(n) ) then + aicen(i,j,n) = aice(i,j) + vicen(i,j,n) = aice(i,j) * hice(i,j) ! m + vsnon(i,j,n) = aice(i,j) * hsno(i,j) ! m + endif + enddo + enddo + + end subroutine put_avg_itd + +!======================================================================= + +! Initialize from single-category state variables and spread ITD +! +! author: FD 2012-11 (from initnc) + + subroutine spread_avg_itd(icells, indxi, indxj, aice, hice, hsno, aicen, vicen, vsnon) + + use ice_constants, only: c0, c1, c2, c3, p2, p5 + use ice_domain_size, only: max_blocks, ncat + use ice_blocks, only: nx_block, ny_block + use ice_arrays_column, only: hin_max + + ! arguments + integer, intent(in) :: icells + + integer (kind=int_kind), dimension(nx_block*ny_block), intent(in) :: & + indxi, indxj ! compressed indices for cells with aicen > puny + + real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: & + aice, hice, hsno ! input array (real, 8-byte) + + real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), intent(out) :: & + aicen, vicen, vsnon ! output array (real, 8-byte) + + ! locals + integer (kind=int_kind) :: & + ij, & ! compact index + i, j, & ! index inside block + n ! ice category index + + real (kind=dbl_kind) :: & + zfbar, zhbar, zhsno, zhinf, & + za, zb, zhm1, zhp1, zjk, zfk, zhk, zhup, zfseuil=0.85 + + aicen(:,:,:) = c0 + vicen(:,:,:) = c0 + vsnon(:,:,:) = c0 + + !----------------------------------------------------------------- + ! Place ice where ocean surface is cold. + ! Note: If SST is not read from a file, then the ocean is assumed + ! to be at its freezing point everywhere, and ice will + ! extend to the prescribed edges. + !----------------------------------------------------------------- + + do ij = 1, icells + i = indxi(ij) + j = indxj(ij) + + !----------------------------------------------------------------- + ! Spread ice through different thickness categories (Mathieu Chevalier's function) +! MCH spread thickness/concentration over all categories + + zfbar = aice(i,j) + zhbar = aice(i,j) * hice(i,j) + zhsno = hsno(i,j) + if ( zfbar <= zfseuil .and. zfbar > 0.01_dbl_kind .and. zhbar > 0.1_dbl_kind) then + ! Assume a linear distribution in the MIZ + zhinf = 3._dbl_kind*zhbar / zfbar + za=2._dbl_kind*zfbar / zhinf + zb=za / zhinf + do n=1,ncat + zhm1=hin_max(n-1) + zhp1=hin_max(n) + if (zhm1 .ge. zhinf ) then + zfk=0._dbl_kind + if (n.eq.ncat) then + zhk=(zhm1+1.)*zfk + else + zhk=0.5_dbl_kind*(zhp1+zhm1)*zfk + end if + else + if (zhinf.le.zhp1) then + zhup=zhinf + else + zhup=zhp1 + end if + zfk = za*(zhup - zhm1)-zb*(zhup*zhup-zhm1*zhm1)/2._dbl_kind + zhk = za*(zhup*zhup - zhm1*zhm1)/2._dbl_kind-zb*(zhup*zhup*zhup - zhm1*zhm1*zhm1)/3._dbl_kind + end if + aicen(i,j,n)=zfk + vicen(i,j,n)=zhk + end do + elseif ( zfbar > zfseuil .and. zhbar > 0.1_dbl_kind ) then +! Assume a parabolic distribution in compact pack ice + zhinf = 2._dbl_kind*zhbar / zfbar + za = 6._dbl_kind * zfbar / (zhinf*zhinf*zhinf) + do n=1,ncat + zhm1=hin_max(n-1) + zhp1=hin_max(n) + if (zhm1 .ge. zhinf ) then + zfk=0._dbl_kind + if (n.eq.ncat) then + zhk=(zhm1+1.)*zfk + else + zhk=0.5_dbl_kind*(zhp1+zhm1)*zfk + end if + else + if (zhinf.le.zhp1) then + zhup=zhinf + else + zhup=zhp1 + end if + zfk = za*( zhinf*(zhup*zhup - zhm1*zhm1) / 2._dbl_kind - ( zhup*zhup*zhup - zhm1*zhm1*zhm1 ) / 3._dbl_kind ) + zhk = za*( zhinf*(zhup*zhup*zhup - zhm1*zhm1*zhm1) / 3._dbl_kind - ( zhup*zhup*zhup*zhup - zhm1*zhm1*zhm1*zhm1 ) / 4._dbl_kind ) + end if + aicen(i,j,n)=zfk + vicen(i,j,n)=zhk + end do + else + do n = 1, ncat + !----------------------------------------------------------------- + ! Place ice at given thickness category + + if ( hin_max(n-1) <= hice(i,j) .and. & + & hice(i,j) < hin_max(n) ) then + aicen(i,j,n) = zfbar + vicen(i,j,n) = zhbar ! m + endif + enddo + end if + vsnon(i,j,:) = min(zhsno * aicen(i,j,:), vicen(i,j,:)*p2) + enddo + + end subroutine spread_avg_itd +!======================================================================= + end module ice_init !=======================================================================