S20 Diags: geopot thickness

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,11 @@
 	<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"/>
+	</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,9 @@
   <!-- 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_group>
 
 <!-- Integer diagnostic group -->

interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90

@@ -0,0 +1,162 @@
+!-----------------------------------------------------------------------------
+! (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
+  USE field_mod,                  ONLY: field_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 aviation_diags_kernel_mod,  ONLY: aviation_diags_kernel_type
+
+  IMPLICIT NONE
+
+PRIVATE
+PUBLIC :: aviation_diags_alg
+
+CONTAINS
+
+  ! Algorithm to calculate the section 20 aviation diagnostics.
+  SUBROUTINE aviation_diags_alg(plev_geopot)
+
+    IMPLICIT NONE
+
+    ! Arguments
+
+    ! We calculate the result from this input field.
+    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), &
+      aviation_diags_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 aviation_diags_alg
+
+END MODULE aviation_diags_alg_mod

interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90

@@ -43,8 +43,9 @@ contains
   !> @param[in]  exner_w3        Exner pressure
   !> @param[in]  mr              Mixing ratio bundle
   !> @param[in]  moist_dyn       Moist_dynamics factors (for 1+sum(m_x))
+  !> @param[out] plev_geopot     Geopotential height on pressure levels
   subroutine pres_lev_diags_alg(derived_fields, theta_wth, exner_w3, &
-                                mr, moist_dyn)
+                                mr, moist_dyn, plev_geopot)
 
     use psykal_lite_phys_mod, only: invoke_pres_interp_kernel_type,  &
                                     invoke_heaviside_kernel_type,    &
@@ -59,6 +60,7 @@ contains
     type( field_type ), intent(in) :: theta_wth, exner_w3
     type( field_type ), intent(in) :: mr(nummr)
     type( field_type ), intent(in) :: moist_dyn(num_moist_factors)
+    type( field_type ), intent(out) :: plev_geopot
     ! Internal variables
     type( field_type ), pointer     :: exner_wth
     type( field_type ), pointer     :: u_in_w3
@@ -70,7 +72,7 @@ contains
 
     ! Define fields
     type( field_type ) :: plev_heaviside, plev_temp, plev_u, plev_v, &
-         plev_w, plev_geopot, temp, omega, plev_omega, plev_mv, plev_qv, qv, &
+         plev_w, temp, omega, plev_omega, plev_mv, plev_qv, qv, &
          plev_thetaw, plev_temp_save
     integer(i_def) :: nplev
     real(r_def) :: minus_g
@@ -81,9 +83,13 @@ contains
          plev_u_flag, plev_u_clim_flag, plev_v_flag, plev_v_clim_flag, &
          plev_w_flag, plev_w_clim_flag, plev_omega_clim_flag,          &
          plev_mv_clim_flag, plev_qv_clim_flag,                         &
-         plev_geopot_flag, plev_geopot_clim_flag, plev_thetaw_flag
+         plev_geopot_flag, plev_geopot_clim_flag, plev_thetaw_flag,    &
+         aviation_thickness_flag
     integer(tik) :: id
 
+    ! For code tidiness, combines all flags requiring plev_geopot.
+    logical(l_def) :: plev_gepot_required
+
     if ( LPROF ) call start_timing( id, 'pres_lev_diags_alg' )
 
     nplev = get_axis_dimension('pressure_levels')
@@ -257,9 +263,16 @@ contains
     end if
 
     ! Geopotential height on pressure levels
+    aviation_thickness_flag = diag_samp('aviation__geopot_thickness_850') .or. &
+                              diag_samp('aviation__geopot_thickness_500')
     plev_geopot_clim_flag = diag_samp('plev__geopot_clim')
-    plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', activate=plev_geopot_clim_flag)
-    if ((plev_geopot_flag .or. plev_geopot_clim_flag) .and. use_xios_io) then
+    plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', &
+            activate=plev_geopot_clim_flag .or. aviation_thickness_flag)
+
+    plev_gepot_required = plev_geopot_flag .or. &
+                          plev_geopot_clim_flag .or. &
+                          aviation_thickness_flag
+    if (plev_gepot_required .and. use_xios_io) then
 
       height_w3 => get_height_fv(W3, exner_w3%get_mesh_id())
       call invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, &

interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90

@@ -0,0 +1,133 @@
+!-------------------------------------------------------------------------------
+! (c) Crown copyright 2021 Met Office. All rights reserved.
+! The file LICENCE, distributed with this code, contains details of the terms
+! under which the code may be used.
+!-------------------------------------------------------------------------------
+!
+! Section 20 aviation diagnostics kernel.
+!
+! 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_kernel_mod
+
+  USE argument_mod,         ONLY: arg_type,            &
+                                  GH_FIELD, GH_SCALAR, GH_LOGICAL, &
+                                  GH_READ, GH_WRITE, GH_INTEGER, &
+                                  GH_REAL, CELL_COLUMN, &
+                                  ANY_DISCONTINUOUS_SPACE_1, &
+                                  ANY_DISCONTINUOUS_SPACE_2
+  USE kernel_mod,           ONLY: kernel_type
+  USE constants_mod,        ONLY: r_def, i_def, l_def
+
+  IMPLICIT NONE
+
+  ! The aviation diagnostics kernel type.
+  TYPE, EXTENDS(kernel_type) :: aviation_diags_kernel_type
+    TYPE(arg_type), DIMENSION(8) :: meta_args = (/ &
+
+      ! Output fields.
+      arg_type(gh_field, gh_real, gh_write, ANY_DISCONTINUOUS_SPACE_1), &
+      arg_type(gh_field, gh_real, gh_write, ANY_DISCONTINUOUS_SPACE_1), &
+
+      ! Source field.
+      arg_type(gh_field, gh_real, gh_read, ANY_DISCONTINUOUS_SPACE_2), &
+
+      ! Request flags.
+      arg_type(gh_scalar, gh_logical, gh_read), &
+      arg_type(gh_scalar, gh_logical, gh_read), &
+
+      ! Level indices.
+      arg_type(gh_scalar, gh_integer, gh_read), &
+      arg_type(gh_scalar, gh_integer, gh_read), &
+      arg_type(gh_scalar, gh_integer, gh_read) &
+      /)
+
+    INTEGER :: operates_on = cell_column
+
+    CONTAINS
+      PROCEDURE, NOPASS :: code => aviation_diags_kernel_code
+    END TYPE aviation_diags_kernel_type
+
+CONTAINS
+
+  SUBROUTINE aviation_diags_kernel_code(nlayers, &
+             ! Output fields.
+             thickness_850, thickness_500, &
+
+             ! Source field.
+             plev_geopot, &
+
+             ! Request flags.
+             thickness_850_flag, thickness_500_flag, &
+
+             ! Level incides.
+             i1000, i850, i500, &
+
+             ! Kernel stuff.
+             result_ndf, result_undf, result_map, &
+             source_ndf, source_undf, source_map)
+
+    ! Subtract geopotential heights at 850 and 500 hPa from that at 1000 hPa
+    ! to calculate thickness fields.
+
+    IMPLICIT NONE
+
+    ! Arguments (kernel)
+
+    ! The number of layers in a column.
+    INTEGER(KIND=i_def), INTENT(IN) :: nlayers
+
+    ! Number of degrees of freedom (columns) in the cell we're processing.
+    INTEGER(KIND=i_def), INTENT(IN) :: result_ndf, source_ndf
+
+    ! Number of unique degrees of freedom in the fields.
+    INTEGER(KIND=i_def), INTENT(IN) :: result_undf, source_undf
+
+    ! Degrees of freedom maps. Offsets to the bottom of each column.
+    INTEGER(KIND=i_def), INTENT(IN), DIMENSION(result_ndf) :: result_map
+    INTEGER(KIND=i_def), INTENT(IN), DIMENSION(source_ndf) :: source_map
+
+
+    ! Arguments (algorithm)
+
+    ! Output thickness fields.
+    REAL(KIND=r_def), INTENT(OUT), DIMENSION(result_undf) :: thickness_850
+    REAL(KIND=r_def), INTENT(OUT), DIMENSION(result_undf) :: thickness_500
+
+    ! Geopotential height at pressure levels.
+    REAL(KIND=r_def), INTENT(IN), DIMENSION(source_undf) :: plev_geopot
+
+    ! Request flags.
+    LOGICAL(KIND=l_def), INTENT(IN) :: thickness_850_flag, thickness_500_flag
+
+    ! Level indices. For i850 and i500, -1 means "not requested".
+    INTEGER(KIND=i_def), INTENT(IN) :: i1000, i850, i500
+
+
+    ! Local variables
+    INTEGER(KIND=i_def) :: df
+    REAL(KIND=r_def) :: gph_1000
+
+
+    ! Process every DOF in this cell.
+    DO df = 1, result_ndf
+
+      gph_1000 = plev_geopot(source_map(df) + i1000-1)
+
+      IF (i850 /= -1) THEN
+        thickness_850(result_map(df)) = &
+          plev_geopot(source_map(df)+i850-1) - gph_1000
+      END IF
+
+      IF (i500 /= -1) THEN
+        thickness_500(result_map(df)) = &
+          plev_geopot(source_map(df)+i500-1) - gph_1000
+      END IF
+
+    END DO
+
+  END SUBROUTINE aviation_diags_kernel_code
+
+END MODULE aviation_diags_kernel_mod