118 ostia ice ancils

applications/lfric_atm/metadata/grid_def_main.xml

@@ -101,6 +101,11 @@
     <axis axis_ref="snow_layers_and_tiles"/>
     <axis axis_ref="one_time_axis"/>
   </grid>
+  <grid id="one_time_seaice_grid">
+    <domain domain_ref="face"/>
+    <axis axis_ref="one_time_axis"/>
+    <axis axis_ref="sea_ice_categories"/>
+  </grid>
   <grid id="one_time_tile_grid">
     <domain domain_ref="face"/>
     <axis axis_ref="land_tiles"/>

interfaces/jules_interface/source/algorithm/init_surface_fields_alg_mod.x90

@@ -154,7 +154,7 @@ contains
                                              canopy_height, tile_temperature, &
                                              tile_lw_grey_albedo ),           &
                  ! Sea ice conductivity will use either this fixed value or
-                 ! it will get overwritten by the coupler
+                 ! it will get overwritten by the coupler or process_inputs
                  setval_c(sea_ice_conductivity, default_conductivity ),       &
                  ! Initialise melt ponds to zero
                  setval_c(melt_pond_fraction, 0.0_r_def ),                    &

interfaces/jules_interface/source/algorithm/process_inputs_alg_mod.x90

@@ -15,7 +15,10 @@ module process_inputs_alg_mod
 
   use initialization_config_mod,  only: sst_source, sst_source_start_dump, &
                                         init_option, init_option_fd_start_dump, &
-                                        snow_source, snow_source_surf
+                                        snow_source, snow_source_surf, &
+                                        sea_ice_source, &
+                                        sea_ice_source_start_dump, &
+                                        sea_ice_source_surf
   use io_config_mod,              only: checkpoint_read
   use section_choice_config_mod,  only: surface, surface_jules, &
                                         radiation, radiation_socrates
@@ -26,6 +29,8 @@ module process_inputs_alg_mod
   use process_ssi_kernel_mod,     only: process_ssi_kernel_type
   use process_soil_kernel_mod,    only: process_soil_kernel_type
   use process_snow_kernel_mod,    only: process_snow_kernel_type
+  use initial_ice_conductivity_kernel_mod, &
+                                  only: initial_ice_conductivity_kernel_type
   use jules_snow_mod,             only: nsmax
   use jules_control_init_mod,     only: n_land_tile, &
                                         first_sea_tile, n_sea_tile, &
@@ -37,6 +42,8 @@ module process_inputs_alg_mod
   use set_topography_kernel_mod,  only: set_topography_kernel_type
   use sci_geometric_constants_mod,      &
                                   only: get_latitude_fv
+  use log_mod,                    only : log_event, log_scratch_space,  &
+                                         LOG_LEVEL_ERROR
 
   implicit none
 
@@ -70,6 +77,7 @@ contains
     type( field_type ), pointer :: land_tile_fraction => null()
     type( field_type ), pointer :: sea_ice_fraction   => null()
     type( field_type ), pointer :: sea_ice_thickness  => null()
+    type( field_type ), pointer :: sea_ice_conductivity  => null()
     type( field_type ), pointer :: tile_fraction      => null()
 
     type( field_type ), pointer :: tstar_sea        => null()
@@ -126,11 +134,19 @@ contains
 
     integer( i_def ) :: n_horizon_angle, n_horizon_layer
 
+    if (checkpoint_read) then
+      ! This algorithm should only be called on a new run
+      ! i.e. checkpoint_read==.false.)
+      ! Add a check here in case anyone adds code that breaks this
+      write(log_scratch_space,'(A)')                                         &
+         'process_input_args: should not be called if checkpoint_read==.true.'
+      call log_event(log_scratch_space, LOG_LEVEL_ERROR)
+    end if
+
     if (surface == surface_jules) then
 
       ! Only needed on cold start from UM
-      if (init_option == init_option_fd_start_dump .and. &
-           .not. checkpoint_read) then
+      if (init_option == init_option_fd_start_dump) then
         call ancil_fields%get_field('land_area_fraction', land_area_fraction)
         call ancil_fields%get_field('land_tile_fraction', land_tile_fraction)
         call fd_fields%get_field('land_tile_temp', land_tile_temp)
@@ -139,10 +155,11 @@ contains
         call surface_fields%get_field('canopy_water', canopy_water)
       end if
 
-      ! For coupled models the sea ice fraction comes from the input
-      ! dump but for atmosphere only models they are provided by an
-      ! ancillary file
-      if (l_couple_sea_ice) then
+      ! For coupled models the sea ice fraction usually comes from the input
+      ! dump. 
+      ! For atmosphere only models or coupled models with inland lakes it is 
+      ! provided by an ancillary file
+      if (sea_ice_source == sea_ice_source_start_dump) then
         call fd_fields%get_field('sea_ice_fraction', sea_ice_fraction)
       else
         call ancil_fields%get_field('sea_ice_fraction', sea_ice_fraction)
@@ -157,8 +174,7 @@ contains
       end if
       call surface_fields%get_field('tile_temperature', tile_temperature)
 
-      if (init_option == init_option_fd_start_dump .and. &
-           .not. checkpoint_read) then
+      if (init_option == init_option_fd_start_dump) then
         call soil_fields%get_field('soil_moist_sat', soil_moist_sat)
         call soil_fields%get_field('soil_moist_wilt', soil_moist_wilt)
         call soil_fields%get_field('mean_topog_index', mean_topog_index)
@@ -176,8 +192,8 @@ contains
         call soil_fields%get_field('frozen_soil_moisture', frozen_soil_moisture)
       end if
 
-      if ( (init_option == init_option_fd_start_dump .and. &
-           .not. checkpoint_read) .or. (snow_source == snow_source_surf)) then
+      if ( (init_option == init_option_fd_start_dump) .or. &
+            (snow_source == snow_source_surf)) then
         ! need to add .or. use_surf_analysis here
         call snow_fields%get_field('tile_snow_mass', tile_snow_mass)
         call snow_fields%get_field('tile_snow_rgrain', tile_snow_rgrain)
@@ -213,8 +229,7 @@ contains
       latitude => get_latitude_fv( W3, twod_mesh%get_id() )
 
       ! Surface fields
-      if (init_option == init_option_fd_start_dump .and. &
-           .not. checkpoint_read) then
+      if (init_option == init_option_fd_start_dump) then
         ! These only need processing when cold-starting from the UM
         call invoke(process_land_kernel_type( land_area_fraction,             &
                                             land_tile_fraction,               &
@@ -235,9 +250,20 @@ contains
                   multi_insert_kernel_type(tile_temperature, tstar_sea,       &
                                            first_sea_tile, n_sea_tile ))
 
+      ! In a coupled model with inland lakes we need to initialise the sea-ice
+      ! conductivity. This will be overwritten by values from SI3 except at
+      ! lake points.
+      if ((sea_ice_source == sea_ice_source_surf ) .and.                     &
+           (l_couple_sea_ice)) then
+           call surface_fields%get_field('sea_ice_conductivity',              &
+                                          sea_ice_conductivity)
+           call invoke(initial_ice_conductivity_kernel_type(sea_ice_fraction, &
+                                                         sea_ice_thickness,   &
+                                                         sea_ice_conductivity))
+      endif
+
       ! Soil fields
-      if (init_option == init_option_fd_start_dump .and. &
-           .not. checkpoint_read) then
+      if (init_option == init_option_fd_start_dump) then
         ! These only need processing when cold-starting from the UM
         call invoke(process_soil_kernel_type( soil_moist_sat,                 &
                                             soil_moist_wilt,                  &
@@ -256,8 +282,8 @@ contains
       end if
 
       ! Snow fields
-      if ((init_option == init_option_fd_start_dump .and. &
-           .not. checkpoint_read) .or. snow_source == snow_source_surf) then
+      if ((init_option == init_option_fd_start_dump) .or. &
+            (snow_source == snow_source_surf)) then
         ! These only need processing when cold-starting from the UM
         call invoke(multi_insert_kernel_type(tile_snow_mass, tile_snow_mass_in,&
                                               1, n_land_tile),                   &
@@ -284,8 +310,7 @@ contains
 
     if (radiation == radiation_socrates .and. &
         topography /= topography_flat .and. &
-        init_option == init_option_fd_start_dump .and. &
-        .not. checkpoint_read) then
+        init_option == init_option_fd_start_dump) then
 
       call ancil_fields%get_field('grad_x_orog', grad_x_orog)
       call ancil_fields%get_field('grad_y_orog', grad_y_orog)

interfaces/jules_interface/source/kernel/initial_ice_conductivity_kernel_mod.F90

@@ -0,0 +1,73 @@
+!-------------------------------------------------------------------------------
+! (c) Crown copyright 2020 Met Office. All rights reserved.
+! The file LICENCE, distributed with this code, contains details of the terms
+! under which the code may be used.
+!-------------------------------------------------------------------------------
+!> @brief Initialise Jules surface fields on tiles
+!> @details Non-standard Surface fields (pseudo-levels) aren't as yet not
+!>  implemented in LFRic. As an interim measure Higher-order W3 fields have
+!>  been used to mimic psuedo-level field behaviour. This code is written
+!>  based on this interim measure and will need to be updated when
+!>  suitable infrastructure is available (Ticket #2081)
+module initial_ice_conductivity_kernel_mod
+  use argument_mod,  only: arg_type,                       &
+                           GH_FIELD, GH_REAL, GH_INTEGER,  &
+                           GH_WRITE, CELL_COLUMN,          &
+                           ANY_DISCONTINUOUS_SPACE_1,      &
+                           ANY_DISCONTINUOUS_SPACE_2
+  use constants_mod, only: r_def, i_def
+  use kernel_mod,    only: kernel_type
+  use jules_control_init_mod, only: n_sea_ice_tile, first_sea_ice_tile
+  use jules_sea_seaice_config_mod,     only: therm_cond_sice => kappai
+  implicit none
+  private
+  !> Kernel metadata for Psyclone
+  type, public, extends(kernel_type) :: initial_ice_conductivity_kernel_type
+    private
+    type(arg_type) :: meta_args(3) = (/                                    &
+         arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), &
+         arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), &
+         arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1) &
+         /)
+    integer :: operates_on = CELL_COLUMN
+  contains
+    procedure, nopass :: initial_ice_conductivity_code
+  end type initial_ice_conductivity_kernel_type
+  public :: initial_ice_conductivity_code
+contains
+  !> @param[in]     nlayers              The number of layers
+  !> @param[in,out] sea_ice_fraction     Sea Ice Fractions on categories
+  !> @param[in,out] sea_ice_thickness    Sea Ice Thickness on categories
+  !> @param[in,out] sea_ice_conductivity Sea Ice Conductivity on categories
+  !> @param[in]     ndf_ice              Number of DOFs per cell for tiles
+  !> @param[in]     undf_ice             Number of total DOFs for tiles
+  !> @param[in]     map_ice              Dofmap for cell for surface tiles
+  subroutine initial_ice_conductivity_code(nlayers,                           &
+                                           sea_ice_fraction,                  &
+                                           sea_ice_thickness,                 &
+                                           sea_ice_conductivity,              &
+                                           ndf_ice, undf_ice, map_ice)
+    implicit none
+    ! Arguments
+    integer(kind=i_def), intent(in) :: nlayers
+    integer(kind=i_def), intent(in) :: ndf_ice, undf_ice
+    integer(kind=i_def), intent(in) :: map_ice(ndf_ice)
+    real(kind=r_def), intent(inout) :: sea_ice_fraction(undf_ice)
+    real(kind=r_def), intent(inout) :: sea_ice_thickness(undf_ice)
+    real(kind=r_def), intent(inout) :: sea_ice_conductivity(undf_ice)
+    ! Internal variables
+    integer(kind=i_def) :: i
+    real(kind=r_def) :: min_ice_thick, max_ice_cond
+    !Taken from UM recon value
+    max_ice_cond = 25.0_r_def
+    min_ice_thick = 8.0_r_def * therm_cond_sice/max_ice_cond
+    do i=0,n_sea_ice_tile-1
+        if (sea_ice_thickness(map_ice(1)+i) >= min_ice_thick) then
+            sea_ice_conductivity(map_ice(1)+i) = (8.0_r_def * therm_cond_sice) &
+                                              / sea_ice_thickness(map_ice(1)+i)
+        else
+            sea_ice_conductivity(map_ice(1)+i) = max_ice_cond
+        endif
+    end do
+  end subroutine initial_ice_conductivity_code
+end module initial_ice_conductivity_kernel_mod

interfaces/jules_interface/source/support/jules_physics_init_mod.f90

@@ -355,6 +355,9 @@ subroutine jules_physics_init()
     ! l_ctile is implicitly true by design of LFRic and should not be changed
     l_ctile              = .true.
     l_iceformdrag_lupkes = l_iceformdrag_lupkes_in
+    ! l_saldep_freeze should always be set to false as it no longer affects 
+    ! the coupled model except at lake points (which aren't coupled).
+    l_saldep_freeze       = .false.
     l_stability_lupkes   = l_stability_lupkes_in
     l_sice_heatflux      = l_sice_heatflux_in
     ! Code has not been included to support this being false as configurations
@@ -376,7 +379,6 @@ subroutine jules_physics_init()
       l_sice_meltponds_cice = .true.
       l_tstar_sice_new      = .false.
       l_cice_alb            = .true.
-      l_saldep_freeze       = .true.
       l_sice_multilayers    = .true.
       l_sice_scattering     = .true.
       l_ssice_albedo        = .true.
@@ -388,7 +390,6 @@ subroutine jules_physics_init()
       l_sice_meltponds_cice = .false.
       l_tstar_sice_new      = .true.
       l_cice_alb            = .false.
-      l_saldep_freeze       = .false.
       l_sice_multilayers    = .false.
       l_sice_scattering     = .false.
       l_ssice_albedo        = .false.

rose-stem/app/lfric_atm/file/file_def_ancil_surf.xml

@@ -8,8 +8,8 @@
 
   <file id="sea_ice_ancil" name="sea_ice_ancil" mode="read" output_freq="1d" convention="UGRID" cyclic="true" enabled=".FALSE.">
     <field id="sea_ice_time" name="time" unit="hours" operation="once" axis_ref="one_time_axis"/>
-    <field id="sea_ice_fraction_data" name="sea_ice_area_fraction" long_name="sea_ice_fraction" unit="1" operation="instant" domain_ref="face" axis_ref="one_time_axis" />
-       <field id="sea_ice_thickness_data" name="sea_ice_thickness" long_name="sea_ice_thickness" unit="m" operation="instant" domain_ref="face" axis_ref="one_time_axis"/>
+    <field id="sea_ice_fraction_data" name="sea_ice_area_fraction" long_name="sea_ice_fraction" unit="1" operation="instant" grid_ref="one_time_seaice_grid" />
+    <field id="sea_ice_thickness_data" name="sea_ice_thickness" long_name="sea_ice_thickness" unit="m" operation="instant" grid_ref="one_time_seaice_grid"/>
   </file>
 
 </file_definition>

rose-stem/app/lfric_atm/opt/rose-app-da.conf

@@ -3,7 +3,7 @@ source=$ROSE_SUITE_DIR/app/lfric_atm/file/iodef_gal_nwp_cycling.xml
 
 [namelist:files]
 iau_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/20210324T0600Z_um2lfric_iau_000001'
-sea_ice_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/seaice_ugrid_postqa_fixed'
+sea_ice_ancil_path='/data/users/tim.graham/LFRIC_SEA_ICE_ANCIL/seaice_ugrid_postqa_fixed'
 snow_analysis_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/20210324T0600Z_glu_snow_um2lfric_fixed'
 sst_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/sst_ugrid_postqa_fixed'
 start_dump_directory='$BIG_DATA_DIR/start_dumps/basic-gal/yak/PR202'

rose-stem/app/lfric_atm/opt/rose-app-eda.conf

@@ -8,7 +8,7 @@ iau_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/20210324T0600Z_um2lfric_iau_00
 iau_pert_path='$BIG_DATA_DIR/IAU/Global/iau_pertinc_start'
 iau_sst_path='$BIG_DATA_DIR/IAU/Global/um2lfric_sstpert'
 iau_surf_path='$BIG_DATA_DIR/IAU/Global/um2lfric_landda'
-sea_ice_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/seaice_ugrid_postqa_fixed'
+sea_ice_ancil_path='/data/users/tim.graham/LFRIC_SEA_ICE_ANCIL/seaice_ugrid_postqa_fixed'
 snow_analysis_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/20210324T0600Z_glu_snow_um2lfric_fixed'
 sst_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/sst_ugrid_postqa_fixed'
 start_dump_directory='$BIG_DATA_DIR/start_dumps/basic-gal/yak/PR202'

rose-stem/app/lfric_atm/opt/rose-app-eda_jada.conf

@@ -4,7 +4,7 @@ source=$ROSE_SUITE_DIR/app/lfric_atm/file/iodef_gal_nwp_cycling.xml
 [namelist:files]
 iau_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/20210324T0600Z_um2lfric_iau_000001'
 iau_pert_path='$BIG_DATA_DIR/IAU/Global/iau_pertinc_start'
-sea_ice_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/seaice_ugrid_postqa_fixed'
+sea_ice_ancil_path='/data/users/tim.graham/LFRIC_SEA_ICE_ANCIL/seaice_ugrid_postqa_fixed'
 snow_analysis_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/20210324T0600Z_glu_snow_um2lfric_fixed'
 sst_ancil_path='$BIG_DATA_DIR/start_dumps/basic-gal/yak/sst_ugrid_postqa_fixed'
 start_dump_directory='$BIG_DATA_DIR/start_dumps/basic-gal/yak/PR202'

rose-stem/app/lfric_coupled_atmosphere/rose-app.conf

@@ -577,7 +577,7 @@ ls_option='none'
 model_eos_height=100
 n_orog_smooth=0
 read_w2h_wind=.false.
-sea_ice_source='ancillary'
+sea_ice_source='start_dump'
 snow_source='start_dump'
 sst_source='ancillary'
 w0_orography_mapping=.false.

science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf

@@ -3231,7 +3231,7 @@ help=The sea ice fraction (and thickness if amip_ice_thick is false) can either
     =The iodef file must be consistent with this setting, specifying
     = seaice from either a static ancillary file or a SURF-derived ancillary file.
 sort-key=01a
-values='ancillary', 'surf'
+values='ancillary','start_dump','surf'
 
 [namelist:initialization=snow_source]
 compulsory=true

science/gungho/source/driver/create_fd_prognostics_mod.f90

@@ -29,7 +29,9 @@ module create_fd_prognostics_mod
                                              ancil_option_updating,       &
                                              read_w2h_wind,               &
                                              sst_source,                  &
-                                             sst_source_start_dump
+                                             sst_source_start_dump,       &
+                                             sea_ice_source,               &
+                                             sea_ice_source_start_dump
   use nlsizes_namelist_mod,           only : sm_levels
   use jules_control_init_mod,         only : n_land_tile, n_sea_ice_tile
   use jules_physics_init_mod,         only : snow_lev_tile
@@ -261,7 +263,7 @@ subroutine create_fd_prognostics( mesh, twod_mesh, fd_field_collection, &
 
       ! For coupled models get the sea ice fraction and thickness from the
       ! dump
-      if (l_couple_sea_ice) then
+      if (l_couple_sea_ice .and. (sea_ice_source == sea_ice_source_start_dump)) then
          call setup_ancil_field("sea_ice_fraction", depository, &
                              fd_field_collection, mesh, twod_mesh, &
                              twod=.true., ndata=n_sea_ice_tile)

science/gungho/source/driver/gungho_init_fields_mod.X90

@@ -917,11 +917,11 @@ subroutine create_model_data( modeldb,      &
             call depository%remove_field("grad_y_orog")
           end if
           if (ancil_option == ancil_option_fixed) then
-            if (sst_source /= sst_source_start_dump) then
-              call ancil_fields%remove_field("tstar_sea")
-              call depository%remove_field("tstar_sea")
-            end if
             if (.not. l_esm_couple) then
+              if (sst_source /= sst_source_start_dump) then
+                call ancil_fields%remove_field("tstar_sea")
+                call depository%remove_field("tstar_sea")
+              end if
               call ancil_fields%remove_field("sea_ice_fraction")
               call depository%remove_field("sea_ice_fraction")
             end if