Merge branch 'add-ice-mask-restoring' into 'concepts/main' (!37)

add a functionality for restoring ice over a mask

Closes #23

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 =', &

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
 
 !=======================================================================

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)

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)