Lfric2lfric vertical regridding infrastructure

applications/lfric2lfric/example/configuration.nml

@@ -28,6 +28,15 @@ stretching_height       = 17507.0,
 stretching_method       = 'smooth',
 /
 
+&extrusion_dst
+domain_height           = 80000.0,
+method                  = 'um_L70_50t_20s_80km',
+number_of_layers        = 70,
+planet_radius           = 6371229.0,
+stretching_height       = 17507.0,
+stretching_method       = 'smooth',
+/
+
 &files
 diag_stem_name          = 'lfric2lfric_diag',
 start_dump_filename     = 'mesh_C12_MG_dump',
@@ -97,3 +106,7 @@ coord_space             = 'Wchi',
 &planet
   scaling_factor        = 125.0,
 /
+
+&orography
+  orog_init_option      = 'none'
+/

applications/lfric2lfric/rose-meta/lfric-lfric2lfric/HEAD/rose-meta.conf

@@ -7,6 +7,121 @@ compulsory=false
 compulsory=false
 !!trigger=namelist:partitioning: .false.;
 
+[namelist:extrusion_dst]
+compulsory=true
+description=Vertical extrusion settings for the destination mesh.
+help=
+ns=namelist/Model/Mesh/Extrusion
+sort-key=Section-A03
+
+[namelist:extrusion_dst=domain_height]
+compulsory=true
+description=Height of domain above the surface [m]
+fail-if=this < 0.0 ;
+help=Height of the model domain above the surface.
+    =Surface is 0m for planar domains and planet_radius for spherical.
+!kind=default
+range=this > 0.0:
+sort-key=Panel-A02
+type=real
+
+[namelist:extrusion_dst=eta_values]
+!bounds=namelist:extrusion_dst=number_of_layers-1
+compulsory=true
+description=Eta values, excluding 0 and 1
+help=Non-dimensional height values defining the model vertical levels.
+    =The length should be smaller by number_of_layers by one.
+    =
+    =The values must strictly be between 0 and 1:
+    =0 < eta_values < 1
+    =because 0 and 1 are appended to the list automatically.
+!kind=default
+length=:
+range=0.0:1.0
+sort-key=Panel-A04
+type=real
+
+[namelist:extrusion_dst=method]
+compulsory=true
+description=Method for generating eta coordinate
+!enumeration=true
+fail-if=this == 'um_L38_29t_9s_40km'  and namelist:extrusion_dst=number_of_layers != 38 ;
+       =this == 'um_L70_50t_20s_80km' and namelist:extrusion_dst=number_of_layers != 70 ;
+       =this == 'um_L85_50t_35s_85km' and namelist:extrusion_dst=number_of_layers != 85 ;
+       =this == 'um_L70_61t_9s_40km'  and namelist:extrusion_dst=number_of_layers != 70 ;
+       =this == 'um_L120_99t_21s_40km' and namelist:extrusion_dst=number_of_layers != 120 ;
+       =this == 'um_L140_122t_18s_40km' and namelist:extrusion_dst=number_of_layers != 140 ;
+help=Available extrusion methods are (\f$n$ is number of layers):
+    =0) Specified by user;
+    =1) Uniform eta spacing (\f$\frac{k}{n}\f$);
+    =2) Quadratic eta spacing (\f$\frac{k}{n}^2\f$);
+    =3) Geometric eta spacing (\f$d\eta = \frac{(s - 1)}{(s^{n} - 1)}$)
+    =    with stretching factor prescribed (\f$s=1.03$);
+    =4) DCMIP eta spacing (Ullrich et al. (2012) DCMIP documentation, Appendix F.2.)
+    =    with flattening parameter prescribed.
+    =5) L38 40km UM specific eta spacing;
+    =6) L70 80km UM specific eta spacing;
+    =7) L85 85km UM specific eta spacing;
+    =8) L70 40km UM specific eta spacing
+    =9) L120 40km UM specific eta spacing
+    =10) L140 40km UM specific eta spacing
+sort-key=Panel-A01
+trigger=namelist:extrusion_dst=eta_values:    'specified_values';
+value-titles=Uniform, Quadratic, Geometric, DCMIP,
+            = um_L38_29t_9s_40km, um_L70_50t_20s_80km, um_L85_50t_35s_85km, um_L70_61t_9s_40km,
+            = um_L120_99t_21s_40km, um_L140_122t_18s_40km
+values='uniform', 'quadratic', 'geometric', 'dcmip',
+      ='um_L38_29t_9s_40km', 'um_L70_50t_20s_80km', 'um_L85_50t_35s_85km', 'um_L70_61t_9s_40km',
+      ='um_L120_99t_21s_40km', 'um_L140_122t_18s_40km'
+
+[namelist:extrusion_dst=number_of_layers]
+compulsory=true
+description=Number of layers in the vertical
+fail-if=this < 1 ;
+help=Setting for number of layers of 3D-mesh in vertical.
+!kind=default
+range=1:
+sort-key=Panel-A03
+type=integer
+
+[namelist:extrusion_dst=planet_radius]
+compulsory=true
+description=Radius of the planet surface [m]
+fail-if=this <= 0.0 ;
+help=Radius of domain bottom. Note: Orography not included.
+!kind=default
+range=this > 0.0:
+sort-key=Panel-A02
+type=real
+
+[namelist:extrusion_dst=stretching_height]
+compulsory=false
+description=Physical height above which surface altitude does not
+           = influence layer height
+fail-if=this < 0.0 ;
+help=Physical height above which surface altitude does not
+    =influence layer height.
+    =Note that in order to reproduce the vertical stretching used in
+    =the UM 'smooth' extrusion, this value should be set to
+    =the value of eta_rho corresponding to the UM level
+    =first_constant_r_rho_level and multiplied by domain_height.
+!kind=default
+range=0:
+sort-key=Panel-A02 ;1
+type=real
+
+[namelist:extrusion_dst=stretching_method]
+compulsory=false
+description=Method of generating stretching
+!enumeration=true
+help=Available stretching methods are:
+    =1) Linear (linear multiple of physical depth)
+    =2) Smooth (quadratic below stretching_height, linear above)
+sort-key=Panel-A01
+trigger=namelist:extrusion_dst=stretching_height: 'smooth' ;
+value-titles=Linear, Smooth
+values='linear', 'smooth'
+
 [namelist:lfric2lfric]
 compulsory=true
 description=
@@ -28,7 +143,6 @@ values='planar', 'spherical'
 [namelist:lfric2lfric=destination_mesh_name]
 compulsory=true
 description=Tag-name for destination-mesh
-fail-if=this == namelist:lfric2lfric=source_mesh_name ;
 help=Mesh topologies are held in UGRID conformant NetCDF files which
     =may contain more than one mesh topology. This tag-name identifies
     =the mesh topology to use from the mesh file namelist:lfric2lfric=destination_meshfile_prefix.
@@ -183,7 +297,6 @@ values='planar', 'spherical'
 [namelist:lfric2lfric=source_mesh_name]
 compulsory=true
 description=Tag-name for source-mesh
-fail-if=this == namelist:lfric2lfric=destination_mesh_name ;
 help=Mesh topologies are held in UGRID conformant NetCDF files which
     =may contain more than one mesh topology. This tag-name identifies
     =the mesh topology to use from the mesh file namelist:lfric2lfric=source_meshfile_prefix.

applications/lfric2lfric/source/algorithm/lfric2lfric_vert_mod.x90

@@ -0,0 +1,94 @@
+!-----------------------------------------------------------------------------
+! (C) Crown copyright 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 Apply vertical regridding using the vertical regridding kernels
+!!        and height values for the source and target vertically extruded
+!!        meshes.
+module lfric2lfric_vert_mod
+
+  use constants_mod,            only: str_def, i_def
+  use driver_modeldb_mod,       only: modeldb_type
+  use field_parent_mod,         only: field_parent_type
+  use field_mod,                only: field_type
+  use field_collection_mod,     only: field_collection_type
+  use field_collection_iterator_mod, &
+                                only: field_collection_iterator_type
+  use sci_geometric_constants_mod, &
+                                only: get_height_fv
+  use mesh_mod,                 only: mesh_type
+  use log_mod,                  only: log_event, &
+                                      log_level_info, &
+                                      log_scratch_space
+
+  implicit none
+
+  private
+  public :: lfric2lfric_vert
+
+contains
+
+  !> @brief Apply vertical regridding
+  !> @details Loop over the fields in the source_fields collection,
+  !!        and apply the vertical regridding kernels to obtain the regridded
+  !!        field in the target_fields collection.
+  !> @param modeldb        Model state object
+  !> @param source_fields  The collection of fields to be regridded
+  !> @param target_fields  The collection of regridded fields
+  subroutine lfric2lfric_vert( modeldb, source_fields, target_fields )
+
+    implicit none
+
+    type(modeldb_type),            intent(inout) :: modeldb
+    type(field_collection_type),   intent(inout) :: source_fields
+    type(field_collection_type),   intent(inout) :: target_fields
+
+    type(field_collection_iterator_type) :: iter
+
+    class(field_parent_type), pointer :: field
+    type(field_type),         pointer :: field_src
+    type(field_type),         pointer :: field_dst
+    type(field_type),         pointer :: height_src
+    type(field_type),         pointer :: height_dst
+    type(mesh_type),          pointer :: mesh_src
+    type(mesh_type),          pointer :: mesh_dst
+
+    character(len=str_def) :: field_name
+    integer(kind=i_def)    :: fs_src, fs_dst
+
+    ! Main loop over fields to be processed
+    call iter%initialise(source_fields)
+    do
+      ! Locate the field to be processed in the field collections
+      if ( .not.iter%has_next() ) exit
+      field => iter%next()
+      field_name = field%get_name()
+
+      call source_fields%get_field(field_name, field_src)
+      call target_fields%get_field(field_name, field_dst)
+
+      write(log_scratch_space, '(A,A)') "Processing lfric field ", &
+                                        trim(field_name)
+      call log_event(log_scratch_space, log_level_info)
+
+      mesh_src => field_src%get_mesh()
+      mesh_dst => field_dst%get_mesh()
+      fs_src = field_src%which_function_space()
+      fs_dst = field_dst%which_function_space()
+
+      ! Get the associated heights, ready for vertical interpolation kernels
+      height_src => get_height_fv(modeldb%config, mesh_src, fs_src)
+      height_dst => get_height_fv(modeldb%config, mesh_dst, fs_dst)
+
+      ! Dummy initialisation - just writing out the heights
+      ! TODO Replaced the dummy initialisation with actual
+      ! vertical regridding from #253. i.e.
+      ! call invoke( lfric2lfric_vert_lin_interp_lin_extrap_code_r_single( &
+      !              field_dst, field_src, height_dst, height_src))
+      call invoke(setval_X(field_dst, height_dst))
+   enddo
+
+  end subroutine lfric2lfric_vert
+
+end module lfric2lfric_vert_mod

applications/lfric2lfric/source/driver/lfric2lfric_driver_mod.F90

@@ -8,7 +8,7 @@
 module lfric2lfric_driver_mod
 
   use constants_mod,            only: str_def, str_max_filename, &
-                                      i_def, l_def, r_second
+                                      i_def, r_def, l_def, r_second
   use driver_fem_mod,           only: final_fem
   use driver_io_mod,            only: final_io
   use driver_modeldb_mod,       only: modeldb_type
@@ -24,6 +24,9 @@ module lfric2lfric_driver_mod
   use sci_checksum_alg_mod,     only: checksum_alg
   use xios,                     only: xios_date, xios_get_current_date, &
                                       xios_date_convert_to_string
+  use orography_config_mod,     only: orog_init_option_analytic, &
+                                      orog_init_option_ancil,    &
+                                      orog_init_option_start_dump
 
   !------------------------------------
   ! lfric2lfric modules
@@ -32,8 +35,11 @@ module lfric2lfric_driver_mod
   use lfric2lfric_infrastructure_mod, only: initialise_infrastructure, &
                                             context_dst, context_src,  &
                                             source_collection_name,    &
-                                            target_collection_name
+                                            target_collection_name,    &
+                                            interm_collection_name,    &
+                                            src, dst
   use lfric2lfric_regrid_mod,         only: lfric2lfric_regrid
+  use lfric2lfric_vert_mod,           only: lfric2lfric_vert
 
   implicit none
 
@@ -98,10 +104,41 @@ subroutine run( modeldb, oasis_clock )
 
     type(field_collection_type),   pointer :: source_fields
     type(field_collection_type),   pointer :: target_fields
+    type(field_collection_type),   pointer :: interm_fields
 
     type(lfric_xios_context_type), pointer :: io_context
 
+    character(len=str_def)  :: mesh_names(2)
+    integer(kind=i_def)     :: src_extrusion_method
+    integer(kind=i_def)     :: dst_extrusion_method
+    integer(kind=i_def)     :: src_number_of_layers
+    integer(kind=i_def)     :: dst_number_of_layers
+    integer(kind=i_def)     :: src_stretching_method
+    integer(kind=i_def)     :: dst_stretching_method
+    integer(kind=i_def)     :: orog_init_option
+    real(kind=r_def)        :: src_domain_height
+    real(kind=r_def)        :: dst_domain_height
+    real(kind=r_def)        :: src_stretching_height
+    real(kind=r_def)        :: dst_stretching_height
+
+    logical(kind=l_def), pointer :: vertical_change
+    logical(kind=l_def), pointer :: horizontal_change
+
     ! Extract configuration variables
+    src_extrusion_method    = modeldb%config%extrusion%method()
+    src_number_of_layers    = modeldb%config%extrusion%number_of_layers()
+    src_domain_height       = modeldb%config%extrusion%domain_height()
+    src_stretching_height   = modeldb%config%extrusion%stretching_height()
+    src_stretching_method   = modeldb%config%extrusion%stretching_method()
+    dst_extrusion_method    = modeldb%config%extrusion_dst%method()
+    dst_number_of_layers    = modeldb%config%extrusion_dst%number_of_layers()
+    dst_domain_height       = modeldb%config%extrusion_dst%domain_height()
+    dst_stretching_height   = modeldb%config%extrusion_dst%stretching_height()
+    dst_stretching_method   = modeldb%config%extrusion_dst%stretching_method()
+    orog_init_option        = modeldb%config%orography%orog_init_option()
+
+    mesh_names(dst)      = modeldb%config%lfric2lfric%destination_mesh_name()
+    mesh_names(src)      = modeldb%config%lfric2lfric%source_mesh_name()
     start_dump_filename  = modeldb%config%files%start_dump_filename()
     checkpoint_stem_name = modeldb%config%files%checkpoint_stem_name()
 
@@ -112,12 +149,30 @@ subroutine run( modeldb, oasis_clock )
     source_fields => modeldb%fields%get_field_collection(source_collection_name)
     target_fields => modeldb%fields%get_field_collection(target_collection_name)
 
+    call modeldb%values%get_value("vertical_change", vertical_change)
+    call modeldb%values%get_value("horizontal_change", horizontal_change)
+
+    if (horizontal_change .and. vertical_change) then
+       interm_fields => modeldb%fields%get_field_collection(interm_collection_name)
+
+    else if (horizontal_change .and. .not. vertical_change) then
+       interm_fields =>  modeldb%fields%get_field_collection(target_collection_name)
+
+    else if (vertical_change .and. .not. horizontal_change) then
+       interm_fields =>  modeldb%fields%get_field_collection(source_collection_name)
+    end if
+
     ! Read fields and perform the regridding
     if (mode == mode_ics) then
       call read_state(source_fields, prefix='restart_')
 
-      call lfric2lfric_regrid(modeldb, oasis_clock, source_fields,   &
-                              target_fields, regrid_method)
+      if (horizontal_change) then
+        call lfric2lfric_regrid(modeldb, oasis_clock, source_fields,   &
+                                interm_fields, regrid_method)
+      end if
+      if (vertical_change) then
+        call lfric2lfric_vert(modeldb, interm_fields, target_fields)
+      end if
 
       ! Write output
       call modeldb%io_contexts%get_io_context(context_dst, io_context)
@@ -139,8 +194,13 @@ subroutine run( modeldb, oasis_clock )
 
         call read_state(source_fields)
 
-        call lfric2lfric_regrid(modeldb, oasis_clock, source_fields, &
-                                target_fields, regrid_method)
+        if (horizontal_change) then
+          call lfric2lfric_regrid(modeldb, oasis_clock, source_fields,   &
+                                  interm_fields, regrid_method)
+        end if
+        if (vertical_change) then
+          call lfric2lfric_vert(modeldb, interm_fields, target_fields)
+        end if
 
         is_running = modeldb%clock%tick()