S20 Diags: icao heights

CONTRIBUTORS.md

@@ -1,13 +1,14 @@
 # Contributors
 | GitHub user     | Real Name         | Affiliation | Date       |
-| --------------- | ----------------- | ----------- | ---------- |
+|-----------------|-------------------| ----------- |------------|
 | james-bruten-mo | James Bruten      | Met Office  | 2025-12-09 |
 | jedbakerMO      | Jed Baker         | Met Office  | 2025-12-29 |
 | jennyhickson    | Jenny Hickson     | Met Office  | 2025-12-10 |
 | mike-hobson     | Mike Hobson       | Met Office  | 2025-12-17 |
 | mo-marqh        | mark Hedley       | Met Office  | 2025-12-11 |
 | yaswant         | Yaswant Pradhan   | Met Office  | 2025-12-16 |
 | oakleybrunt     | Oakley Brunt      | Met Office  | 2025-12-19 |
+| bblay-mo        | Byron Blay        | Met Office  | 2026-01-07 |
 | harry-shepherd  | Harry Shepherd    | Met Office  | 2026-01-08 |
 | DrTVockerodtMO  | Terence Vockerodt | Met Office  | 2026-01-08 |
 | MetBenjaminWent | Benjamin Went     | Met Office  | 2026-01-15 |

applications/lfric_atm/example/iodef.xml

@@ -94,6 +94,13 @@
 	<field field_ref="forcing__dt_force" />
       </file>
 
+	<file id="lf_diag_aviation" name="lf_diag_aviation" output_freq="6h" convention="UGRID" enabled=".TRUE.">
+		<field field_ref="aviation__geopot_thickness_850"/>
+		<field field_ref="aviation__geopot_thickness_500"/>
+		<field field_ref="aviation__conv_cloud_base_icao_height"/>
+		<field field_ref="aviation__conv_cloud_top_icao_height"/>
+	</file>
+
       <file id="lfric_averages" name="lfric_averages" output_freq="12h" convention="UGRID" enabled=".TRUE.">
         <field field_ref="theta" operation="average" />
         <field field_ref="exner" operation="average" />

applications/lfric_atm/metadata/field_def_diags.xml

@@ -989,6 +989,11 @@
   <!-- Vertical integrals -->
   <field id="processed__sw_aer_optical_depth_rts" field_ref="radiation__sw_aer_optical_depth_rts" grid_ref="vert_sum_face_grid"/>
   <field id="processed__lw_aer_optical_depth_rts" field_ref="radiation__lw_aer_optical_depth_rts" grid_ref="vert_sum_face_grid"/>
+  <!-- Aviation diagnostics -->
+  <field id="aviation__geopot_thickness_850" name="aviation_geopot_thickness_850" standard_name="atmosphere_layer_thickness_expressed_as_geopotential_height_difference" long_name="geopotential_height_thickness_between_1000hPa_and_850hPa" unit="m" domain_ref="face"/>
+  <field id="aviation__geopot_thickness_500" name="aviation_geopot_thickness_500" standard_name="atmosphere_layer_thickness_expressed_as_geopotential_height_difference" long_name="geopotential_height_thickness_between_1000hPa_and_500hPa" unit="m" domain_ref="face"/>
+  <field id="aviation__conv_cloud_base_icao_height" name="aviation_conv_cloud_base_icao_height" long_name="height of convective cloud base" unit="kft" domain_ref="face" />
+  <field id="aviation__conv_cloud_top_icao_height" name="aviation_conv_cloud_top_icao_height" long_name="height of convective cloud top" unit="kft" domain_ref="face" />
 </field_group>
 
 <!-- Integer diagnostic group -->

interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90

@@ -0,0 +1,253 @@
+!-----------------------------------------------------------------------------
+! (C) Crown copyright 2025 Met Office. All rights reserved.
+! The file LICENCE, distributed with this code, contains details of the terms
+! under which the code may be used.
+!-----------------------------------------------------------------------------
+!
+! The main entry point for calculating section 20 aviation diagnostics.
+!
+! Code Owner: Please refer to the UM file CodeOwners.txt
+! This file currently belongs in section: physics_schemes_interface
+! whilst discussions are ongoing about its final location.
+!
+MODULE aviation_diags_alg_mod
+
+  USE constants_mod,              ONLY: r_def, i_def, l_def, RMDI
+  USE field_mod,                  ONLY: field_type
+  USE field_collection_mod,       ONLY: field_collection_type
+  USE io_config_mod,              ONLY: USE_xios_io
+  USE initialise_diagnostics_mod, ONLY: init_diag => init_diagnostic_field
+
+  USE lfric_xios_diag_mod,        ONLY: get_axis_dimension, get_axis_values
+  USE log_mod,                    ONLY: log_event, log_scratch_space, &
+                                        LOG_LEVEL_ALWAYS, LOG_LEVEL_ERROR, &
+                                        LOG_LEVEL_INFO, LOG_LEVEL_DEBUG
+
+  USE diags_geopot_kernel_mod,  ONLY: diags_geopot_kernel_type
+  USE diags_icao_heights_kernel_mod,  ONLY: diags_icao_heights_kernel_type
+
+  IMPLICIT NONE
+
+PRIVATE
+PUBLIC :: aviation_diags_alg
+
+CONTAINS
+
+  ! Algorithm to calculate the section 20 aviation diagnostics.
+  SUBROUTINE aviation_diags_alg(plev_geopot, convection_fields)
+
+    IMPLICIT NONE
+
+    ! Arguments
+    TYPE(field_type), INTENT(IN) :: plev_geopot
+    TYPE(field_collection_type), INTENT(INOUT) :: convection_fields
+
+
+    ! Calculate thickness from geopotential height at pressure levels.
+    CALL from_plev_geopot(plev_geopot)
+
+
+    ! Calculate icao heights of top and base of cloud from their pressures.
+    CALL icao_heights(convection_fields)
+
+
+  END SUBROUTINE aviation_diags_alg
+
+  SUBROUTINE from_plev_geopot(plev_geopot)
+    ! Calculate diagnostics from plvel_geopot.
+    ! - Thickness from geopotential height at pressure levels.
+
+    ! Arguments
+    TYPE(field_type), INTENT(IN) :: plev_geopot
+
+    ! Local variables
+
+    ! These flags tell us which results are requested at this time step.
+    LOGICAL(l_def) :: aviation_thick_850_flag, aviation_thick_500_flag
+
+    ! The array of pressure levels.
+    INTEGER(I_DEF) :: nplev
+    REAL(R_DEF), ALLOCATABLE :: plevs(:)
+    INTEGER(I_DEF) :: plevs_alloc_stat
+
+    ! Level indices for the three pressure levels we want to read.
+    INTEGER(I_DEF) :: i1000, i850, i500
+
+    ! Approximate equality tolerance - is there a common variable somewhere?
+    REAL(R_DEF), PARAMETER :: plev_tol = 0.1_r_def
+
+    ! The two output fields we produce.
+    TYPE( field_type ) :: aviation_thick_850, aviation_thick_500
+
+    INTEGER(I_DEF) :: i
+
+
+    ! Check the request flags.
+    aviation_thick_850_flag = &
+      init_diag(aviation_thick_850, 'aviation__geopot_thickness_850')
+    aviation_thick_500_flag = &
+      init_diag(aviation_thick_500, 'aviation__geopot_thickness_500')
+    IF ( aviation_thick_850_flag .OR. aviation_thick_500_flag ) THEN
+      WRITE(log_scratch_space, '(A)') 'Section 20 thickness is on'
+      CALL log_event(log_scratch_space, LOG_LEVEL_DEBUG)
+    ELSE
+      ! Nothing requested at this time step.
+      RETURN
+    END IF
+
+    ! Get the array of pressure levels.
+    nplev = get_axis_dimension('pressure_levels')
+    IF (nplev <= 0) THEN
+      WRITE(log_scratch_space, '(A, I0)') 'No pressure levels, nplev=', nplev
+      CALL log_event(log_scratch_space, LOG_LEVEL_ERROR)
+      RETURN
+    END IF
+
+    ALLOCATE(plevs(nplev), stat=plevs_alloc_stat)
+    IF (plevs_alloc_stat /= 0) THEN
+      WRITE(log_scratch_space, '(A, I0)') 'allocate(plevs) failed, stat=', &
+        plevs_alloc_stat
+      CALL log_event(log_scratch_space, LOG_LEVEL_ERROR)
+      RETURN
+    END IF
+    plevs = get_axis_values('pressure_levels',nplev)
+
+    ! Find the level indices for our three pressures of interest.
+    ! Leave as -1 if not requested, telling the kernel not to calculate it.
+    ! Assumes hPa. Should we worry about other units?
+    i1000 = -1
+    i850  = -1
+    i500  = -1
+    DO i = 1, nplev
+
+      ! 1000?
+      IF ( abs(plevs(i) - 100000.0_r_def) < plev_tol ) THEN
+        i1000 = i
+
+      ! 850?
+      ELSE IF ( abs(plevs(i) - 85000.0_r_def) < plev_tol ) THEN
+        i850 = i
+
+      ! 500?
+      ELSE IF ( abs(plevs(i) - 50000.0_r_def) < plev_tol ) THEN
+        i500 = i
+
+      END IF
+    END DO
+
+    ! Check we found the required levels.
+    IF (i1000 == -1) THEN
+        WRITE(log_scratch_space, '(A)') 'could not find 1000hPa'
+        CALL log_event(log_scratch_space, LOG_LEVEL_ERROR)
+        RETURN
+    END IF
+
+    IF (aviation_thick_850_flag .AND. i850 == -1) THEN
+      WRITE(log_scratch_space, '(A)') 'could not find 850hPa'
+      CALL log_event(log_scratch_space, LOG_LEVEL_ERROR)
+      RETURN
+    END IF
+
+    IF (aviation_thick_500_flag .AND. i500 == -1) THEN
+      WRITE(log_scratch_space, '(A)') 'could not find 500hPa'
+      CALL log_event(log_scratch_space, LOG_LEVEL_ERROR)
+      RETURN
+    END IF
+
+    ! Call the kernel to subtract the levels.
+    ! Todo: Comment why we initialise to 1.0 and not 0.0. I can't remember.
+    CALL invoke( &
+      setval_c(aviation_thick_850, 1.0_r_def), &
+      setval_c(aviation_thick_500, 1.0_r_def), &
+      diags_geopot_kernel_type( &
+        aviation_thick_850, aviation_thick_500, &
+        plev_geopot, &
+        aviation_thick_850_flag, aviation_thick_500_flag, &
+        i1000, i850, i500) &
+    )
+
+    ! Write the fields
+    IF (aviation_thick_850_flag) CALL aviation_thick_850%write_field()
+    IF (aviation_thick_500_flag) CALL aviation_thick_500%write_field()
+
+    ! Clean up.
+    IF (allocated(plevs)) THEN
+      DEALLOCATE(plevs)
+    END IF
+
+  END SUBROUTINE from_plev_geopot
+
+
+  SUBROUTINE icao_heights(convection_fields)
+    ! Calculate iaco heights of top and base of cloud from their pressures.
+
+    ! Arguments
+    TYPE(field_collection_type), INTENT(INOUT) :: convection_fields
+
+    ! Locals
+
+    ! Pressure fields to unpack.
+    type(field_type), pointer :: conv_cloud_base_pressure
+    type(field_type), pointer :: conv_cloud_top_pressure
+
+    ! Height fields to create and write.
+    TYPE(field_type) :: conv_cloud_base_icao_height
+    TYPE(field_type) :: conv_cloud_top_icao_height
+
+    ! Request flags.
+    LOGICAL(l_def) :: base_height_flag
+    LOGICAL(l_def) :: top_height_flag
+
+    ! which heights are requested this time step?
+    top_height_flag = init_diag( &
+      conv_cloud_top_icao_height, 'aviation__conv_cloud_top_icao_height')
+
+    base_height_flag = init_diag( &
+      conv_cloud_base_icao_height, 'aviation__conv_cloud_base_icao_height')
+
+    if (base_height_flag .or. top_height_flag) then
+      WRITE(log_scratch_space, '(A)') 'Section 20 cloud ICAO heights on'
+      CALL log_event(log_scratch_space, LOG_LEVEL_DEBUG)
+    else
+      ! nothing requested at this time step
+      return
+    end if
+
+
+    ! cloud base
+    IF (base_height_flag) THEN
+      call convection_fields%get_field( &
+              'convection__pres_cv_base', conv_cloud_base_pressure)
+      call invoke( &
+        setval_c(conv_cloud_base_icao_height, RMDI), &
+        diags_icao_heights_kernel_type(conv_cloud_base_icao_height, &
+          conv_cloud_base_pressure))
+
+      call conv_cloud_base_icao_height%write_field()
+
+      ! Is it necessary/safe to remove this field now?
+      ! Did the get trigger an unpack which we need to now remove?
+      ! Based on https://code.metoffice.gov.uk/trac/lfric_apps/ticket/189
+      ! call convection_fields%remove_field('convection__pres_cv_base')
+
+    END IF
+
+    ! could top
+    IF (top_height_flag) THEN
+      call convection_fields%get_field( &
+              'convection__pres_cv_top', conv_cloud_top_pressure)
+      call invoke( &
+        setval_c(conv_cloud_top_icao_height, RMDI), &
+        diags_icao_heights_kernel_type(conv_cloud_top_icao_height, &
+          conv_cloud_top_pressure))
+
+      call conv_cloud_top_icao_height%write_field()
+
+      ! Is it necessary/safe to remove this field now?
+      ! call convection_fields%remove_field('convection__pres_cv_top')
+
+    END IF
+
+  END SUBROUTINE icao_heights
+
+END MODULE aviation_diags_alg_mod

interfaces/physics_schemes_interface/source/algorithm/conv_gr_alg_mod.x90

@@ -27,6 +27,8 @@ module conv_gr_alg_mod
   use log_mod,              only: log_event, LOG_LEVEL_DEBUG
   use mesh_mod,             only: mesh_type
 
+  use initialise_diagnostics_mod, only : diag_samp => diagnostic_to_be_sampled
+
   implicit none
 
   private
@@ -62,14 +64,14 @@ contains
     type( field_type ), intent( inout )     :: mr(nummr)
     type( field_type ), intent( in )        :: rho, exner
 
-    type( field_collection_type ), intent(in) :: derived_fields
-    type( field_collection_type ), intent(in) :: turbulence_fields
-    type( field_collection_type ), intent(in) :: convection_fields
-    type( field_collection_type ), intent(in) :: cloud_fields
-    type( field_collection_type ), intent(in) :: surface_fields
-    type( field_collection_type ), intent(in) :: chemistry_fields
-    type( field_collection_type ), intent(in) :: aerosol_fields
-    type( field_collection_type ), intent(in) :: microphysics_fields
+    type( field_collection_type ), intent(in)    :: derived_fields
+    type( field_collection_type ), intent(in)    :: turbulence_fields
+    type( field_collection_type ), intent(inout) :: convection_fields
+    type( field_collection_type ), intent(in)    :: cloud_fields
+    type( field_collection_type ), intent(in)    :: surface_fields
+    type( field_collection_type ), intent(in)    :: chemistry_fields
+    type( field_collection_type ), intent(in)    :: aerosol_fields
+    type( field_collection_type ), intent(in)    :: microphysics_fields
 
     integer(kind=i_def), intent(in)         :: outer
 
@@ -767,6 +769,20 @@ contains
                                   dth_conv_noshal,       &
                                   dmv_conv_noshal)
 
+
+      ! Aviation dependencies to activate pres_cv_base and pres_cv_top
+      if (diag_samp('aviation__conv_cloud_base_icao_height')) then
+        if (.not. convection_fields%field_exists('convection__pres_cv_base')) then
+          call convection_fields%add_field(pres_cv_base)
+        end if
+      end if
+
+      if (diag_samp('aviation__conv_cloud_top_icao_height')) then
+        if (.not. convection_fields%field_exists('convection__pres_cv_top')) then
+          call convection_fields%add_field(pres_cv_top)
+        end if
+      end if
+
     end if
 
     nullify( mesh )