Iss615/read ancils in time windows

interfaces/jules_interface/source/algorithm/init_surface_fields_alg_mod.x90

@@ -15,7 +15,9 @@ module init_surface_fields_alg_mod
   use c_kappai,                  only: kappai, de
 
   use initialization_config_mod, only: init_option, init_option_analytic, &
-                                       ancil_option, ancil_option_none
+                                       ancil_option, ancil_option_none,   &
+                                       ancil_option_updating
+  use jules_radiation_config_mod, only: l_sea_alb_var_chl
   use specified_surface_config_mod, only: internal_flux_method, &
                                        internal_flux_method_uniform, &
                                        internal_flux_value, &
@@ -129,9 +131,14 @@ contains
     call surface_fields%get_field('ocn_cpl_point',ocn_cpl_point)
     call surface_fields%get_field('wspd10m', wspd10m)
 
+    ! In updating mode the sea chlorophyll is set by the ancillary file object
+    if ( .not. ( ancil_option == ancil_option_updating .and. &
+                 l_sea_alb_var_chl ) ) then
+      call invoke( setval_c(chloro_sea, 5.0e-7_r_def) )
+    end if
+
     call invoke(                                                              &
                  ! Need initialising for aquaplanet
-                 setval_c(chloro_sea,          5.0e-7_r_def),                 &
                  setval_c(sea_ice_thickness,   1.0_r_def),                    &
                  setval_c(sea_ice_temperature, 260.0_r_def),                  &
                  ! Diagnostics

interfaces/jules_interface/source/algorithm/update_ancils_alg_mod.x90

@@ -76,7 +76,7 @@ contains
       if (.not. l_couple_sea_ice) then
         call ancil_fields%get_field('sea_ice_fraction', sea_ice_fraction)
       end if
-      if (.not. l_couple_ocean) then
+      if (.not. l_couple_ocean .and. ancil_fields%field_exists('tstar_sea')) then
         call ancil_fields%get_field('tstar_sea', tstar_sea)
       end if
 

interfaces/jules_interface/source/kernel/initial_surface_kernel_mod.F90

@@ -21,14 +21,15 @@ module initial_surface_kernel_mod
 
   use jules_control_init_mod, only: n_land_tile, &
                                     n_sea_tile, n_sea_ice_tile, &
-       first_sea_tile, first_sea_ice_tile, n_surf_tile
+       first_sea_tile, first_sea_ice_tile, n_surf_tile, pft_ancil_file_owned
   use ideal_surface_config_mod, only: surf_tile_fracs, surf_tile_temps, &
                                       canopy_height_in => canopy_height, leaf_area_index_in => leaf_area_index
   use jules_surface_types_mod, only: npft
   use nvegparm, only: emis_nvg
   use pftparm, only: emis_pft
   use jules_sea_seaice_mod, only: emis_sea, emis_sice
   use section_choice_config_mod, only: surface, surface_jules
+  use initialization_config_mod, only: ancil_option, ancil_option_updating
 
   implicit none
 
@@ -100,11 +101,14 @@ subroutine initial_surface_code(nlayers,                       &
 
     if (surface == surface_jules) then
 
-      ! Set from a namelist variable when no ancil read
-      do i = 1, npft
-        leaf_area_index(map_pft(1)+i-1) = leaf_area_index_in(i)
-        canopy_height(map_pft(1)+i-1) = canopy_height_in(i)
-      end do
+      ! Set from a namelist variable when no ancil read. When the ancillary
+      ! file object owns these fields they have already been set.
+      if (.not. (ancil_option == ancil_option_updating .and. pft_ancil_file_owned)) then
+        do i = 1, npft
+          leaf_area_index(map_pft(1)+i-1) = leaf_area_index_in(i)
+          canopy_height(map_pft(1)+i-1) = canopy_height_in(i)
+        end do
+      end if
 
       ! Initialise grey_albedos for all tiles
       do i = 1, n_land_tile

interfaces/jules_interface/source/support/jules_control_init_mod.f90

@@ -12,7 +12,7 @@ module jules_control_init_mod
   use jules_sea_seaice_config_mod,          only : nice_in => nice
 
   ! Other LFRic modules used
-  use constants_mod,        only : i_def, imdi
+  use constants_mod,        only : i_def, l_def, imdi
 
   implicit none
 
@@ -23,6 +23,11 @@ module jules_control_init_mod
   ! Has nothing to do with number of surface types
   integer(kind=i_def), parameter :: n_surf_interp = 11
 
+  ! Temporary development switch. When true and ancil_option is 'updating',
+  ! leaf_area_index and canopy_height are read by an ancillary file object
+  ! rather than the legacy time axis.
+  logical(kind=l_def), parameter :: pft_ancil_file_owned = .false.
+
   integer(kind=i_def), protected :: n_land_tile
   integer(kind=i_def), protected :: n_sea_ice_tile
   integer(kind=i_def), protected :: n_surf_tile
@@ -33,7 +38,7 @@ module jules_control_init_mod
   private
   public :: n_land_tile, n_sea_tile, n_sea_ice_tile, n_surf_tile, &
        first_sea_tile, first_sea_ice_tile, jules_control_init,    &
-       soil_lev_tile, n_surf_interp
+       soil_lev_tile, n_surf_interp, pft_ancil_file_owned
 
 contains
 

interfaces/physics_schemes_interface/source/algorithm/init_aerosol_fields_alg_mod.x90

@@ -29,6 +29,7 @@ contains
   subroutine init_aerosol_fields_alg(aerosol_fields)
 
     use section_choice_config_mod, only: aerosol, aerosol_um
+    use initialization_config_mod, only: ancil_option, ancil_option_updating
     use aerosol_config_mod, only: glomap_mode,               &
                                   glomap_mode_climatology,   &
                                   glomap_mode_dust_and_clim, &
@@ -330,24 +331,27 @@ contains
            glomap_mode == glomap_mode_dust_and_clim ) ) then
 
       ! GLOMAP-mode climatology scheme selected.
-      call invoke( setval_c(n_ait_sol,             0.0_r_def), &
-                   setval_c(ait_sol_su,            0.0_r_def), &
-                   setval_c(ait_sol_bc,            0.0_r_def), &
-                   setval_c(ait_sol_om,            0.0_r_def), &
-                   setval_c(n_acc_sol,             0.0_r_def), &
-                   setval_c(acc_sol_su,            0.0_r_def), &
-                   setval_c(acc_sol_bc,            0.0_r_def), &
-                   setval_c(acc_sol_om,            0.0_r_def), &
-                   setval_c(acc_sol_ss,            0.0_r_def), &
-                   setval_c(n_cor_sol,             0.0_r_def), &
-                   setval_c(cor_sol_su,            0.0_r_def), &
-                   setval_c(cor_sol_bc,            0.0_r_def), &
-                   setval_c(cor_sol_om,            0.0_r_def), &
-                   setval_c(cor_sol_ss,            0.0_r_def), &
-                   setval_c(n_ait_ins,             0.0_r_def), &
-                   setval_c(ait_ins_bc,            0.0_r_def), &
-                   setval_c(ait_ins_om,            0.0_r_def), &
-                   setval_c(n_acc_ins,             0.0_r_def), &
+      ! In updating mode these fields are set by the ancillary file object.
+      if (ancil_option /= ancil_option_updating) then
+        call invoke( setval_c(n_ait_sol,             0.0_r_def), &
+                     setval_c(ait_sol_su,            0.0_r_def), &
+                     setval_c(ait_sol_bc,            0.0_r_def), &
+                     setval_c(ait_sol_om,            0.0_r_def), &
+                     setval_c(n_acc_sol,             0.0_r_def), &
+                     setval_c(acc_sol_su,            0.0_r_def), &
+                     setval_c(acc_sol_bc,            0.0_r_def), &
+                     setval_c(acc_sol_om,            0.0_r_def), &
+                     setval_c(acc_sol_ss,            0.0_r_def), &
+                     setval_c(n_cor_sol,             0.0_r_def), &
+                     setval_c(cor_sol_su,            0.0_r_def), &
+                     setval_c(cor_sol_bc,            0.0_r_def), &
+                     setval_c(cor_sol_om,            0.0_r_def), &
+                     setval_c(cor_sol_ss,            0.0_r_def), &
+                     setval_c(n_ait_ins,             0.0_r_def), &
+                     setval_c(ait_ins_bc,            0.0_r_def), &
+                     setval_c(ait_ins_om,            0.0_r_def))
+      end if
+      call invoke( setval_c(n_acc_ins,             0.0_r_def), &
                    setval_c(acc_ins_du,            0.0_r_def), &
                    setval_c(n_cor_ins,             0.0_r_def), &
                    setval_c(cor_ins_du,            0.0_r_def), &

science/gungho/source/driver/field_maker_mod.F90

@@ -342,22 +342,30 @@ subroutine add_real_field(field_collection, &
     type(function_space_type), pointer             :: window_size_space
     logical(l_def)                                 :: checkpointed
     logical(l_def)                                 :: advected
+    logical(l_def)                                 :: preexisting
 
     ! pointers for xios write interface
     procedure(write_interface), pointer :: write_behaviour => null()
     procedure(read_interface),  pointer :: read_behaviour => null()
     procedure(checkpoint_write_interface), pointer :: checkpoint_write_behaviour => null()
     procedure(checkpoint_read_interface), pointer  :: checkpoint_read_behaviour => null()
 
+    ! Fields read by file objects are created before the I/O context closes
+    preexisting = .false.
+    if (.not. associated(external_field)) preexisting = depository%field_exists(name)
+
     ! Create the new field
     if (associated(external_field)) then
       new_field_ptr => external_field
+    else if (preexisting) then
+      call depository%get_field(name, field_ptr)
+      new_field_ptr => field_ptr
     else
       new_field_ptr => new_field
     end if
 
     ! deal with time axis
-    if (associated(time_axis)) then
+    if (associated(time_axis) .and. .not. preexisting) then
       ! pre-initialise field with window size, add to time axis;
       ! cf. init_time_axis_mod / setup_field
       if (associated(vector_space)) then
@@ -386,7 +394,9 @@ subroutine add_real_field(field_collection, &
     end if
 
     ! regular field initialisation
-    if (associated(vector_space)) then
+    if (preexisting) then
+      ! Already initialised when it was created
+    else if (associated(vector_space)) then
       if (empty) then
         call new_field_ptr%initialise( vector_space, name=trim(name), &
           override_data = empty_real_data )
@@ -434,6 +444,8 @@ subroutine add_real_field(field_collection, &
     ! Add the field to the depository, unless it is external
     if (associated(external_field)) then
       tmp_ptr => external_field
+    else if (preexisting) then
+      tmp_ptr => field_ptr
     else
       call depository%add_field(new_field)
       call depository%get_field(name, field_ptr)

science/gungho/source/driver/gungho_driver_mod.F90

@@ -76,6 +76,9 @@ module gungho_driver_mod
   use lfric_xios_time_axis_mod,    only : regridder
   use intermesh_mappings_alg_mod,  only : map_scalar_intermesh
   use update_ancils_alg_mod,       only : update_ancils_alg
+  use field_mod,                   only : field_type
+  use init_ancils_mod,             only : aerosol_ancil_names
+  use sci_field_minmax_alg_mod,    only : log_field_minmax
   use gas_calc_all_mod,            only : gas_calc_all
   use multires_coupling_config_mod, &
                                    only : lowest_order_aero_flag
@@ -384,6 +387,9 @@ subroutine step( modeldb )
 
     type( field_collection_type ), pointer :: surface_fields
     type( field_collection_type ), pointer :: ancil_fields
+    type( field_type ),            pointer :: aerosol_fld => null()
+    integer(i_def)                         :: i_aer
+    character(len=str_def)                 :: aerosol_label
 #endif
 
     integer(i_def) :: theta_forcing
@@ -526,6 +532,21 @@ subroutine step( modeldb )
                               modeldb%clock,               &
                               ancil_fields,                &
                               surface_fields)
+
+      ! Log min/max of the aerosol ancillary fields each timestep so that
+      ! runs using the file-owned read can be compared with the legacy path
+      depository => modeldb%fields%get_field_collection("depository")
+      do i_aer = 1, size(aerosol_ancil_names, 2)
+        if ( depository%field_exists( trim(aerosol_ancil_names(1,i_aer)) ) ) then
+          call depository%get_field( trim(aerosol_ancil_names(1,i_aer)), &
+                                     aerosol_fld )
+          write(aerosol_label, '(A,I0,A,A)') 'aerosol ancil step ',     &
+                modeldb%clock%get_step(), ' ',                          &
+                trim(aerosol_ancil_names(1,i_aer))
+          call log_field_minmax( LOG_LEVEL_INFO, trim(aerosol_label), aerosol_fld )
+        end if
+      end do
+      nullify(aerosol_fld)
     end if
 
     ! Update the time varying trace gases

science/gungho/source/driver/gungho_model_mod.F90

@@ -98,6 +98,9 @@ module gungho_model_mod
 
 #ifdef UM_PHYSICS
   use gas_calc_all_mod,            only : gas_calc_all
+  use init_ancils_mod,             only : create_file_owned_ancils, &
+                                          easy_aerosol_names,       &
+                                          easy_aerosol_collection_name
   use jules_control_init_mod,      only : jules_control_init
   use jules_physics_init_mod,      only : jules_physics_init
   use planet_constants_mod,        only : set_planet_constants
@@ -569,6 +572,7 @@ subroutine initialise_infrastructure( io_context_name, modeldb )
 #ifdef UM_PHYSICS
     type(field_collection_type),  pointer :: radiation_fields
 #endif
+    integer :: ancil_index
     integer :: start_index, end_index
 
     character(str_def) :: prime_mesh_name
@@ -997,6 +1001,25 @@ subroutine initialise_infrastructure( io_context_name, modeldb )
 
     call log_event("Initialising I/O context", LOG_LEVEL_INFO)
 
+    ! Collections of fields read by the ancillary file objects. They are needed
+    ! by the file list, so must exist before the context is closed.
+    call modeldb%fields%add_empty_field_collection("aerosol_ancil_fields", table_len = 40)
+    call modeldb%fields%add_empty_field_collection("albedo_vis_ancil_fields", table_len = 4)
+    call modeldb%fields%add_empty_field_collection("albedo_nir_ancil_fields", table_len = 4)
+    call modeldb%fields%add_empty_field_collection("cdnc_ancil_fields", table_len = 4)
+    call modeldb%fields%add_empty_field_collection("plant_func_ancil_fields", table_len = 4)
+    call modeldb%fields%add_empty_field_collection("sea_ancil_fields", table_len = 4)
+    call modeldb%fields%add_empty_field_collection("ozone_ancil_fields", table_len = 4)
+    do ancil_index = 1, size(easy_aerosol_names)
+      call modeldb%fields%add_empty_field_collection( &
+             easy_aerosol_collection_name(easy_aerosol_names(ancil_index)), table_len = 4 )
+    end do
+#ifdef UM_PHYSICS
+    if ( use_physics ) then
+      call create_file_owned_ancils( modeldb, mesh, twod_mesh )
+    end if
+#endif
+
     files_init_ptr => init_gungho_files
 
     if ( use_multires_coupling ) then