diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index 61936809c..17cc68b42 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -121,7 +121,7 @@ subroutine input_data algo_nonlin, fpfunc_andacc, dim_andacc, reltol_andacc, & damping_andacc, start_andacc, use_mean_vrel, ortho_type use ice_transport_driver, only: advection, conserv_check - use ice_restoring, only: restore_ice, restore_ice_to + use ice_restoring, only: restore_ice, restore_ice_to, restore_ice_mask use ice_timers, only: timer_stats use ice_memusage, only: memory_stats use ice_fileunits, only: goto_nml @@ -280,7 +280,7 @@ subroutine input_data fyear_init, ycycle, wave_spec_file,restart_coszen, & atm_data_dir, ocn_data_dir, bgc_data_dir, & atm_data_format, ocn_data_format, rotate_wind, & - oceanmixed_file, atm_data_version + oceanmixed_file, atm_data_version, restore_ice_mask !----------------------------------------------------------------- ! default values @@ -549,6 +549,7 @@ subroutine input_data trestore = 90 ! restoring timescale, days (0 instantaneous) restore_ice = .false. ! restore ice state on grid edges if true restore_ice_to = 'initial' ! restore ice state to initial ice state ('initial'), hardcoded values ('defined') or zero ('zero') + restore_ice_mask= .false. ! by default no ice restoring depending on mask, set to true in namelist if you want to restore depending on mask read from netcdf file debug_forcing = .false. ! true writes diagnostics for input forcing latpnt(1) = 90._dbl_kind ! latitude of diagnostic point 1 (deg) @@ -1156,6 +1157,7 @@ subroutine input_data call broadcast_scalar(trestore, master_task) call broadcast_scalar(restore_ice, master_task) call broadcast_scalar(restore_ice_to, master_task) + call broadcast_scalar(restore_ice_mask, master_task) call broadcast_scalar(debug_forcing, master_task) call broadcast_array (latpnt(1:2), master_task) call broadcast_array (lonpnt(1:2), master_task) @@ -2496,6 +2498,13 @@ subroutine input_data write(nu_diag,1031) ' restore_ice_to = ', restore_ice_to if (restore_ice .or. restore_ocn) & write(nu_diag,1021) ' trestore = ', trestore + if (restore_ice_mask) then + if ( trim(restore_ice_to) /= 'zero' ) then + if (my_task == master_task) write(nu_diag,*) subname//' ERROR: inconsistency between restore_ice_to ',trim(restore_ice_to), ' and restore_ice_mask=T' + abort_list = trim(abort_list)//":26" + endif + write(nu_diag,1011) ' restore_ice_mask = ', restore_ice_mask + endif write(nu_diag,*) ' ' write(nu_diag,'(a31,2f8.2)') 'Diagnostic point 1: lat, lon =', & diff --git a/cicecore/cicedyn/infrastructure/ice_restoring.F90 b/cicecore/cicedyn/infrastructure/ice_restoring.F90 index dd09ad863..ee72ad434 100644 --- a/cicecore/cicedyn/infrastructure/ice_restoring.F90 +++ b/cicecore/cicedyn/infrastructure/ice_restoring.F90 @@ -23,7 +23,7 @@ module ice_restoring implicit none private - public :: ice_HaloRestore_init, ice_HaloRestore + public :: ice_HaloRestore_init, ice_HaloRestore, ice_MaskRestore, ice_read_restoring_mask logical (kind=log_kind), public :: & restore_ice ! restore ice state if true @@ -34,6 +34,9 @@ module ice_restoring ! 'defined' => restore to hardcoded values ! 'zero' => restore to zero + logical (kind=log_kind), public :: & + restore_ice_mask ! restore ice over mask if true + !----------------------------------------------------------------- ! state of the ice for each category !----------------------------------------------------------------- @@ -46,6 +49,9 @@ module ice_restoring real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable, public :: & trcrn_rest ! tracers + real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & + iceRestMask ! ice restoring mask + !======================================================================= contains @@ -81,7 +87,7 @@ subroutine ice_HaloRestore_init character(len=*), parameter :: subname = '(ice_HaloRestore_init)' - if (.not. restore_ice) return + if (.not. restore_ice .and. .not.restore_ice_mask) return call icepack_query_tracer_sizes(ntrcr_out=ntrcr) call icepack_warnings_flush(nu_diag) @@ -741,6 +747,156 @@ end subroutine ice_HaloRestore !======================================================================= +! This subroutine is intended for restoring the ice state to desired +! values in all cells where restoring mask is true +! Note: This routine will need to be modified for nghost > 1. +! We assume padding occurs only on east and north edges. + + subroutine ice_MaskRestore + + use ice_blocks, only: block, get_block, nblocks_x, nblocks_y + use ice_calendar, only: dt + use ice_domain, only: nblocks, blocks_ice + +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + integer (int_kind) :: & + i,j,iblk,nt,n, &! dummy loop indices + ilo,ihi,jlo,jhi, &! beginning and end of physical domain + ntrcr + + type (block) :: & + this_block ! block info for current block + + real (dbl_kind) :: & + secday, &! + ctime ! dt/trest + + character(len=*), parameter :: subname = '(ice_MaskRestore)' + + call ice_timer_start(timer_bound) + call icepack_query_parameters(secday_out=secday) + call icepack_query_tracer_sizes(ntrcr_out=ntrcr) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + +!----------------------------------------------------------------------- +! +! Initialize +! +!----------------------------------------------------------------------- + + ! for now, use same restoring constant as for SST + if (trestore == 0) then + trest = dt ! use data instantaneously + else + trest = real(trestore,kind=dbl_kind) * secday ! seconds + endif + ctime = dt/trest + +!----------------------------------------------------------------------- +! +! Restore values in cells surrounding the grid +! +!----------------------------------------------------------------------- + + !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block, & + !$OMP i,j,n,nt,ibc,ice_restore_masknpad) + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + do n = 1, ncat + do j = 1, ny_block + do i = 1, nx_block + aicen(i,j,n,iblk) = aicen(i,j,n,iblk) & + + (aicen_rest(i,j,n,iblk)-aicen(i,j,n,iblk))*ctime * iceRestMask(i,j,iblk) + vicen(i,j,n,iblk) = vicen(i,j,n,iblk) & + + (vicen_rest(i,j,n,iblk)-vicen(i,j,n,iblk))*ctime * iceRestMask(i,j,iblk) + vsnon(i,j,n,iblk) = vsnon(i,j,n,iblk) & + + (vsnon_rest(i,j,n,iblk)-vsnon(i,j,n,iblk))*ctime * iceRestMask(i,j,iblk) +! the following lines are kept commented until a physical value for tracers are given for the interior restoring fields +! do nt = 1, ntrcr +! trcrn(i,j,nt,n,iblk) = trcrn(i,j,nt,n,iblk) & +! + (trcrn_rest(i,j,nt,n,iblk)-trcrn(i,j,nt,n,iblk))*ctime * iceRestMask(i,j,iblk) +! enddo + enddo + enddo + enddo + enddo ! iblk + !$OMP END PARALLEL DO + + call ice_timer_stop(timer_bound) + + end subroutine ice_MaskRestore + +!======================================================================= + +! Read restoring mask from netcdf +! +! author: FD 2023-09 + + subroutine ice_read_restoring_mask + + use ice_communicate, only: my_task, master_task + use ice_fileunits, only: nu_diag + use ice_read_write + use ice_constants, only: field_loc_center, field_type_scalar + use ice_blocks, only: nx_block, ny_block + use ice_grid, only: hm + + ! locals + integer (kind=int_kind) :: & + fid_mask ! file id for netCDF file + + character (char_len) :: & + fieldname ! field name in netcdf file + + character (char_len_long) :: & ! input data file names + mask_file = 'ice_mask_restore.nc' + + logical (kind=log_kind) :: diag=.true. + + + allocate (iceRestMask(nx_block,ny_block,max_blocks)) + + + if (my_task == master_task) then + + write (nu_diag,*) ' ' + write (nu_diag,*) 'Restoring mask file: ', trim(mask_file) + call flush(nu_diag) + + call ice_open_nc(mask_file,fid_mask) + + endif + + fieldname='mask' + + call ice_read_nc(fid_mask,1,fieldname,iceRestMask,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + if (my_task == master_task) then + call ice_close_nc(fid_mask) + write(nu_diag,*) 'closing file ',TRIM(mask_file) + call flush(nu_diag) + endif + + ! ensure mask is only one over wet points + iceRestMask(:,:,:) = iceRestMask(:,:,:) * hm(:,:,:) + + end subroutine ice_read_restoring_mask +!======================================================================= + end module ice_restoring !======================================================================= diff --git a/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 b/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 index dbf1f82cc..5f5775430 100644 --- a/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 +++ b/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 @@ -86,7 +86,7 @@ subroutine cice_init use ice_init, only: input_data, init_state use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers use ice_kinds_mod - use ice_restoring, only: ice_HaloRestore_init + use ice_restoring, only: ice_HaloRestore_init, restore_ice_mask, ice_read_restoring_mask use ice_timers, only: timer_total, init_ice_timers, ice_timer_start use ice_transport_driver, only: init_transport use lib_mpp, only: mpi_comm_opa ! NEMO MPI communicator @@ -162,6 +162,8 @@ subroutine cice_init call init_state ! initialize the ice state call init_transport ! initialize horizontal transport call ice_HaloRestore_init ! restored boundary conditions + if (restore_ice_mask) & + call ice_read_restoring_mask! read restoring mask if restore_ice_mask is true call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & wave_spec_out=wave_spec) diff --git a/cicecore/drivers/direct/nemo_concepts/CICE_RunMod.F90 b/cicecore/drivers/direct/nemo_concepts/CICE_RunMod.F90 index cf0b55960..6623d893b 100644 --- a/cicecore/drivers/direct/nemo_concepts/CICE_RunMod.F90 +++ b/cicecore/drivers/direct/nemo_concepts/CICE_RunMod.F90 @@ -153,7 +153,7 @@ subroutine ice_step write_restart_bgc, write_restart_hbrine use ice_restart_driver, only: dumpfile use ice_restart_shared, only: restart, runtype - use ice_restoring, only: restore_ice, ice_HaloRestore + use ice_restoring, only: restore_ice, ice_HaloRestore, restore_ice_mask, ice_MaskRestore use ice_state, only: trcrn use ice_step_mod, only: prep_radiation, step_therm1, step_therm2, & update_state, step_dyn_horiz, step_dyn_ridge, step_radiation, & @@ -203,6 +203,7 @@ subroutine ice_step !----------------------------------------------------------------- if (restore_ice) call ice_HaloRestore + if (restore_ice_mask) call ice_MaskRestore !--------------------------------------------------------------- ! 1st time step: initialize shortwave again (with correct forcing)