From f0a35d64be73f4d707e16d510e42a196dcca6e3b Mon Sep 17 00:00:00 2001 From: bblay Date: Fri, 12 Dec 2025 16:33:20 +0000 Subject: [PATCH 01/26] s20 diag thickness --- CONTRIBUTORS.md | 3 +- applications/lfric_atm/example/iodef.xml | 5 + .../lfric_atm/metadata/field_def_diags.xml | 5 + .../algorithm/aviation_diags_alg_mod.x90 | 176 ++++++++++++++++++ .../algorithm/pres_lev_diags_alg_mod.x90 | 23 ++- .../kernel/aviation_diags_kernel_mod.F90 | 134 +++++++++++++ .../file/file_def_diags_oper_nwp_gl.xml | 9 +- .../driver/gungho_diagnostics_driver_mod.F90 | 10 +- 8 files changed, 357 insertions(+), 8 deletions(-) create mode 100644 interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 create mode 100644 interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 00c6e6edb..7f3d1c142 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -9,6 +9,7 @@ | mo-rickywong | Ricky Wong | Met Office | 2026-24-02 | | 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 | | ricky-lv426 | Ricky Olivier | University of Exeter | 2026-01-12 | @@ -34,4 +35,4 @@ | thomasmelvin | Thomas Melvin | Met Office | 2026-01-15 | | tinyendian | Wolfgang Hayek | Earth Sciences New Zealand | 2026-02-02 | | DanStoneMO | Daniel Stone | Met Office | 2026-02-26 | -| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | \ No newline at end of file +| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | diff --git a/applications/lfric_atm/example/iodef.xml b/applications/lfric_atm/example/iodef.xml index 96653326d..adb66d722 100644 --- a/applications/lfric_atm/example/iodef.xml +++ b/applications/lfric_atm/example/iodef.xml @@ -94,6 +94,11 @@ + + + + + diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index a17d884f9..9f941bcb6 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -991,6 +991,11 @@ + + + + + diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 new file mode 100644 index 000000000..a9b24ce4a --- /dev/null +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -0,0 +1,176 @@ +!----------------------------------------------------------------------------- +! (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, str_long + 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 + + character(len=str_long) :: message + + + ! 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') + + ! Anything to do? + IF ( .NOT. (aviation_thick_850_flag .OR. aviation_thick_500_flag) ) THEN + ! Nothing requested at this time step. + RETURN + END IF + + IF ( .NOT. use_xios_io ) THEN + RETURN + ENDIF + + ! Initialise required fields to 1. Todo: Why 1 and not 0? + IF ( aviation_thick_850_flag ) THEN + message = 'Section 20: geopot thick 850 requested' + CALL log_event( message, LOG_LEVEL_DEBUG ) + CALL invoke( setval_c(aviation_thick_850, 1.0_r_def)) + END IF + + IF ( aviation_thick_500_flag ) THEN + message = 'Section 20: geopot thick 500 requested' + CALL log_event( message, LOG_LEVEL_DEBUG ) + CALL invoke( setval_c(aviation_thick_500, 1.0_r_def)) + END IF + + + ! Get the array of pressure levels. + nplev = get_axis_dimension('pressure_levels') + IF (nplev <= 0) THEN + message = 'Section 20: No pressure levels' + CALL log_event( message, LOG_LEVEL_DEBUG ) + RETURN + END IF + + ALLOCATE(plevs(nplev), stat=plevs_alloc_stat) + IF (plevs_alloc_stat /= 0) THEN + message = 'Section 20: allocate(plevs) failed' + CALL log_event( message, LOG_LEVEL_DEBUG ) + + RETURN + END IF + plevs = get_axis_values('pressure_levels',nplev) + + ! Find the level indices for our three pressures of interest. + ! Assumes hPa. + 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 + message = 'Section 20: could not find 1000hPa' + CALL log_event( message, LOG_LEVEL_INFO ) + RETURN + END IF + + IF (aviation_thick_850_flag .AND. i850 == -1) THEN + message = 'Section 20: could not find 850hPa' + CALL log_event( message, LOG_LEVEL_INFO ) + RETURN + END IF + + IF (aviation_thick_500_flag .AND. i500 == -1) THEN + message = 'Section 20: could not find 500hPa' + CALL log_event( message, LOG_LEVEL_INFO ) + 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( 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 diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index 0907509b0..62bbd9038 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/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, 'diags.pressure_lev' ) 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, & diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 new file mode 100644 index 000000000..1cde71126 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 @@ -0,0 +1,134 @@ +!------------------------------------------------------------------------------- +! (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. +!------------------------------------------------------------------------------- +! +! 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 (thickness_850_flag .and. i850 /= -1) THEN + thickness_850(result_map(df)) = & + plev_geopot(source_map(df)+i850-1) - gph_1000 + END IF + + IF (thickness_500_flag .and. 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 diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index fb5ee2571..d9b7c7d7c 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -190,7 +190,14 @@ - + + + + + + + + diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index 1bfca57a2..6fb2ff97e 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -59,6 +59,7 @@ module gungho_diagnostics_driver_mod use pmsl_alg_mod, only : pmsl_alg use rh_diag_alg_mod, only : rh_diag_alg use freeze_lev_alg_mod, only : freeze_lev_alg + use aviation_diags_alg_mod, only : aviation_diags_alg #endif implicit none @@ -95,6 +96,11 @@ subroutine gungho_diagnostics_driver( modeldb, & type(field_type), pointer :: moist_dyn(:) type(field_collection_type), pointer :: derived_fields + ! For S20 Aviation diagnostics +#ifdef UM_PHYSICS + type( field_type ) :: plev_geopot ! Set by pres_lev_diags_alg(). +#endif + type(field_type), pointer :: theta type(field_type), pointer :: u type(field_type), pointer :: h_u @@ -332,9 +338,11 @@ subroutine gungho_diagnostics_driver( modeldb, & ! Call PMSL algorithm call pmsl_alg(exner, derived_fields, theta, twod_mesh) ! Pressure level diagnostics - call pres_lev_diags_alg(derived_fields, theta, exner, mr, moist_dyn) + call pres_lev_diags_alg(derived_fields, theta, exner, mr, moist_dyn, plev_geopot) ! Wet bulb freezing level call freeze_lev_alg(theta, mr, moist_dyn, exner_in_wth) + ! Aviation diagnostics + call aviation_diags_alg(plev_geopot) #endif temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') From 6c2c3b7c8f297a877f442f83182fefee1d4ec523 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 30 Apr 2026 12:54:01 +0100 Subject: [PATCH 02/26] tidy --- .../source/algorithm/aviation_diags_alg_mod.x90 | 15 +++++++-------- .../source/algorithm/pres_lev_diags_alg_mod.x90 | 6 +++--- .../source/kernel/aviation_diags_kernel_mod.F90 | 2 +- 3 files changed, 11 insertions(+), 12 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index a9b24ce4a..300efba9c 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -14,11 +14,11 @@ MODULE aviation_diags_alg_mod USE constants_mod, ONLY: r_def, i_def, l_def, str_long USE field_mod, ONLY: field_type - USE io_config_mod, ONLY: USE_xios_io + 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, & + USE log_mod, ONLY: log_event, & LOG_LEVEL_ALWAYS, LOG_LEVEL_ERROR, & LOG_LEVEL_INFO, LOG_LEVEL_DEBUG @@ -100,15 +100,14 @@ CONTAINS nplev = get_axis_dimension('pressure_levels') IF (nplev <= 0) THEN message = 'Section 20: No pressure levels' - CALL log_event( message, LOG_LEVEL_DEBUG ) + CALL log_event( message, LOG_LEVEL_ERROR ) RETURN END IF ALLOCATE(plevs(nplev), stat=plevs_alloc_stat) IF (plevs_alloc_stat /= 0) THEN message = 'Section 20: allocate(plevs) failed' - CALL log_event( message, LOG_LEVEL_DEBUG ) - + CALL log_event( message, LOG_LEVEL_ERROR ) RETURN END IF plevs = get_axis_values('pressure_levels',nplev) @@ -138,19 +137,19 @@ CONTAINS ! Check we found the required levels. IF (i1000 == -1) THEN message = 'Section 20: could not find 1000hPa' - CALL log_event( message, LOG_LEVEL_INFO ) + CALL log_event( message, LOG_LEVEL_ERROR ) RETURN END IF IF (aviation_thick_850_flag .AND. i850 == -1) THEN message = 'Section 20: could not find 850hPa' - CALL log_event( message, LOG_LEVEL_INFO ) + CALL log_event( message, LOG_LEVEL_ERROR ) RETURN END IF IF (aviation_thick_500_flag .AND. i500 == -1) THEN message = 'Section 20: could not find 500hPa' - CALL log_event( message, LOG_LEVEL_INFO ) + CALL log_event( message, LOG_LEVEL_ERROR ) RETURN END IF diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index 62bbd9038..a2fe4353a 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -57,9 +57,9 @@ contains ! Arguments type( field_collection_type ), intent(in) :: derived_fields - 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(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 diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 index 1cde71126..ca19c7bf3 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 @@ -103,7 +103,7 @@ SUBROUTINE aviation_diags_kernel_code(nlayers, & ! Request flags. LOGICAL(KIND=l_def), INTENT(IN) :: thickness_850_flag, thickness_500_flag - ! Level indices. For i850 and i500, -1 means "not requested". + ! Level indices. INTEGER(KIND=i_def), INTENT(IN) :: i1000, i850, i500 From e3cdd7d8a7be404433a53e2430013caac0ff9a0e Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 30 Apr 2026 15:56:51 +0100 Subject: [PATCH 03/26] lower case --- .../algorithm/aviation_diags_alg_mod.x90 | 148 +++++++++--------- .../kernel/aviation_diags_kernel_mod.F90 | 94 +++++------ 2 files changed, 121 insertions(+), 121 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index 300efba9c..d6dce9ed4 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -10,58 +10,58 @@ ! This file currently belongs in section: physics_schemes_interface ! whilst discussions are ongoing about its final location. ! -MODULE aviation_diags_alg_mod +module aviation_diags_alg_mod - USE constants_mod, ONLY: r_def, i_def, l_def, str_long - 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 constants_mod, only: r_def, i_def, l_def, str_long + 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_LEVEL_ALWAYS, LOG_LEVEL_ERROR, & - LOG_LEVEL_INFO, LOG_LEVEL_DEBUG + use lfric_xios_diag_mod, only: get_axis_dimension, get_axis_values + use log_mod, only: log_event, & + log_level_always, log_level_error, & + log_level_info, log_level_debug - USE aviation_diags_kernel_mod, ONLY: aviation_diags_kernel_type + use aviation_diags_kernel_mod, only: aviation_diags_kernel_type - IMPLICIT NONE + implicit none -PRIVATE -PUBLIC :: aviation_diags_alg +private +public :: aviation_diags_alg -CONTAINS +contains ! Algorithm to calculate the section 20 aviation diagnostics. - SUBROUTINE aviation_diags_alg(plev_geopot) + subroutine aviation_diags_alg(plev_geopot) - IMPLICIT NONE + implicit none ! Arguments ! We calculate the result from this input field. - TYPE(field_type), INTENT(IN) :: plev_geopot + 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 + 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 + 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 + 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 + 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 + type( field_type ) :: aviation_thick_850, aviation_thick_500 - INTEGER(I_DEF) :: i + integer(i_def) :: i character(len=str_long) :: message @@ -73,43 +73,43 @@ CONTAINS init_diag(aviation_thick_500, 'aviation__geopot_thickness_500') ! Anything to do? - IF ( .NOT. (aviation_thick_850_flag .OR. aviation_thick_500_flag) ) THEN + if ( .not. (aviation_thick_850_flag .or. aviation_thick_500_flag) ) then ! Nothing requested at this time step. - RETURN - END IF + return + end if - IF ( .NOT. use_xios_io ) THEN - RETURN - ENDIF + if ( .not. use_xios_io ) then + return + endif ! Initialise required fields to 1. Todo: Why 1 and not 0? - IF ( aviation_thick_850_flag ) THEN + if ( aviation_thick_850_flag ) then message = 'Section 20: geopot thick 850 requested' - CALL log_event( message, LOG_LEVEL_DEBUG ) - CALL invoke( setval_c(aviation_thick_850, 1.0_r_def)) - END IF + call log_event( message, log_level_debug ) + call invoke( setval_c(aviation_thick_850, 1.0_r_def)) + end if - IF ( aviation_thick_500_flag ) THEN + if ( aviation_thick_500_flag ) then message = 'Section 20: geopot thick 500 requested' - CALL log_event( message, LOG_LEVEL_DEBUG ) - CALL invoke( setval_c(aviation_thick_500, 1.0_r_def)) - END IF + call log_event( message, log_level_debug ) + call invoke( setval_c(aviation_thick_500, 1.0_r_def)) + end if ! Get the array of pressure levels. nplev = get_axis_dimension('pressure_levels') - IF (nplev <= 0) THEN + if (nplev <= 0) then message = 'Section 20: No pressure levels' - CALL log_event( message, LOG_LEVEL_ERROR ) - RETURN - END IF + call log_event( message, log_level_error ) + return + end if - ALLOCATE(plevs(nplev), stat=plevs_alloc_stat) - IF (plevs_alloc_stat /= 0) THEN + allocate(plevs(nplev), stat=plevs_alloc_stat) + if (plevs_alloc_stat /= 0) then message = 'Section 20: allocate(plevs) failed' - CALL log_event( message, LOG_LEVEL_ERROR ) - RETURN - END IF + call log_event( message, log_level_error ) + return + end if plevs = get_axis_values('pressure_levels',nplev) ! Find the level indices for our three pressures of interest. @@ -117,59 +117,59 @@ CONTAINS i1000 = -1 i850 = -1 i500 = -1 - DO i = 1, nplev + do i = 1, nplev ! 1000? - IF ( abs(plevs(i) - 100000.0_r_def) < plev_tol ) THEN + 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 + 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 + else if ( abs(plevs(i) - 50000.0_r_def) < plev_tol ) then i500 = i - END IF - END DO + end if + end do ! Check we found the required levels. - IF (i1000 == -1) THEN + if (i1000 == -1) then message = 'Section 20: could not find 1000hPa' - CALL log_event( message, LOG_LEVEL_ERROR ) - RETURN - END IF + call log_event( message, log_level_error ) + return + end if - IF (aviation_thick_850_flag .AND. i850 == -1) THEN + if (aviation_thick_850_flag .and. i850 == -1) then message = 'Section 20: could not find 850hPa' - CALL log_event( message, LOG_LEVEL_ERROR ) - RETURN - END IF + call log_event( message, log_level_error ) + return + end if - IF (aviation_thick_500_flag .AND. i500 == -1) THEN + if (aviation_thick_500_flag .and. i500 == -1) then message = 'Section 20: could not find 500hPa' - CALL log_event( message, LOG_LEVEL_ERROR ) - RETURN - END IF + call log_event( message, 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( aviation_diags_kernel_type( & + call invoke( 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() + 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 + if (allocated(plevs)) then + deallocate(plevs) + end if - END SUBROUTINE aviation_diags_alg + end subroutine aviation_diags_alg -END MODULE aviation_diags_alg_mod +end module aviation_diags_alg_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 index ca19c7bf3..8cb79805e 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 @@ -10,30 +10,30 @@ ! This file currently belongs in section: physics_schemes_interface ! whilst discussions are ongoing about its final location. ! -MODULE aviation_diags_kernel_mod +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 - + 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 + + implicit none ! The aviation diagnostics kernel type. - TYPE, EXTENDS(kernel_type) :: aviation_diags_kernel_type - TYPE(arg_type), DIMENSION(8) :: meta_args = (/ & + 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), & + 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), & + arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_2), & ! Request flags. arg_type(gh_scalar, gh_logical, gh_read), & @@ -45,90 +45,90 @@ MODULE aviation_diags_kernel_mod arg_type(gh_scalar, gh_integer, gh_read) & /) - INTEGER :: operates_on = cell_column + integer :: operates_on = cell_column - CONTAINS - PROCEDURE, NOPASS :: code => aviation_diags_kernel_code - END TYPE aviation_diags_kernel_type + contains + procedure, nopass :: code => aviation_diags_kernel_code + end type aviation_diags_kernel_type -CONTAINS +contains - SUBROUTINE aviation_diags_kernel_code(nlayers, & + subroutine aviation_diags_kernel_code(nlayers, & ! Output fields. - thickness_850, thickness_500, & + thickness_850, thickness_500, & ! Source field. - plev_geopot, & + plev_geopot, & ! Request flags. thickness_850_flag, thickness_500_flag, & ! Level incides. - i1000, i850, i500, & + i1000, i850, i500, & ! Kernel stuff. - result_ndf, result_undf, result_map, & + 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 + implicit none ! Arguments (kernel) ! The number of layers in a column. - INTEGER(KIND=i_def), INTENT(IN) :: nlayers + 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 + 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 + 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 + 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 + 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 + 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 + logical(kind=l_def), intent(in) :: thickness_850_flag, thickness_500_flag ! Level indices. - INTEGER(KIND=i_def), INTENT(IN) :: i1000, i850, i500 + integer(kind=i_def), intent(in) :: i1000, i850, i500 ! Local variables - INTEGER(KIND=i_def) :: df - REAL(KIND=r_def) :: gph_1000 + integer(kind=i_def) :: df + real(kind=r_def) :: gph_1000 ! Process every DOF in this cell. - DO df = 1, result_ndf + do df = 1, result_ndf gph_1000 = plev_geopot(source_map(df) + i1000-1) - IF (thickness_850_flag .and. i850 /= -1) THEN + if (thickness_850_flag .and. i850 /= -1) then thickness_850(result_map(df)) = & plev_geopot(source_map(df)+i850-1) - gph_1000 - END IF + end if - IF (thickness_500_flag .and. i500 /= -1) THEN + if (thickness_500_flag .and. i500 /= -1) then thickness_500(result_map(df)) = & plev_geopot(source_map(df)+i500-1) - gph_1000 - END IF + end if - END DO + end do - END SUBROUTINE aviation_diags_kernel_code + end subroutine aviation_diags_kernel_code -END MODULE aviation_diags_kernel_mod +end module aviation_diags_kernel_mod From 474430725d811f1855157473f3bbf384928d50fb Mon Sep 17 00:00:00 2001 From: bblay Date: Fri, 12 Dec 2025 16:33:20 +0000 Subject: [PATCH 04/26] S20 icao heights --- applications/lfric_atm/example/iodef.xml | 2 + .../lfric_atm/metadata/field_def_diags.xml | 3 +- .../algorithm/aviation_diags_alg_mod.x90 | 105 +++++++++++-- .../source/algorithm/conv_gr_alg_mod.x90 | 32 +++- .../source/diagnostics/conv_diags_mod.x90 | 11 +- .../source/kernel/diags_geopot_kernel_mod.F90 | 133 +++++++++++++++++ .../kernel/diags_icao_heights_kernel_mod.F90 | 139 ++++++++++++++++++ .../file/file_def_diags_oper_nwp_gl.xml | 2 + .../driver/gungho_diagnostics_driver_mod.F90 | 2 +- .../science_aviation_constants_mod.F90 | 55 +++++++ 10 files changed, 461 insertions(+), 23 deletions(-) create mode 100644 interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 create mode 100644 interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 create mode 100644 science/shared/source/constants/science_aviation_constants_mod.F90 diff --git a/applications/lfric_atm/example/iodef.xml b/applications/lfric_atm/example/iodef.xml index adb66d722..e472a05d4 100644 --- a/applications/lfric_atm/example/iodef.xml +++ b/applications/lfric_atm/example/iodef.xml @@ -97,6 +97,8 @@ + + diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 9f941bcb6..e07083425 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -991,10 +991,11 @@ - + + diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index d6dce9ed4..ebb8f622a 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -12,8 +12,9 @@ ! module aviation_diags_alg_mod - use constants_mod, only: r_def, i_def, l_def, str_long + use constants_mod, only: r_def, i_def, l_def, str_long, 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 @@ -22,7 +23,8 @@ module aviation_diags_alg_mod log_level_always, log_level_error, & log_level_info, log_level_debug - use aviation_diags_kernel_mod, only: aviation_diags_kernel_type + use diags_geopot_kernel_mod, only: diags_geopot_kernel_type + use diags_icao_heights_kernel_mod, only: diags_icao_heights_kernel_type implicit none @@ -32,15 +34,27 @@ public :: aviation_diags_alg contains ! Algorithm to calculate the section 20 aviation diagnostics. - subroutine aviation_diags_alg(plev_geopot) - + subroutine aviation_diags_alg(plev_geopot, convection_fields) implicit none ! Arguments - - ! We calculate the result from this input field. 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 @@ -140,7 +154,7 @@ contains call log_event( message, log_level_error ) return end if - + if (aviation_thick_850_flag .and. i850 == -1) then message = 'Section 20: could not find 850hPa' call log_event( message, log_level_error ) @@ -155,7 +169,7 @@ contains ! 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( aviation_diags_kernel_type( & + call invoke( diags_geopot_kernel_type( & aviation_thick_850, aviation_thick_500, & plev_geopot, & aviation_thick_850_flag, aviation_thick_500_flag, & @@ -170,6 +184,79 @@ contains deallocate(plevs) end if - end subroutine aviation_diags_alg + 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_base_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 diff --git a/interfaces/physics_schemes_interface/source/algorithm/conv_gr_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/conv_gr_alg_mod.x90 index 68a06a6ab..6fd735a27 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/conv_gr_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/conv_gr_alg_mod.x90 @@ -28,6 +28,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 @@ -63,14 +65,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 @@ -776,6 +778,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 ) diff --git a/interfaces/physics_schemes_interface/source/diagnostics/conv_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/conv_diags_mod.x90 index 705a22118..24dd6c318 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/conv_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/conv_diags_mod.x90 @@ -10,8 +10,8 @@ module conv_diags_mod use constants_mod, only: l_def use field_mod, only: field_type use timing_mod, only: start_timing, stop_timing, tik, LPROF - use initialise_diagnostics_mod, only : init_diag & - => init_diagnostic_field + use initialise_diagnostics_mod, only : init_diag => init_diagnostic_field, & + diag_samp => diagnostic_to_be_sampled implicit none @@ -197,13 +197,16 @@ contains ! not zeroed because overwritten ! cv_base and _top are dependencies for chemistry, hence need to be 'active' + ! They are also used for aviation diagnostics. ! cv_base is zeroed because value tested in kernel call invoke( setval_c(cv_base, 0.0_r_def) ) - pres_cv_base_flag = init_diag(pres_cv_base, 'convection__pres_cv_base') + pres_cv_base_flag = init_diag(pres_cv_base, 'convection__pres_cv_base', & + activate=diag_samp('aviation__conv_cloud_base_icao_height')) ! not zeroed because overwritten - pres_cv_top_flag = init_diag(pres_cv_top, 'convection__pres_cv_top') + pres_cv_top_flag = init_diag(pres_cv_top, 'convection__pres_cv_top', & + activate=diag_samp('aviation__conv_cloud_top_icao_height')) ! not zeroed because overwritten pres_lowest_cv_base_flag = init_diag(pres_lowest_cv_base, & diff --git a/interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 new file mode 100644 index 000000000..4489b53f8 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/diags_geopot_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 diags_geopot_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) :: diags_geopot_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 => diags_geopot_kernel_code + END TYPE diags_geopot_kernel_type + +CONTAINS + + SUBROUTINE diags_geopot_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 diags_geopot_kernel_code + +END MODULE diags_geopot_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 new file mode 100644 index 000000000..ebcbc5ca1 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 @@ -0,0 +1,139 @@ +!------------------------------------------------------------------------------- +! (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 diags_icao_heights_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 constants_mod, ONLY: r_def, i_def, l_def + USE kernel_mod, ONLY: kernel_type + + IMPLICIT NONE + + ! The aviation diagnostics kernel type. + TYPE, EXTENDS(kernel_type) :: diags_icao_heights_kernel_type + TYPE(arg_type), DIMENSION(2) :: meta_args = (/ & + + ! Output icao height. + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & + + ! Input pressure + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1) & + + /) + + INTEGER :: operates_on = cell_column + + CONTAINS + PROCEDURE, NOPASS :: code => diags_icao_heights_kernel_code + END TYPE diags_icao_heights_kernel_type + +CONTAINS + + SUBROUTINE diags_icao_heights_kernel_code(nlayers, & + icao_height, & ! Output icao height. + pressure_field, & ! Pressure. + ndf, undf, map) + + ! Calculate iaco height from the pressure field. + ! Assumes lowest order W3 data, where ndf is always 1. + + USE constants_mod, ONLY: RMDI + USE planet_config_mod, ONLY: gravity, rd + USE science_conversions_mod, ONLY: feet_to_metres + USE science_aviation_constants_mod, ONLY: mtokft, & + isa_Lapse_RateL, isa_Lapse_RateU, isa_Press_Bot, & + isa_Press_Mid, isa_Press_Top, isa_Temp_Bot, isa_Temp_Top, & + Gpm1, Gpm2 + + 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) :: ndf + + ! Number of unique degrees of freedom in the fields. + INTEGER(KIND=i_def), INTENT(IN) :: undf + + ! Degrees of freedom maps. Offsets to the bottom of each column. + INTEGER(KIND=i_def), INTENT(IN), DIMENSION(ndf) :: map + + + ! Arguments (algorithm) + + ! Output icao height. + REAL(KIND=r_def), INTENT(INOUT), DIMENSION(undf) :: icao_height + + ! Pressure in Pa. + REAL(KIND=r_def), INTENT(IN), DIMENSION(undf) :: pressure_field + + + ! Local variables + INTEGER(KIND=i_def) :: df + + REAL(KIND=r_def) :: g_over_rd + REAL(KIND=r_def) :: zp1 + REAL(KIND=r_def) :: zp2 + + REAL(KIND=r_def) :: pressure + + + g_over_rd = gravity / rd + zp1 = isa_Lapse_RateL / g_over_rd + zp2 = isa_Lapse_RateU / g_over_rd + + pressure = pressure_field(map(1)) + + ! Setting a safeguard limit to the lowest pressure to prevent + ! extremely large icao height values near the top at the + ! atmosphere. + if ( (pressure >= 0.0_r_def) .and. (pressure <= 1000.0_r_def) ) then + pressure = 1000.0_r_def + end if + + ! Pressure must not be greater than surface pressure. + pressure = MIN(isa_Press_bot, pressure) + + ! Missing or invalid data? + if (pressure == RMDI .or. pressure <= 0.0_r_def) then + icao_height(map(1)) = RMDI + + ! Heights up to 11,000 GPM + else if ( pressure > isa_Press_Mid ) then + pressure = pressure / isa_Press_Bot + pressure = 1.0_r_def - pressure**zp1 + icao_height(map(1)) = pressure * isa_Temp_Bot / isa_Lapse_RateL + + ! Heights between 11,000 GPM and 20,000 GPM + else if ( pressure > isa_Press_Top ) then + pressure = pressure / isa_Press_Mid + pressure = -LOG(pressure) + icao_height(map(1)) = Gpm1 + pressure * isa_Temp_Top / g_over_rd + + ! Heights above 20,000 GPM + else + pressure = pressure / isa_Press_Top + pressure = 1.0_r_def - pressure**zp2 + icao_height(map(1)) = Gpm2 + pressure * isa_Temp_Top / isa_Lapse_RateU + end if + + END SUBROUTINE diags_icao_heights_kernel_code + +END MODULE diags_icao_heights_kernel_mod diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index d9b7c7d7c..6a4d5fa0f 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -196,6 +196,8 @@ + + diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index 6fb2ff97e..2f710b093 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -342,7 +342,7 @@ subroutine gungho_diagnostics_driver( modeldb, & ! Wet bulb freezing level call freeze_lev_alg(theta, mr, moist_dyn, exner_in_wth) ! Aviation diagnostics - call aviation_diags_alg(plev_geopot) + call aviation_diags_alg(plev_geopot, modeldb%fields%get_field_collection("convection_fields")) #endif temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') diff --git a/science/shared/source/constants/science_aviation_constants_mod.F90 b/science/shared/source/constants/science_aviation_constants_mod.F90 new file mode 100644 index 000000000..0176d9878 --- /dev/null +++ b/science/shared/source/constants/science_aviation_constants_mod.F90 @@ -0,0 +1,55 @@ +!---------------------------------------------------------------------------- +! (c) Crown copyright 2018 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 LFRic aviation constants module +!---------------------------------------------------------------------------- +! +MODULE science_aviation_constants_mod + + USE constants_mod, ONLY: r_def + USE science_conversions_mod, ONLY: feet_to_metres + + IMPLICIT NONE + + PRIVATE + PUBLIC :: mtokft, isa_Lapse_RateL, isa_Lapse_RateU, & + isa_Press_Bot, isa_Press_Mid, isa_Press_Top, & + isa_Temp_Bot, isa_Temp_Top, Gpm1, Gpm2 + + ! meters to thousands of feet + REAL(r_def), PARAMETER :: mtokft = 1.0_r_def / (feet_to_metres * 1000.0_r_def) + + ! International Standard Atmosphere (isa) within troposhpere for: + ! Values taken from: + ! https://en.wikipedia.org/wiki/International_Standard_Atmosphere + ! Lapse rate (degreeC/km) for levels below 11,000 gpm + REAL(r_def), PARAMETER :: isa_Lapse_RateL = 6.5e-03_r_def + + ! Lapse rate (degreeC/km) for levels above 11,000 gpm + REAL(r_def), PARAMETER :: isa_Lapse_RateU = -1.0e-03_r_def + + ! surface pressure (Pa) + REAL(r_def), PARAMETER :: isa_Press_Bot = 101325.0_r_def + + ! pressure (Pa) at 11,000 gpm + REAL(r_def), PARAMETER :: isa_Press_Mid = 22632.0_r_def + + ! pressure (Pa) at 20,000 gpm + REAL(r_def), PARAMETER :: isa_Press_Top = 5474.87_r_def + + ! Surface temperature (K) + REAL(r_def), PARAMETER :: isa_Temp_Bot = 288.15_r_def + + ! Temperature (K) of isothermal layer - at tropopause + REAL(r_def), PARAMETER :: isa_Temp_Top = 216.65_r_def + + ! Height limit (gpm) for std lower lapse rate + REAL(r_def), PARAMETER :: Gpm1 = 11000.0_r_def + + ! Height (gpm) of top of isothermal layer + REAL(r_def), PARAMETER :: Gpm2 = 20000.0_r_def + +END MODULE science_aviation_constants_mod + From dbb50d118ad00cf98fbb455ec1c59822ec0f58a2 Mon Sep 17 00:00:00 2001 From: bblay Date: Fri, 23 Jan 2026 11:55:01 +0000 Subject: [PATCH 05/26] running to completion locally --- .../source/algorithm/aviation_diags_alg_mod.x90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index ebb8f622a..3bac72ccc 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -250,7 +250,7 @@ contains diags_icao_heights_kernel_type(conv_cloud_top_icao_height, & conv_cloud_top_pressure)) - call conv_cloud_base_icao_height%write_field() + 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') From f753421be4c6a3f867d0c392b6ef67c326370624 Mon Sep 17 00:00:00 2001 From: bblay Date: Fri, 12 Dec 2025 16:33:20 +0000 Subject: [PATCH 06/26] S20 snow prob --- applications/lfric_atm/example/iodef.xml | 1 + .../lfric_atm/metadata/field_def_diags.xml | 1 + .../algorithm/aviation_diags_alg_mod.x90 | 37 +++++++++++++------ .../kernel/aviation_diags_kernel_mod.F90 | 35 +++++++++++++----- .../file/file_def_diags_oper_nwp_gl.xml | 1 + 5 files changed, 53 insertions(+), 22 deletions(-) diff --git a/applications/lfric_atm/example/iodef.xml b/applications/lfric_atm/example/iodef.xml index e472a05d4..919b71e25 100644 --- a/applications/lfric_atm/example/iodef.xml +++ b/applications/lfric_atm/example/iodef.xml @@ -97,6 +97,7 @@ + diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index e07083425..2e2f8a10f 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -994,6 +994,7 @@ + diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index 3bac72ccc..c244daca7 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -59,7 +59,8 @@ contains ! 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 + logical(l_def) :: aviation_thick_850_flag, aviation_thick_500_flag, & + aviation_snow_probability_flag ! The array of pressure levels. integer(i_def) :: nplev @@ -72,8 +73,9 @@ contains ! 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 + ! The output fields we produce. + type( field_type ) :: aviation_thick_850, aviation_thick_500, & + aviation_snow_probability integer(i_def) :: i @@ -81,14 +83,18 @@ contains ! Check the request flags. + aviation_snow_probability_flag = & + init_diag(aviation_snow_probability, 'aviation__snow_probability') aviation_thick_850_flag = & - init_diag(aviation_thick_850, 'aviation__geopot_thickness_850') + init_diag(aviation_thick_850, 'aviation__geopot_thickness_850', & + activate=aviation_snow_probability_flag) aviation_thick_500_flag = & init_diag(aviation_thick_500, 'aviation__geopot_thickness_500') ! Anything to do? - if ( .not. (aviation_thick_850_flag .or. aviation_thick_500_flag) ) then - ! Nothing requested at this time step. + if ( .not. (aviation_thick_850_flag .or. & + aviation_thick_500_flag .or. & + aviation_snow_probability_flag) ) then return end if @@ -96,10 +102,17 @@ contains return endif - ! Initialise required fields to 1. Todo: Why 1 and not 0? + ! Initialise required fields to 1. + if ( aviation_snow_probability_flag ) then + message = 'Section 20: snow probability requested' + call log_event( message, log_level_debug ) + call invoke( setval_c(aviation_snow_probability, 0.0_r_def)) + end if + if ( aviation_thick_850_flag ) then message = 'Section 20: geopot thick 850 requested' call log_event( message, log_level_debug ) + ! Todo: Why 1 and not 0? call invoke( setval_c(aviation_thick_850, 1.0_r_def)) end if @@ -109,7 +122,6 @@ contains call invoke( setval_c(aviation_thick_500, 1.0_r_def)) end if - ! Get the array of pressure levels. nplev = get_axis_dimension('pressure_levels') if (nplev <= 0) then @@ -134,15 +146,15 @@ contains do i = 1, nplev ! 1000? - if ( abs(plevs(i) - 100000.0_r_def) < plev_tol ) then + 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 + 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 + else if (abs(plevs(i) - 50000.0_r_def) < plev_tol) then i500 = i end if @@ -173,11 +185,12 @@ contains aviation_thick_850, aviation_thick_500, & plev_geopot, & aviation_thick_850_flag, aviation_thick_500_flag, & - i1000, i850, i500)) + aviation_snow_probability_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() + if (aviation_snow_probability_flag) call aviation_snow_probability%write_field() ! Clean up. if (allocated(plevs)) then diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 index 8cb79805e..9464a6ee0 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 @@ -26,11 +26,12 @@ module aviation_diags_kernel_mod ! The aviation diagnostics kernel type. type, extends(kernel_type) :: aviation_diags_kernel_type - type(arg_type), dimension(8) :: meta_args = (/ & + type(arg_type), dimension(10) :: 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), & + 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), & @@ -38,6 +39,7 @@ module aviation_diags_kernel_mod ! Request flags. arg_type(gh_scalar, gh_logical, gh_read), & 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), & @@ -55,23 +57,27 @@ module aviation_diags_kernel_mod subroutine aviation_diags_kernel_code(nlayers, & ! Output fields. - thickness_850, thickness_500, & + thickness_850, thickness_500, snow_probability, & ! Source field. plev_geopot, & ! Request flags. - thickness_850_flag, thickness_500_flag, & + thickness_850_flag, thickness_500_flag, snow_probability_flag, & - ! Level incides. + ! Level indices. 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. + ! Thickness: + ! Subtract geopotential heights at 850 and 500 hPa from that at 1000 hPa. + ! + ! Snow probability: + ! Implement Boyden (1964), using 850 and 1000 hPa. + ! See https://github.com/MetOffice/Section20/issues/21 implicit none @@ -93,15 +99,18 @@ subroutine aviation_diags_kernel_code(nlayers, & ! Arguments (algorithm) - ! Output thickness fields. + ! Output fields. real(kind=r_def), intent(out), dimension(result_undf) :: thickness_850 real(kind=r_def), intent(out), dimension(result_undf) :: thickness_500 + real(kind=r_def), intent(out), dimension(result_undf) :: snow_probability ! 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 + logical(kind=l_def), intent(in) :: thickness_850_flag + logical(kind=l_def), intent(in) :: thickness_500_flag + logical(kind=l_def), intent(in) :: snow_probability_flag ! Level indices. integer(kind=i_def), intent(in) :: i1000, i850, i500 @@ -109,13 +118,14 @@ subroutine aviation_diags_kernel_code(nlayers, & ! Local variables integer(kind=i_def) :: df - real(kind=r_def) :: gph_1000 - + real(kind=r_def) :: gph_1000, gph_850, gph_500 ! Process every DOF in this cell. do df = 1, result_ndf gph_1000 = plev_geopot(source_map(df) + i1000-1) + gph_850 = plev_geopot(source_map(df) + i850-1) + gph_500 = plev_geopot(source_map(df) + i500-1) if (thickness_850_flag .and. i850 /= -1) then thickness_850(result_map(df)) = & @@ -127,6 +137,11 @@ subroutine aviation_diags_kernel_code(nlayers, & plev_geopot(source_map(df)+i500-1) - gph_1000 end if + if(snow_probability_flag) then + snow_probability(result_map(df)) = & + 5220.0 + 3.86666*gph_1000 - 4.0*gph_850 + end if + end do end subroutine aviation_diags_kernel_code diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index 6a4d5fa0f..8d70bf348 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -196,6 +196,7 @@ + From 436d9876c7262b8171a48450462fbb6039e57e52 Mon Sep 17 00:00:00 2001 From: bblay Date: Mon, 11 May 2026 16:57:31 +0100 Subject: [PATCH 07/26] running locally --- applications/lfric_atm/example/iodef.xml | 2 +- .../algorithm/aviation_diags_alg_mod.x90 | 18 ++- .../algorithm/pres_lev_diags_alg_mod.x90 | 16 +- .../aviation_icao_heights_kernel_mod.F90 | 139 ++++++++++++++++++ ...=> aviation_thickness_snow_kernel_mod.F90} | 61 ++++---- .../source/kernel/diags_geopot_kernel_mod.F90 | 133 ----------------- .../kernel/diags_icao_heights_kernel_mod.F90 | 139 ------------------ 7 files changed, 192 insertions(+), 316 deletions(-) create mode 100644 interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 rename interfaces/physics_schemes_interface/source/kernel/{aviation_diags_kernel_mod.F90 => aviation_thickness_snow_kernel_mod.F90} (75%) delete mode 100644 interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 delete mode 100644 interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 diff --git a/applications/lfric_atm/example/iodef.xml b/applications/lfric_atm/example/iodef.xml index 919b71e25..2440fbbbf 100644 --- a/applications/lfric_atm/example/iodef.xml +++ b/applications/lfric_atm/example/iodef.xml @@ -94,7 +94,7 @@ - + diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index c244daca7..2ea23df5c 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -23,8 +23,10 @@ module aviation_diags_alg_mod 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 + use aviation_thickness_snow_kernel_mod, & + only: aviation_thickness_snow_kernel_type + use aviation_icao_heights_kernel_mod, & + only: aviation_icao_heights_kernel_type implicit none @@ -181,8 +183,9 @@ contains ! 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( diags_geopot_kernel_type( & + call invoke( aviation_thickness_snow_kernel_type( & aviation_thick_850, aviation_thick_500, & + aviation_snow_probability, & plev_geopot, & aviation_thick_850_flag, aviation_thick_500_flag, & aviation_snow_probability_flag, i1000, i850, i500)) @@ -207,6 +210,7 @@ contains TYPE(field_collection_type), INTENT(INOUT) :: convection_fields ! Locals + character(len=str_long) :: message ! Pressure fields to unpack. type(field_type), pointer :: conv_cloud_base_pressure @@ -228,8 +232,8 @@ contains 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) + message = 'Section 20: cloud ICAO heights on' + call log_event( message, log_level_debug ) else ! nothing requested at this time step return @@ -242,7 +246,7 @@ contains '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, & + aviation_icao_heights_kernel_type(conv_cloud_base_icao_height, & conv_cloud_base_pressure)) call conv_cloud_base_icao_height%write_field() @@ -260,7 +264,7 @@ contains '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, & + aviation_icao_heights_kernel_type(conv_cloud_top_icao_height, & conv_cloud_top_pressure)) call conv_cloud_top_icao_height%write_field() diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index a2fe4353a..d9d3c828a 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -83,8 +83,11 @@ 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, & - aviation_thickness_flag + plev_geopot_flag, plev_geopot_clim_flag, plev_thetaw_flag + + ! Aviation diagnostics require plev_geopot. + logical(l_def) :: aviation_geopot_flag + integer(tik) :: id ! For code tidiness, combines all flags requiring plev_geopot. @@ -263,15 +266,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') + aviation_geopot_flag = diag_samp('aviation__geopot_thickness_850') .or. & + diag_samp('aviation__geopot_thickness_500') .or. & + diag_samp('aviation__snow_probability') plev_geopot_clim_flag = diag_samp('plev__geopot_clim') plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', & - activate=plev_geopot_clim_flag .or. aviation_thickness_flag) + activate=plev_geopot_clim_flag .or. aviation_geopot_flag) plev_gepot_required = plev_geopot_flag .or. & plev_geopot_clim_flag .or. & - aviation_thickness_flag + aviation_geopot_flag if (plev_gepot_required .and. use_xios_io) then height_w3 => get_height_fv(W3, exner_w3%get_mesh_id()) diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 new file mode 100644 index 000000000..d252d68fc --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -0,0 +1,139 @@ +!------------------------------------------------------------------------------- +! (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. +!------------------------------------------------------------------------------- +! +! 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_icao_heights_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 constants_mod, only: r_def, i_def, l_def + use kernel_mod, only: kernel_type + + implicit none + + ! the aviation diagnostics kernel type. + type, extends(kernel_type) :: aviation_icao_heights_kernel_type + type(arg_type), dimension(2) :: meta_args = (/ & + + ! output icao height. + arg_type(gh_field, gh_real, gh_write, any_discontinuous_space_1), & + + ! input pressure + arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_1) & + + /) + + integer :: operates_on = cell_column + + contains + procedure, nopass :: code => aviation_icao_heights_kernel_code + end type aviation_icao_heights_kernel_type + +contains + + subroutine aviation_icao_heights_kernel_code(nlayers, & + icao_height, & ! output icao height. + pressure_field, & ! pressure. + ndf, undf, map) + + ! calculate iaco height from the pressure field. + ! assumes lowest order w3 data, where ndf is always 1. + + use constants_mod, only: rmdi + use planet_config_mod, only: gravity, rd + use science_conversions_mod, only: feet_to_metres + use science_aviation_constants_mod, only: mtokft, & + isa_lapse_ratel, isa_lapse_rateu, isa_press_bot, & + isa_press_mid, isa_press_top, isa_temp_bot, isa_temp_top, & + gpm1, gpm2 + + 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) :: ndf + + ! number of unique degrees of freedom in the fields. + integer(kind=i_def), intent(in) :: undf + + ! degrees of freedom maps. offsets to the bottom of each column. + integer(kind=i_def), intent(in), dimension(ndf) :: map + + + ! arguments (algorithm) + + ! output icao height. + real(kind=r_def), intent(inout), dimension(undf) :: icao_height + + ! pressure in pa. + real(kind=r_def), intent(in), dimension(undf) :: pressure_field + + + ! local variables + integer(kind=i_def) :: df + + real(kind=r_def) :: g_over_rd + real(kind=r_def) :: zp1 + real(kind=r_def) :: zp2 + + real(kind=r_def) :: pressure + + + g_over_rd = gravity / rd + zp1 = isa_lapse_ratel / g_over_rd + zp2 = isa_lapse_rateu / g_over_rd + + pressure = pressure_field(map(1)) + + ! setting a safeguard limit to the lowest pressure to prevent + ! extremely large icao height values near the top at the + ! atmosphere. + if ( (pressure >= 0.0_r_def) .and. (pressure <= 1000.0_r_def) ) then + pressure = 1000.0_r_def + end if + + ! pressure must not be greater than surface pressure. + pressure = min(isa_press_bot, pressure) + + ! missing or invalid data? + if (pressure == rmdi .or. pressure <= 0.0_r_def) then + icao_height(map(1)) = rmdi + + ! heights up to 11,000 gpm + else if ( pressure > isa_press_mid ) then + pressure = pressure / isa_press_bot + pressure = 1.0_r_def - pressure**zp1 + icao_height(map(1)) = pressure * isa_temp_bot / isa_lapse_ratel + + ! heights between 11,000 gpm and 20,000 gpm + else if ( pressure > isa_press_top ) then + pressure = pressure / isa_press_mid + pressure = -log(pressure) + icao_height(map(1)) = gpm1 + pressure * isa_temp_top / g_over_rd + + ! heights above 20,000 gpm + else + pressure = pressure / isa_press_top + pressure = 1.0_r_def - pressure**zp2 + icao_height(map(1)) = gpm2 + pressure * isa_temp_top / isa_lapse_rateu + end if + + end subroutine aviation_icao_heights_kernel_code + +end module aviation_icao_heights_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 similarity index 75% rename from interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 rename to interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 index 9464a6ee0..fe5ea655b 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 @@ -5,43 +5,43 @@ !------------------------------------------------------------------------------- ! ! Section 20 aviation diagnostics kernel. +! This kernel calculates geopotential thickness and snow probabbility. ! ! 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 +module aviation_thickness_snow_kernel_mod - use argument_mod, only: arg_type, & + 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, & + 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, extends(kernel_type) :: aviation_thickness_snow_kernel_type type(arg_type), dimension(10) :: meta_args = (/ & - ! Output fields. + ! Output fields: thickness 850 & 500, and snow probability. 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), & - ! Source field. + ! Source field plev_geopot. arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_2), & - ! Request flags. + ! Field request flags. arg_type(gh_scalar, gh_logical, gh_read), & arg_type(gh_scalar, gh_logical, gh_read), & arg_type(gh_scalar, gh_logical, gh_read), & - ! Level indices. + ! Level indices: 1000, 850 & 500. arg_type(gh_scalar, gh_integer, gh_read), & arg_type(gh_scalar, gh_integer, gh_read), & arg_type(gh_scalar, gh_integer, gh_read) & @@ -50,26 +50,26 @@ module aviation_diags_kernel_mod integer :: operates_on = cell_column contains - procedure, nopass :: code => aviation_diags_kernel_code - end type aviation_diags_kernel_type + procedure, nopass :: code => aviation_thickness_snow_kernel_code + end type aviation_thickness_snow_kernel_type contains - subroutine aviation_diags_kernel_code(nlayers, & - ! Output fields. + subroutine aviation_thickness_snow_kernel_code(nlayers, & + ! output fields. thickness_850, thickness_500, snow_probability, & - ! Source field. - plev_geopot, & + ! source field. + plev_geopot, & - ! Request flags. + ! request flags. thickness_850_flag, thickness_500_flag, snow_probability_flag, & - ! Level indices. - i1000, i850, i500, & + ! level incides. + i1000, i850, i500, & - ! Kernel stuff. - result_ndf, result_undf, result_map, & + ! kernel stuff. + result_ndf, result_undf, result_map, & source_ndf, source_undf, source_map) ! Thickness: @@ -81,7 +81,7 @@ subroutine aviation_diags_kernel_code(nlayers, & implicit none - ! Arguments (kernel) + ! Arguments (kernel). ! The number of layers in a column. integer(kind=i_def), intent(in) :: nlayers @@ -92,14 +92,14 @@ subroutine aviation_diags_kernel_code(nlayers, & ! 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. + ! 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) + ! Arguments (algorithm). - ! Output fields. + ! 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 real(kind=r_def), intent(out), dimension(result_undf) :: snow_probability @@ -112,14 +112,15 @@ subroutine aviation_diags_kernel_code(nlayers, & logical(kind=l_def), intent(in) :: thickness_500_flag logical(kind=l_def), intent(in) :: snow_probability_flag - ! Level indices. + ! Level indices. for i850 and i500, -1 means "not requested". integer(kind=i_def), intent(in) :: i1000, i850, i500 - ! Local variables + ! Local variables. integer(kind=i_def) :: df real(kind=r_def) :: gph_1000, gph_850, gph_500 + ! Process every DOF in this cell. do df = 1, result_ndf @@ -144,6 +145,6 @@ subroutine aviation_diags_kernel_code(nlayers, & end do - end subroutine aviation_diags_kernel_code + end subroutine aviation_thickness_snow_kernel_code -end module aviation_diags_kernel_mod +end module aviation_thickness_snow_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 deleted file mode 100644 index 4489b53f8..000000000 --- a/interfaces/physics_schemes_interface/source/kernel/diags_geopot_kernel_mod.F90 +++ /dev/null @@ -1,133 +0,0 @@ -!------------------------------------------------------------------------------- -! (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 diags_geopot_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) :: diags_geopot_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 => diags_geopot_kernel_code - END TYPE diags_geopot_kernel_type - -CONTAINS - - SUBROUTINE diags_geopot_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 diags_geopot_kernel_code - -END MODULE diags_geopot_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 deleted file mode 100644 index ebcbc5ca1..000000000 --- a/interfaces/physics_schemes_interface/source/kernel/diags_icao_heights_kernel_mod.F90 +++ /dev/null @@ -1,139 +0,0 @@ -!------------------------------------------------------------------------------- -! (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 diags_icao_heights_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 constants_mod, ONLY: r_def, i_def, l_def - USE kernel_mod, ONLY: kernel_type - - IMPLICIT NONE - - ! The aviation diagnostics kernel type. - TYPE, EXTENDS(kernel_type) :: diags_icao_heights_kernel_type - TYPE(arg_type), DIMENSION(2) :: meta_args = (/ & - - ! Output icao height. - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & - - ! Input pressure - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1) & - - /) - - INTEGER :: operates_on = cell_column - - CONTAINS - PROCEDURE, NOPASS :: code => diags_icao_heights_kernel_code - END TYPE diags_icao_heights_kernel_type - -CONTAINS - - SUBROUTINE diags_icao_heights_kernel_code(nlayers, & - icao_height, & ! Output icao height. - pressure_field, & ! Pressure. - ndf, undf, map) - - ! Calculate iaco height from the pressure field. - ! Assumes lowest order W3 data, where ndf is always 1. - - USE constants_mod, ONLY: RMDI - USE planet_config_mod, ONLY: gravity, rd - USE science_conversions_mod, ONLY: feet_to_metres - USE science_aviation_constants_mod, ONLY: mtokft, & - isa_Lapse_RateL, isa_Lapse_RateU, isa_Press_Bot, & - isa_Press_Mid, isa_Press_Top, isa_Temp_Bot, isa_Temp_Top, & - Gpm1, Gpm2 - - 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) :: ndf - - ! Number of unique degrees of freedom in the fields. - INTEGER(KIND=i_def), INTENT(IN) :: undf - - ! Degrees of freedom maps. Offsets to the bottom of each column. - INTEGER(KIND=i_def), INTENT(IN), DIMENSION(ndf) :: map - - - ! Arguments (algorithm) - - ! Output icao height. - REAL(KIND=r_def), INTENT(INOUT), DIMENSION(undf) :: icao_height - - ! Pressure in Pa. - REAL(KIND=r_def), INTENT(IN), DIMENSION(undf) :: pressure_field - - - ! Local variables - INTEGER(KIND=i_def) :: df - - REAL(KIND=r_def) :: g_over_rd - REAL(KIND=r_def) :: zp1 - REAL(KIND=r_def) :: zp2 - - REAL(KIND=r_def) :: pressure - - - g_over_rd = gravity / rd - zp1 = isa_Lapse_RateL / g_over_rd - zp2 = isa_Lapse_RateU / g_over_rd - - pressure = pressure_field(map(1)) - - ! Setting a safeguard limit to the lowest pressure to prevent - ! extremely large icao height values near the top at the - ! atmosphere. - if ( (pressure >= 0.0_r_def) .and. (pressure <= 1000.0_r_def) ) then - pressure = 1000.0_r_def - end if - - ! Pressure must not be greater than surface pressure. - pressure = MIN(isa_Press_bot, pressure) - - ! Missing or invalid data? - if (pressure == RMDI .or. pressure <= 0.0_r_def) then - icao_height(map(1)) = RMDI - - ! Heights up to 11,000 GPM - else if ( pressure > isa_Press_Mid ) then - pressure = pressure / isa_Press_Bot - pressure = 1.0_r_def - pressure**zp1 - icao_height(map(1)) = pressure * isa_Temp_Bot / isa_Lapse_RateL - - ! Heights between 11,000 GPM and 20,000 GPM - else if ( pressure > isa_Press_Top ) then - pressure = pressure / isa_Press_Mid - pressure = -LOG(pressure) - icao_height(map(1)) = Gpm1 + pressure * isa_Temp_Top / g_over_rd - - ! Heights above 20,000 GPM - else - pressure = pressure / isa_Press_Top - pressure = 1.0_r_def - pressure**zp2 - icao_height(map(1)) = Gpm2 + pressure * isa_Temp_Top / isa_Lapse_RateU - end if - - END SUBROUTINE diags_icao_heights_kernel_code - -END MODULE diags_icao_heights_kernel_mod From a2e5556e63d45a185b127e28051ea63981521817 Mon Sep 17 00:00:00 2001 From: bblay Date: Tue, 12 May 2026 12:44:08 +0100 Subject: [PATCH 08/26] turn things off until it runs --- science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index 2f710b093..ac46c2175 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -342,7 +342,7 @@ subroutine gungho_diagnostics_driver( modeldb, & ! Wet bulb freezing level call freeze_lev_alg(theta, mr, moist_dyn, exner_in_wth) ! Aviation diagnostics - call aviation_diags_alg(plev_geopot, modeldb%fields%get_field_collection("convection_fields")) + !call aviation_diags_alg(plev_geopot, modeldb%fields%get_field_collection("convection_fields")) #endif temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') From 8f31833e09d70aeb4be34ec597452ebee8f626db Mon Sep 17 00:00:00 2001 From: bblay Date: Tue, 12 May 2026 15:51:02 +0100 Subject: [PATCH 09/26] put stuff back on --- science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index ac46c2175..2f710b093 100644 --- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 +++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 @@ -342,7 +342,7 @@ subroutine gungho_diagnostics_driver( modeldb, & ! Wet bulb freezing level call freeze_lev_alg(theta, mr, moist_dyn, exner_in_wth) ! Aviation diagnostics - !call aviation_diags_alg(plev_geopot, modeldb%fields%get_field_collection("convection_fields")) + call aviation_diags_alg(plev_geopot, modeldb%fields%get_field_collection("convection_fields")) #endif temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') From efb9da4b7760541dc2dde7d3a0e39e354378ce80 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 14 May 2026 09:51:07 +0100 Subject: [PATCH 10/26] icao heights gpm to kft --- applications/lfric_atm/metadata/field_def_diags.xml | 4 ++-- .../source/kernel/aviation_icao_heights_kernel_mod.F90 | 5 +++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 2e2f8a10f..9ea29a432 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -995,8 +995,8 @@ - - + + diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 index d252d68fc..8a420c7dc 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -134,6 +134,11 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & icao_height(map(1)) = gpm2 + pressure * isa_temp_top / isa_lapse_rateu end if + ! convert to kft + if (icao_height(map(1)) /= rmdi) then + icao_height(map(1)) = icao_height(map(1)) * mtokft + end if + end subroutine aviation_icao_heights_kernel_code end module aviation_icao_heights_kernel_mod From 68d0af0589a93cd049d45f19739af8b0d876d551 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 14 May 2026 11:24:57 +0100 Subject: [PATCH 11/26] name fix --- applications/lfric_atm/metadata/field_def_diags.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 9ea29a432..c1d193212 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -996,7 +996,7 @@ - + From 6c0d13269af75429d0af298fa1bb1f4012b89014 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 14 May 2026 15:26:53 +0100 Subject: [PATCH 12/26] finishing touches --- .../algorithm/aviation_diags_alg_mod.x90 | 100 +++++++++--------- .../aviation_icao_heights_kernel_mod.F90 | 53 +++++----- .../aviation_thickness_snow_kernel_mod.F90 | 9 +- .../science_aviation_constants_mod.F90 | 62 +++++------ 4 files changed, 109 insertions(+), 115 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index 2ea23df5c..245ca156b 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -61,8 +61,7 @@ contains ! 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, & - aviation_snow_probability_flag + logical(l_def) :: thick_850_flag, thick_500_flag, snow_probability_flag ! The array of pressure levels. integer(i_def) :: nplev @@ -76,8 +75,7 @@ contains real(r_def), parameter :: plev_tol = 0.1_r_def ! The output fields we produce. - type( field_type ) :: aviation_thick_850, aviation_thick_500, & - aviation_snow_probability + type( field_type ) :: thick_850, thick_500, snow_probability integer(i_def) :: i @@ -85,18 +83,17 @@ contains ! Check the request flags. - aviation_snow_probability_flag = & - init_diag(aviation_snow_probability, 'aviation__snow_probability') - aviation_thick_850_flag = & - init_diag(aviation_thick_850, 'aviation__geopot_thickness_850', & - activate=aviation_snow_probability_flag) - aviation_thick_500_flag = & - init_diag(aviation_thick_500, 'aviation__geopot_thickness_500') + snow_probability_flag = & + init_diag(snow_probability, 'aviation__snow_probability') + thick_850_flag = & + init_diag(thick_850, 'aviation__geopot_thickness_850') + thick_500_flag = & + init_diag(thick_500, 'aviation__geopot_thickness_500') ! Anything to do? - if ( .not. (aviation_thick_850_flag .or. & - aviation_thick_500_flag .or. & - aviation_snow_probability_flag) ) then + if ( .not. (thick_850_flag .or. & + thick_500_flag .or. & + snow_probability_flag) ) then return end if @@ -105,23 +102,23 @@ contains endif ! Initialise required fields to 1. - if ( aviation_snow_probability_flag ) then + if ( snow_probability_flag ) then message = 'Section 20: snow probability requested' call log_event( message, log_level_debug ) - call invoke( setval_c(aviation_snow_probability, 0.0_r_def)) + call invoke( setval_c(snow_probability, 0.0_r_def)) end if - if ( aviation_thick_850_flag ) then + if ( thick_850_flag ) then message = 'Section 20: geopot thick 850 requested' call log_event( message, log_level_debug ) ! Todo: Why 1 and not 0? - call invoke( setval_c(aviation_thick_850, 1.0_r_def)) + call invoke( setval_c(thick_850, 1.0_r_def)) end if - if ( aviation_thick_500_flag ) then + if ( thick_500_flag ) then message = 'Section 20: geopot thick 500 requested' call log_event( message, log_level_debug ) - call invoke( setval_c(aviation_thick_500, 1.0_r_def)) + call invoke( setval_c(thick_500, 1.0_r_def)) end if ! Get the array of pressure levels. @@ -169,13 +166,13 @@ contains return end if - if (aviation_thick_850_flag .and. i850 == -1) then + if ((thick_850_flag .or. snow_probability_flag) .and. i850 == -1) then message = 'Section 20: could not find 850hPa' call log_event( message, log_level_error ) return end if - if (aviation_thick_500_flag .and. i500 == -1) then + if (thick_500_flag .and. i500 == -1) then message = 'Section 20: could not find 500hPa' call log_event( message, log_level_error ) return @@ -183,17 +180,16 @@ contains ! 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( aviation_thickness_snow_kernel_type( & - aviation_thick_850, aviation_thick_500, & - aviation_snow_probability, & - plev_geopot, & - aviation_thick_850_flag, aviation_thick_500_flag, & - aviation_snow_probability_flag, i1000, i850, i500)) + call invoke( aviation_thickness_snow_kernel_type( & + thick_850, thick_500, snow_probability, & + plev_geopot, & + thick_850_flag, thick_500_flag, snow_probability_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() - if (aviation_snow_probability_flag) call aviation_snow_probability%write_field() + if (thick_850_flag) call thick_850%write_field() + if (thick_500_flag) call thick_500%write_field() + if (snow_probability_flag) call snow_probability%write_field() ! Clean up. if (allocated(plevs)) then @@ -203,13 +199,13 @@ contains end subroutine from_plev_geopot - SUBROUTINE icao_heights(convection_fields) - ! Calculate iaco heights of top and base of cloud from their pressures. + subroutine icao_heights(convection_fields) + ! Calculate ICAO heights of top and base of cloud from their pressures. - ! Arguments - TYPE(field_collection_type), INTENT(INOUT) :: convection_fields + ! Arguments. + type(field_collection_type), intent(inout) :: convection_fields - ! Locals + ! Locals. character(len=str_long) :: message ! Pressure fields to unpack. @@ -217,14 +213,14 @@ contains 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 + 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 + logical(l_def) :: base_height_flag + logical(l_def) :: top_height_flag - ! which heights are requested this time step? + ! Which heights are requested this time step? top_height_flag = init_diag( & conv_cloud_top_icao_height, 'aviation__conv_cloud_top_icao_height') @@ -232,20 +228,20 @@ contains conv_cloud_base_icao_height, 'aviation__conv_cloud_base_icao_height') if (base_height_flag .or. top_height_flag) then - message = 'Section 20: cloud ICAO heights on' + message = 'section 20: cloud icao heights on' call log_event( message, log_level_debug ) else - ! nothing requested at this time step + ! Nothing requested at this time step. return end if - ! cloud base - IF (base_height_flag) THEN + ! 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), & + setval_c(conv_cloud_base_icao_height, rmdi), & aviation_icao_heights_kernel_type(conv_cloud_base_icao_height, & conv_cloud_base_pressure)) @@ -256,14 +252,14 @@ contains ! Based on https://code.metoffice.gov.uk/trac/lfric_apps/ticket/189 ! call convection_fields%remove_field('convection__pres_cv_base') - END IF + end if - ! could top - IF (top_height_flag) THEN + ! Cloud 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), & + setval_c(conv_cloud_top_icao_height, rmdi), & aviation_icao_heights_kernel_type(conv_cloud_top_icao_height, & conv_cloud_top_pressure)) @@ -272,8 +268,8 @@ contains ! Is it necessary/safe to remove this field now? ! call convection_fields%remove_field('convection__pres_cv_top') - END IF + end if - END SUBROUTINE icao_heights + end subroutine icao_heights end module aviation_diags_alg_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 index 8a420c7dc..e5a34ac03 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -23,14 +23,14 @@ module aviation_icao_heights_kernel_mod implicit none - ! the aviation diagnostics kernel type. + ! The aviation diagnostics kernel type. type, extends(kernel_type) :: aviation_icao_heights_kernel_type type(arg_type), dimension(2) :: meta_args = (/ & - ! output icao height. + ! Output icao height. arg_type(gh_field, gh_real, gh_write, any_discontinuous_space_1), & - ! input pressure + ! Input pressure arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_1) & /) @@ -44,48 +44,47 @@ module aviation_icao_heights_kernel_mod contains subroutine aviation_icao_heights_kernel_code(nlayers, & - icao_height, & ! output icao height. - pressure_field, & ! pressure. + icao_height, & ! Output icao height. + pressure_field, & ! Pressure. ndf, undf, map) - ! calculate iaco height from the pressure field. - ! assumes lowest order w3 data, where ndf is always 1. + ! Calculate iaco height from the pressure field. + ! Assumes lowest order w3 data, where ndf is always 1. use constants_mod, only: rmdi use planet_config_mod, only: gravity, rd - use science_conversions_mod, only: feet_to_metres - use science_aviation_constants_mod, only: mtokft, & - isa_lapse_ratel, isa_lapse_rateu, isa_press_bot, & + use science_aviation_constants_mod, only: mtokft, & + isa_lapse_ratel, isa_lapse_rateu, isa_press_bot, & isa_press_mid, isa_press_top, isa_temp_bot, isa_temp_top, & gpm1, gpm2 implicit none - ! arguments (kernel) + ! Arguments (kernel). - ! the number of layers in a column. + ! 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. + ! Number of degrees of freedom (columns) in the cell we're processing. integer(kind=i_def), intent(in) :: ndf - ! number of unique degrees of freedom in the fields. + ! Number of unique degrees of freedom in the fields. integer(kind=i_def), intent(in) :: undf - ! degrees of freedom maps. offsets to the bottom of each column. + ! Degrees of freedom maps. offsets to the bottom of each column. integer(kind=i_def), intent(in), dimension(ndf) :: map - ! arguments (algorithm) + ! Arguments (algorithm) - ! output icao height. - real(kind=r_def), intent(inout), dimension(undf) :: icao_height + ! Output icao height. + real(kind=r_def), intent(out), dimension(undf) :: icao_height - ! pressure in pa. + ! Pressure in pa. real(kind=r_def), intent(in), dimension(undf) :: pressure_field - ! local variables + ! Local variables integer(kind=i_def) :: df real(kind=r_def) :: g_over_rd @@ -101,40 +100,40 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & pressure = pressure_field(map(1)) - ! setting a safeguard limit to the lowest pressure to prevent + ! Setting a safeguard limit to the lowest pressure to prevent ! extremely large icao height values near the top at the ! atmosphere. if ( (pressure >= 0.0_r_def) .and. (pressure <= 1000.0_r_def) ) then pressure = 1000.0_r_def end if - ! pressure must not be greater than surface pressure. + ! Pressure must not be greater than surface pressure. pressure = min(isa_press_bot, pressure) - ! missing or invalid data? + ! Missing or invalid data? if (pressure == rmdi .or. pressure <= 0.0_r_def) then icao_height(map(1)) = rmdi - ! heights up to 11,000 gpm + ! Heights up to 11,000 gpm else if ( pressure > isa_press_mid ) then pressure = pressure / isa_press_bot pressure = 1.0_r_def - pressure**zp1 icao_height(map(1)) = pressure * isa_temp_bot / isa_lapse_ratel - ! heights between 11,000 gpm and 20,000 gpm + ! Heights between 11,000 gpm and 20,000 gpm else if ( pressure > isa_press_top ) then pressure = pressure / isa_press_mid pressure = -log(pressure) icao_height(map(1)) = gpm1 + pressure * isa_temp_top / g_over_rd - ! heights above 20,000 gpm + ! Heights above 20,000 gpm else pressure = pressure / isa_press_top pressure = 1.0_r_def - pressure**zp2 icao_height(map(1)) = gpm2 + pressure * isa_temp_top / isa_lapse_rateu end if - ! convert to kft + ! Convert to kft if (icao_height(map(1)) /= rmdi) then icao_height(map(1)) = icao_height(map(1)) * mtokft end if diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 index fe5ea655b..d8eef3734 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 @@ -118,15 +118,13 @@ subroutine aviation_thickness_snow_kernel_code(nlayers, & ! Local variables. integer(kind=i_def) :: df - real(kind=r_def) :: gph_1000, gph_850, gph_500 + real(kind=r_def) :: gph_1000, gph_850 ! Process every DOF in this cell. do df = 1, result_ndf gph_1000 = plev_geopot(source_map(df) + i1000-1) - gph_850 = plev_geopot(source_map(df) + i850-1) - gph_500 = plev_geopot(source_map(df) + i500-1) if (thickness_850_flag .and. i850 /= -1) then thickness_850(result_map(df)) = & @@ -138,9 +136,10 @@ subroutine aviation_thickness_snow_kernel_code(nlayers, & plev_geopot(source_map(df)+i500-1) - gph_1000 end if - if(snow_probability_flag) then + if(snow_probability_flag .and. i850 /= -1) then + gph_850 = plev_geopot(source_map(df) + i850-1) snow_probability(result_map(df)) = & - 5220.0 + 3.86666*gph_1000 - 4.0*gph_850 + 5220.0_r_def + 3.86666_r_def*gph_1000 - 4.0_r_def*gph_850 end if end do diff --git a/science/shared/source/constants/science_aviation_constants_mod.F90 b/science/shared/source/constants/science_aviation_constants_mod.F90 index 0176d9878..bddda2365 100644 --- a/science/shared/source/constants/science_aviation_constants_mod.F90 +++ b/science/shared/source/constants/science_aviation_constants_mod.F90 @@ -6,50 +6,50 @@ !> @brief LFRic aviation constants module !---------------------------------------------------------------------------- ! -MODULE science_aviation_constants_mod +module science_aviation_constants_mod - USE constants_mod, ONLY: r_def - USE science_conversions_mod, ONLY: feet_to_metres + use constants_mod, only: r_def + use science_conversions_mod, only: feet_to_metres - IMPLICIT NONE + implicit none - PRIVATE - PUBLIC :: mtokft, isa_Lapse_RateL, isa_Lapse_RateU, & - isa_Press_Bot, isa_Press_Mid, isa_Press_Top, & - isa_Temp_Bot, isa_Temp_Top, Gpm1, Gpm2 + private + public :: mtokft, isa_lapse_ratel, isa_lapse_rateu, & + isa_press_bot, isa_press_mid, isa_press_top, & + isa_temp_bot, isa_temp_top, gpm1, gpm2 - ! meters to thousands of feet - REAL(r_def), PARAMETER :: mtokft = 1.0_r_def / (feet_to_metres * 1000.0_r_def) + ! Meters to thousands of feet. + real(r_def), parameter :: mtokft = 1.0_r_def / (feet_to_metres * 1000.0_r_def) - ! International Standard Atmosphere (isa) within troposhpere for: + ! International Standard Atmosphere (ISA) within troposphere for: ! Values taken from: - ! https://en.wikipedia.org/wiki/International_Standard_Atmosphere - ! Lapse rate (degreeC/km) for levels below 11,000 gpm - REAL(r_def), PARAMETER :: isa_Lapse_RateL = 6.5e-03_r_def + ! https://en.wikipedia.org/wiki/international_standard_atmosphere + ! Lapse rate (degreeC/km) for levels below 11,000 gpm. + real(r_def), parameter :: isa_lapse_ratel = 6.5e-03_r_def - ! Lapse rate (degreeC/km) for levels above 11,000 gpm - REAL(r_def), PARAMETER :: isa_Lapse_RateU = -1.0e-03_r_def + ! Lapse rate (degreeC/km) for levels above 20,000 gpm. + real(r_def), parameter :: isa_lapse_rateu = -1.0e-03_r_def - ! surface pressure (Pa) - REAL(r_def), PARAMETER :: isa_Press_Bot = 101325.0_r_def + ! Surface pressure (Pa). + real(r_def), parameter :: isa_press_bot = 101325.0_r_def - ! pressure (Pa) at 11,000 gpm - REAL(r_def), PARAMETER :: isa_Press_Mid = 22632.0_r_def + ! Pressure (Pa) at 11,000 gpm. + real(r_def), parameter :: isa_press_mid = 22632.0_r_def - ! pressure (Pa) at 20,000 gpm - REAL(r_def), PARAMETER :: isa_Press_Top = 5474.87_r_def + ! Pressure (Pa) at 20,000 gpm. + real(r_def), parameter :: isa_press_top = 5474.87_r_def - ! Surface temperature (K) - REAL(r_def), PARAMETER :: isa_Temp_Bot = 288.15_r_def + ! Surface temperature (K). + real(r_def), parameter :: isa_temp_bot = 288.15_r_def - ! Temperature (K) of isothermal layer - at tropopause - REAL(r_def), PARAMETER :: isa_Temp_Top = 216.65_r_def + ! Temperature (K) of isothermal layer - at tropopause. + real(r_def), parameter :: isa_temp_top = 216.65_r_def - ! Height limit (gpm) for std lower lapse rate - REAL(r_def), PARAMETER :: Gpm1 = 11000.0_r_def + ! Height limit (gpm) for standard lower lapse rate. + real(r_def), parameter :: gpm1 = 11000.0_r_def - ! Height (gpm) of top of isothermal layer - REAL(r_def), PARAMETER :: Gpm2 = 20000.0_r_def + ! Height (gpm) of top of isothermal layer. + real(r_def), parameter :: gpm2 = 20000.0_r_def -END MODULE science_aviation_constants_mod +end module science_aviation_constants_mod From 12b13dd03b03fd1b70d3cb6ba848c462fef29d60 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 14 May 2026 15:40:47 +0100 Subject: [PATCH 13/26] finishing touches --- .../source/algorithm/aviation_diags_alg_mod.x90 | 1 - 1 file changed, 1 deletion(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index 245ca156b..636b51936 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -249,7 +249,6 @@ contains ! 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 From 8372788f42f5a5ae03bd548a457443449cf28e25 Mon Sep 17 00:00:00 2001 From: bblay Date: Tue, 19 May 2026 09:41:28 +0100 Subject: [PATCH 14/26] snow prob clamp --- .../kernel/aviation_thickness_snow_kernel_mod.F90 | 12 ++++++++++++ .../lfric_atm/file/file_def_diags_oper_nwp_gl.xml | 2 +- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 index d8eef3734..6beb01af4 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 @@ -140,6 +140,18 @@ subroutine aviation_thickness_snow_kernel_code(nlayers, & gph_850 = plev_geopot(source_map(df) + i850-1) snow_probability(result_map(df)) = & 5220.0_r_def + 3.86666_r_def*gph_1000 - 4.0_r_def*gph_850 + + if (gph_1000 > 1.0e8) then + snow_probability(result_map(df)) = 0.0_r_def + end if + + ! Limit to percentage. + if (snow_probability(result_map(df)) < 0.0_r_def) then + snow_probability(result_map(df)) = 0.0_r_def + else if (snow_probability(result_map(df)) > 100.0_r_def) then + snow_probability(result_map(df)) = 100.0_r_def + end if + end if end do diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml index 8d70bf348..5318da9e6 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml @@ -193,7 +193,7 @@ - + From f834cdfa4fad7b61d5c0fdd6f364f6d9c706ad15 Mon Sep 17 00:00:00 2001 From: bblay Date: Tue, 9 Jun 2026 15:56:25 +0100 Subject: [PATCH 15/26] review actions wip --- .../aviation_icao_heights_kernel_mod.F90 | 39 +++--- .../aviation_thickness_snow_kernel_mod.F90 | 47 ++++---- .../source/kernel/conv_comorph_kernel_mod.F90 | 8 +- .../source/kernel/conv_gr_kernel_mod.F90 | 8 +- .../aviation_icao_heights_kernel_mod_test.pf | 114 ++++++++++++++++++ 5 files changed, 162 insertions(+), 54 deletions(-) create mode 100644 interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 index e5a34ac03..b618ab2ea 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -43,16 +43,20 @@ module aviation_icao_heights_kernel_mod contains + !> @brief Calculate icao height from the pressure field. + !> Assumes lowerst order W3 data, where ndf is always 1. + !> @param[in] nlayers The number of layers in a column. + !> @param[out] icao_height Output icao height in kft. + !> @param[in] pressure_field Pressure in pa. + !> @param[in] ndf Number of DOFs in the cell. + !> @param[in] undf Number of DOFs in the field. + !> @param[in] map Dofmap to the bottom cell. subroutine aviation_icao_heights_kernel_code(nlayers, & - icao_height, & ! Output icao height. - pressure_field, & ! Pressure. + icao_height, pressure_field, & ndf, undf, map) - ! Calculate iaco height from the pressure field. - ! Assumes lowest order w3 data, where ndf is always 1. - use constants_mod, only: rmdi - use planet_config_mod, only: gravity, rd + use planet_config_mod, only: g_over_r use science_aviation_constants_mod, only: mtokft, & isa_lapse_ratel, isa_lapse_rateu, isa_press_bot, & isa_press_mid, isa_press_top, isa_temp_bot, isa_temp_top, & @@ -61,42 +65,27 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & 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) :: ndf - - ! Number of unique degrees of freedom in the fields. integer(kind=i_def), intent(in) :: undf - - ! Degrees of freedom maps. offsets to the bottom of each column. integer(kind=i_def), intent(in), dimension(ndf) :: map - ! Arguments (algorithm) - - ! Output icao height. real(kind=r_def), intent(out), dimension(undf) :: icao_height - - ! Pressure in pa. real(kind=r_def), intent(in), dimension(undf) :: pressure_field ! Local variables integer(kind=i_def) :: df - real(kind=r_def) :: g_over_rd real(kind=r_def) :: zp1 real(kind=r_def) :: zp2 real(kind=r_def) :: pressure - g_over_rd = gravity / rd - zp1 = isa_lapse_ratel / g_over_rd - zp2 = isa_lapse_rateu / g_over_rd + zp1 = isa_lapse_ratel / g_over_r + zp2 = isa_lapse_rateu / g_over_r pressure = pressure_field(map(1)) @@ -111,7 +100,7 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & pressure = min(isa_press_bot, pressure) ! Missing or invalid data? - if (pressure == rmdi .or. pressure <= 0.0_r_def) then + if (pressure == rmdi) then icao_height(map(1)) = rmdi ! Heights up to 11,000 gpm @@ -124,7 +113,7 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & else if ( pressure > isa_press_top ) then pressure = pressure / isa_press_mid pressure = -log(pressure) - icao_height(map(1)) = gpm1 + pressure * isa_temp_top / g_over_rd + icao_height(map(1)) = gpm1 + pressure * isa_temp_top / g_over_r ! Heights above 20,000 gpm else diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 index 6beb01af4..e2e0d50d3 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_thickness_snow_kernel_mod.F90 @@ -55,6 +55,31 @@ module aviation_thickness_snow_kernel_mod contains + + !> @brief Calculate thickness and snow probability from the geopotential height at pressure levels. + !> Assumes lowest order W3 data, where ndf is always 1. + !> Thickness: + !> Subtract geopotential heights at 850 and 500 hPa from that at 1000 hPa. + !> Snow probability: + !> Implement Boyden (1964), using 850 and 1000 hPa. + !> See https://github.com/MetOffice/Section20/issues/21 + !> @param[in] nlayers The number of layers in a column. + !> @param[out] thickness_850 Output geopot thickness between 1000 and 850 hPa in m. + !> @param[out] thickness_500 Output geopot thickness between 1000 and 500 hPa in m. + !> @param[out] snow_probability Output snow probability in percentage. + !> @param[in] plev_geopot Geopotential height at pressure levels in m. + !> @param[in] thickness_850_flag thickness_850 request flag. + !> @param[in] thickness_500_flag thickness_500 request flag. + !> @param[in] snow_probability_flag snow_probability request flag. + !> @param[in] i1000 Index of 1000 hPa level in plev_geopot. + !> @param[in] i850 Index of 850 hPa level in plev_geopot. + !> @param[in] i500 Index of 500 hPa level in plev_geopot. + !> @param[in] result_ndf Number of DOFs in the result cell. + !> @param[in] result_undf Number of DOFs in the result field. + !> @param[in] result_map Dofmap to the result's bottom cell. + !> @param[in] source_ndf Number of DOFs in the source cell. + !> @param[in] source_undf Number of DOFs in the source field. + !> @param[in] source_map Dofmap to the source's bottom cell. subroutine aviation_thickness_snow_kernel_code(nlayers, & ! output fields. thickness_850, thickness_500, snow_probability, & @@ -72,39 +97,19 @@ subroutine aviation_thickness_snow_kernel_code(nlayers, & result_ndf, result_undf, result_map, & source_ndf, source_undf, source_map) - ! Thickness: - ! Subtract geopotential heights at 850 and 500 hPa from that at 1000 hPa. - ! - ! Snow probability: - ! Implement Boyden (1964), using 850 and 1000 hPa. - ! See https://github.com/MetOffice/Section20/issues/21 - 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 real(kind=r_def), intent(out), dimension(result_undf) :: snow_probability - - ! Geopotential height at pressure levels. real(kind=r_def), intent(in), dimension(source_undf) :: plev_geopot ! Request flags. @@ -112,7 +117,7 @@ subroutine aviation_thickness_snow_kernel_code(nlayers, & logical(kind=l_def), intent(in) :: thickness_500_flag logical(kind=l_def), intent(in) :: snow_probability_flag - ! Level indices. for i850 and i500, -1 means "not requested". + ! Level indices. integer(kind=i_def), intent(in) :: i1000, i850, i500 diff --git a/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 index ba4cb4bac..b8ad67f5f 100644 --- a/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 @@ -2716,7 +2716,7 @@ subroutine conv_comorph_code(nlayers, seg_len, & if (cct(i,1) > 0) then pres_cv_top(map_2d(1,i)) = p_rho_levels(i,1,cct(i,1)) else - pres_cv_top(map_2d(1,i)) = 0.0_r_def + pres_cv_top(map_2d(1,i)) = rmdi end if end do end if @@ -2725,7 +2725,7 @@ subroutine conv_comorph_code(nlayers, seg_len, & if (ccb(i,1) > 0) then pres_cv_base(map_2d(1,i)) = p_rho_levels(i,1,ccb(i,1)) else - pres_cv_base(map_2d(1,i))= 0.0_r_def + pres_cv_base(map_2d(1,i))= rmdi end if end do end if @@ -2736,7 +2736,7 @@ subroutine conv_comorph_code(nlayers, seg_len, & if (lctop(i,1) > 0) then pres_lowest_cv_top(map_2d(1,i)) = p_rho_levels(i,1,lctop(i,1)) else - pres_lowest_cv_top(map_2d(1,i)) = 0.0_r_def + pres_lowest_cv_top(map_2d(1,i)) = rmdi end if end do end if @@ -2745,7 +2745,7 @@ subroutine conv_comorph_code(nlayers, seg_len, & if (lcbase(i,1) > 0) then pres_lowest_cv_base(map_2d(1,i)) = p_rho_levels(i,1,lcbase(i,1)) else - pres_lowest_cv_base(map_2d(1,i))= 0.0_r_def + pres_lowest_cv_base(map_2d(1,i))= rmdi end if end do end if diff --git a/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 index 2fcf1d420..1b73cfb93 100644 --- a/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 @@ -2728,7 +2728,7 @@ subroutine conv_gr_code(nlayers, & if (cct(i,1) > 0) then pres_cv_top(map_2d(1,i)) = p_rho_levels(i,1,cct(i,1)) else - pres_cv_top(map_2d(1,i)) = 0.0_r_def + pres_cv_top(map_2d(1,i)) = rmdi end if end do end if @@ -2737,7 +2737,7 @@ subroutine conv_gr_code(nlayers, & if (ccb(i,1) > 0) then pres_cv_base(map_2d(1,i)) = p_rho_levels(i,1,ccb(i,1)) else - pres_cv_base(map_2d(1,i))= 0.0_r_def + pres_cv_base(map_2d(1,i))= rmdi end if end do end if @@ -2748,7 +2748,7 @@ subroutine conv_gr_code(nlayers, & if (lctop(i,1) > 0) then pres_lowest_cv_top(map_2d(1,i)) = p_rho_levels(i,1,lctop(i,1)) else - pres_lowest_cv_top(map_2d(1,i)) = 0.0_r_def + pres_lowest_cv_top(map_2d(1,i)) = rmdi end if end do end if @@ -2757,7 +2757,7 @@ subroutine conv_gr_code(nlayers, & if (lcbase(i,1) > 0) then pres_lowest_cv_base(map_2d(1,i)) = p_rho_levels(i,1,lcbase(i,1)) else - pres_lowest_cv_base(map_2d(1,i))= 0.0_r_def + pres_lowest_cv_base(map_2d(1,i))= rmdi end if end do end if diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf new file mode 100644 index 000000000..dba100df4 --- /dev/null +++ b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf @@ -0,0 +1,114 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- +!> Test aviation_icao_heights_kernel +!> + +module aviation_icao_heights_kernel_mod_test + + use funit + use constants_mod, only: rmdi + use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_kernel_code + + implicit none + + private + public :: test_rmdi, test_low_safeguard + + @TestCase + type, extends(TestCase), public :: aviation_icao_heights_test_type + private + real(kind=r_def), allocatable :: pressure(:) + real(kind=r_def), allocatable :: icao_heights(:) + real(kind=r_def), allocatable :: expect(:) + + contains + procedure setUp + procedure tearDown + procedure test_rmdi + procedure test_low_safeguard + + end type aviation_icao_heights_test_type + + integer(kind=i_def), parameter :: nlayers = -999_i_def + integer(kind=i_def), parameter :: ndf = 1_i_def + integer(kind=i_def), parameter :: undf = 2_i_def + integer(kind=i_def), parameter :: map(2) = (/1_i_def, -999_i_def /) + + integer(kind=i_def), parameter :: arr_len = undf + +! integer(kind=i_def), parameter :: mode_dimen = n_radaer_mode + 1_i_def +! integer(kind=i_def), parameter :: mmr_arr_len = (nlayers + 1) * (n_radaer_mode + 1) + + +contains + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @before + subroutine setUp( this ) + + implicit none + + class(aviation_icao_heights_test_type), intent(inout) :: this + + allocate(this%icao_heights(arr_len), & + this%pressure(arr_len), & + this%expect(arr_len) ) + + end subroutine setUp + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @after + subroutine tearDown( this ) + + implicit none + + class(aviation_icao_heights_test_type), intent(inout) :: this + + deallocate(this%icao_heights, this%pressure, this%expect) + + end subroutine tearDown + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @test + subroutine test_rmdi ( this ) + use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_code + implicit none + + class(aviation_icao_heights_test_type), intent(inout) :: this + + this%pressure = (/ rmdi, -999.0_r_def /) + this%expect = (/ rmdi, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_code( nlayers, & + this%icao_heights, this%pressure, & + ndf, undf, map ) + + @assertEqual( this%icao_heights, this%expect ) + + end subroutine test_rmdi + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @test + subroutine test_low_safeguard ( this ) + use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_code + implicit none + + class(aviation_icao_heights_test_type), intent(inout) :: this + + this%pressure = (/ 500.0_r_def, -999.0_r_def /) + this%expect = (/ 1234.0_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_code( nlayers, & + this%icao_heights, this%pressure, & + ndf, undf, map ) + + @assertEqual( this%icao_heights, this%expect ) + + end subroutine test_low_safeguard + +end module aviation_icao_heights_kernel_mod_test From 4760ec75c1f674ea33049eab0afb4a398ebbb797 Mon Sep 17 00:00:00 2001 From: bblay Date: Tue, 9 Jun 2026 16:39:13 +0100 Subject: [PATCH 16/26] review actions wip --- .../source/constants/planet_constants_mod.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/constants/planet_constants_mod.F90 b/interfaces/physics_schemes_interface/source/constants/planet_constants_mod.F90 index 707ad3367..7f80894b2 100644 --- a/interfaces/physics_schemes_interface/source/constants/planet_constants_mod.F90 +++ b/interfaces/physics_schemes_interface/source/constants/planet_constants_mod.F90 @@ -26,7 +26,7 @@ module planet_constants_mod c_virtual_bl, etar_bl, repsilon_bl, ls_bl, r_32b, c_virtual_32b, & etar_32b, lcrcp_32b, ls_32b, lsrcp_32b, planet_radius_bl, & recip_kappa_bl, power, ex_power, lcrcp_def, lsrcp_def, & - recip_kappa_def + recip_kappa_def, g_over_r ! The following variables have been hidden as they are not currently ! required to build the extracted UM code. They have been left in @@ -36,7 +36,6 @@ module planet_constants_mod ! Disabled variables: ! sclht, omega, two_omega, recip_p_zero, - ! g_over_r !---------------------------------------------------------------------- From eb2ad4c18a0affd8e653e4c630162aec1a6935e2 Mon Sep 17 00:00:00 2001 From: bblay Date: Wed, 10 Jun 2026 09:10:19 +0100 Subject: [PATCH 17/26] gotta commit every little change just to run cylc --- .../source/kernel/aviation_icao_heights_kernel_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 index b618ab2ea..89c205a61 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -56,7 +56,7 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & ndf, undf, map) use constants_mod, only: rmdi - use planet_config_mod, only: g_over_r + use planet_constants_mod, only: g_over_r use science_aviation_constants_mod, only: mtokft, & isa_lapse_ratel, isa_lapse_rateu, isa_press_bot, & isa_press_mid, isa_press_top, isa_temp_bot, isa_temp_top, & From a4efd9c6680b4ae593503d8d6e4a09e9867fe0e8 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 11 Jun 2026 10:21:33 +0100 Subject: [PATCH 18/26] trying to get unit tests to compile --- dependencies.yaml | 16 ++++++++-------- .../source/algorithm/conv_comorph_alg_mod.x90 | 2 +- .../source/kernel/conv_comorph_kernel_mod.F90 | 2 +- .../source/kernel/conv_gr_kernel_mod.F90 | 2 +- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/dependencies.yaml b/dependencies.yaml index 9bd6b263a..a6432034b 100644 --- a/dependencies.yaml +++ b/dependencies.yaml @@ -19,11 +19,11 @@ casim: source: git@github.com:MetOffice/casim.git - ref: 2026.03.2 + ref: stable jules: source: git@github.com:MetOffice/jules.git - ref: 2026.03.2 + ref: stable lfric_apps: source: @@ -31,24 +31,24 @@ lfric_apps: lfric_core: source: git@github.com:MetOffice/lfric_core.git - ref: 2026.03.2 + ref: stable moci: source: git@github.com:MetOffice/moci.git - ref: 2026.03.2 + ref: stable SimSys_Scripts: source: git@github.com:MetOffice/SimSys_Scripts.git - ref: 2026.03.2 + ref: stable socrates: source: git@github.com:MetOffice/socrates.git - ref: 2026.03.2 + ref: stable socrates-spectral: source: git@github.com:MetOffice/socrates-spectral.git - ref: 2026.03.2 + ref: stable ukca: source: git@github.com:MetOffice/ukca.git - ref: 2026.03.2 + ref: stable diff --git a/interfaces/physics_schemes_interface/source/algorithm/conv_comorph_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/conv_comorph_alg_mod.x90 index 466e96012..ba86b561d 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/conv_comorph_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/conv_comorph_alg_mod.x90 @@ -7,7 +7,7 @@ module conv_comorph_alg_mod - use constants_mod, only: i_def + use constants_mod, only: i_def, rmdi use field_mod, only: field_type use integer_field_mod, only: integer_field_type use field_collection_mod, only: field_collection_type diff --git a/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 index b8ad67f5f..75f19f7cf 100644 --- a/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/conv_comorph_kernel_mod.F90 @@ -15,7 +15,7 @@ module conv_comorph_kernel_mod ANY_DISCONTINUOUS_SPACE_1, & ANY_DISCONTINUOUS_SPACE_2, & ANY_DISCONTINUOUS_SPACE_3 - use constants_mod, only : i_def, i_um, r_def, r_um + use constants_mod, only : i_def, i_um, r_def, r_um, rmdi use empty_data_mod, only : empty_real_data use fs_continuity_mod, only : W3, Wtheta use kernel_mod, only : kernel_type diff --git a/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 index 1b73cfb93..6d80c4779 100644 --- a/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/conv_gr_kernel_mod.F90 @@ -14,7 +14,7 @@ module conv_gr_kernel_mod GH_READWRITE, DOMAIN, & ANY_DISCONTINUOUS_SPACE_1, & ANY_DISCONTINUOUS_SPACE_2 - use constants_mod, only : i_def, i_um, r_def, r_um + use constants_mod, only : i_def, i_um, r_def, r_um, rmdi use tuning_segments_mod, only : conv_gr_segment_size use empty_data_mod, only : empty_real_data use fs_continuity_mod, only : W3, Wtheta From fa8138b543ae9d1906cee146fd7a51a4655b59aa Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 11 Jun 2026 11:13:46 +0100 Subject: [PATCH 19/26] omg i keep forgetting to commit before i run, so inconvenient --- .../aviation_icao_heights_kernel_mod_test.pf | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf index dba100df4..4147e5b30 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf @@ -10,7 +10,7 @@ module aviation_icao_heights_kernel_mod_test use funit use constants_mod, only: rmdi - use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_kernel_code +! use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_kernel_code implicit none @@ -83,32 +83,32 @@ contains this%expect = (/ rmdi, -999.0_r_def /) this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) - call aviation_icao_heights_code( nlayers, & - this%icao_heights, this%pressure, & - ndf, undf, map ) +! call aviation_icao_heights_code( nlayers, & +! this%icao_heights, this%pressure, & +! ndf, undf, map ) @assertEqual( this%icao_heights, this%expect ) end subroutine test_rmdi !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @test - subroutine test_low_safeguard ( this ) - use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_code - implicit none - - class(aviation_icao_heights_test_type), intent(inout) :: this - - this%pressure = (/ 500.0_r_def, -999.0_r_def /) - this%expect = (/ 1234.0_r_def, -999.0_r_def /) - - this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) - call aviation_icao_heights_code( nlayers, & - this%icao_heights, this%pressure, & - ndf, undf, map ) - - @assertEqual( this%icao_heights, this%expect ) - - end subroutine test_low_safeguard +! @test +! subroutine test_low_safeguard ( this ) +! use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_code +! implicit none +! +! class(aviation_icao_heights_test_type), intent(inout) :: this +! +! this%pressure = (/ 500.0_r_def, -999.0_r_def /) +! this%expect = (/ 1234.0_r_def, -999.0_r_def /) +! +! this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) +! call aviation_icao_heights_code( nlayers, & +! this%icao_heights, this%pressure, & +! ndf, undf, map ) +! +! @assertEqual( this%icao_heights, this%expect ) +! +! end subroutine test_low_safeguard end module aviation_icao_heights_kernel_mod_test From ee719fbdbcb63d42a38c373fb1610a7172a550b2 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 11 Jun 2026 15:05:00 +0100 Subject: [PATCH 20/26] omg i keep forgetting to commit before i run, so inconvenient --- dependencies.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dependencies.yaml b/dependencies.yaml index a6432034b..01d8565ab 100644 --- a/dependencies.yaml +++ b/dependencies.yaml @@ -39,7 +39,7 @@ moci: SimSys_Scripts: source: git@github.com:MetOffice/SimSys_Scripts.git - ref: stable + ref: 2026.03.2 socrates: source: git@github.com:MetOffice/socrates.git From 20cb5881d061417012509755cadc9f2f125a093a Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 11 Jun 2026 15:14:05 +0100 Subject: [PATCH 21/26] omg i keep forgetting to commit before i run, so inconvenient --- dependencies.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dependencies.yaml b/dependencies.yaml index 01d8565ab..e3a28c9f9 100644 --- a/dependencies.yaml +++ b/dependencies.yaml @@ -47,7 +47,7 @@ socrates: socrates-spectral: source: git@github.com:MetOffice/socrates-spectral.git - ref: stable + ref: 2026.03.2 ukca: source: git@github.com:MetOffice/ukca.git From 8ebd1767048e8666b07eb35737895ba002d3b6c8 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 11 Jun 2026 18:37:07 +0100 Subject: [PATCH 22/26] can't get the unit test to compile --- .../aviation_icao_heights_kernel_mod_test.pf | 22 +++++++++++-------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf index 4147e5b30..c78c28302 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf @@ -9,13 +9,13 @@ module aviation_icao_heights_kernel_mod_test use funit - use constants_mod, only: rmdi -! use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_kernel_code + use constants_mod, only: i_def, r_def, rmdi +! use feign_config_mod, only: feign_extrusion_config implicit none private - public :: test_rmdi, test_low_safeguard + public :: test_rmdi @TestCase type, extends(TestCase), public :: aviation_icao_heights_test_type @@ -28,7 +28,6 @@ module aviation_icao_heights_kernel_mod_test procedure setUp procedure tearDown procedure test_rmdi - procedure test_low_safeguard end type aviation_icao_heights_test_type @@ -39,10 +38,6 @@ module aviation_icao_heights_kernel_mod_test integer(kind=i_def), parameter :: arr_len = undf -! integer(kind=i_def), parameter :: mode_dimen = n_radaer_mode + 1_i_def -! integer(kind=i_def), parameter :: mmr_arr_len = (nlayers + 1) * (n_radaer_mode + 1) - - contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -53,6 +48,14 @@ contains class(aviation_icao_heights_test_type), intent(inout) :: this +! call feign_extrusion_config( planet_radius=6000000.0_r_def, & +! domain_height=10.0_r_def, & +! method=method_uniform, & +! number_of_layers=2_i_def, & +! stretching_method=stretching_method_linear, & +! stretching_height=1.0_r_def, & +! eta_values=(/0.5_r_def/) ) + allocate(this%icao_heights(arr_len), & this%pressure(arr_len), & this%expect(arr_len) ) @@ -74,7 +77,8 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @test subroutine test_rmdi ( this ) - use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_code + use aviation_icao_heights_kernel_mod, only: aviation_icao_heights_kernel_code + implicit none class(aviation_icao_heights_test_type), intent(inout) :: this From c9eeb6cdbc15c5dae6eddbe1fc0c9861f86c7b67 Mon Sep 17 00:00:00 2001 From: bblay Date: Thu, 11 Jun 2026 18:52:31 +0100 Subject: [PATCH 23/26] need to ask for help --- .../unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf index c78c28302..27778b088 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf @@ -77,7 +77,7 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @test subroutine test_rmdi ( this ) - use aviation_icao_heights_kernel_mod, only: aviation_icao_heights_kernel_code +! use aviation_icao_heights_kernel_mod, only: aviation_icao_heights_kernel_code implicit none From f9e332a6197afbd85054b673bd2a5243fdd9291e Mon Sep 17 00:00:00 2001 From: bblay Date: Fri, 12 Jun 2026 15:06:38 +0100 Subject: [PATCH 24/26] finally got the unit tests compiling --- .../algorithm/aviation_diags_alg_mod.x90 | 9 ++- .../aviation_icao_heights_kernel_mod.F90 | 12 ++-- .../aviation_icao_heights_kernel_mod_test.pf | 65 +++++++++---------- 3 files changed, 44 insertions(+), 42 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 index 636b51936..2c60cca7e 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -202,6 +202,10 @@ contains subroutine icao_heights(convection_fields) ! Calculate ICAO heights of top and base of cloud from their pressures. + ! we can't access this from the kernel + use planet_constants_mod, only: g_over_r + + ! Arguments. type(field_collection_type), intent(inout) :: convection_fields @@ -220,6 +224,7 @@ contains 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') @@ -243,7 +248,7 @@ contains call invoke( & setval_c(conv_cloud_base_icao_height, rmdi), & aviation_icao_heights_kernel_type(conv_cloud_base_icao_height, & - conv_cloud_base_pressure)) + conv_cloud_base_pressure, real(g_over_r, kind=r_def))) call conv_cloud_base_icao_height%write_field() @@ -260,7 +265,7 @@ contains call invoke( & setval_c(conv_cloud_top_icao_height, rmdi), & aviation_icao_heights_kernel_type(conv_cloud_top_icao_height, & - conv_cloud_top_pressure)) + conv_cloud_top_pressure, real(g_over_r, kind=r_def))) call conv_cloud_top_icao_height%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 index 89c205a61..c05d3d120 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -25,13 +25,16 @@ module aviation_icao_heights_kernel_mod ! The aviation diagnostics kernel type. type, extends(kernel_type) :: aviation_icao_heights_kernel_type - type(arg_type), dimension(2) :: meta_args = (/ & + type(arg_type), dimension(3) :: meta_args = (/ & ! Output icao height. arg_type(gh_field, gh_real, gh_write, any_discontinuous_space_1), & ! Input pressure - arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_1) & + arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_1), & + + ! g_over_r + arg_type(gh_scalar, gh_real, gh_read) & /) @@ -48,15 +51,15 @@ module aviation_icao_heights_kernel_mod !> @param[in] nlayers The number of layers in a column. !> @param[out] icao_height Output icao height in kft. !> @param[in] pressure_field Pressure in pa. + !> @param[in] g_over_r Constant must be passed in. !> @param[in] ndf Number of DOFs in the cell. !> @param[in] undf Number of DOFs in the field. !> @param[in] map Dofmap to the bottom cell. subroutine aviation_icao_heights_kernel_code(nlayers, & - icao_height, pressure_field, & + icao_height, pressure_field, g_over_r, & ndf, undf, map) use constants_mod, only: rmdi - use planet_constants_mod, only: g_over_r use science_aviation_constants_mod, only: mtokft, & isa_lapse_ratel, isa_lapse_rateu, isa_press_bot, & isa_press_mid, isa_press_top, isa_temp_bot, isa_temp_top, & @@ -73,6 +76,7 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & ! Arguments (algorithm) real(kind=r_def), intent(out), dimension(undf) :: icao_height real(kind=r_def), intent(in), dimension(undf) :: pressure_field + real(kind=r_def), intent(in) :: g_over_r ! Local variables diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf index 27778b088..df8514cbd 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf @@ -10,12 +10,12 @@ module aviation_icao_heights_kernel_mod_test use funit use constants_mod, only: i_def, r_def, rmdi -! use feign_config_mod, only: feign_extrusion_config + use aviation_icao_heights_kernel_mod, only: aviation_icao_heights_kernel_code implicit none private - public :: test_rmdi + public :: test_rmdi, test_low_safeguard @TestCase type, extends(TestCase), public :: aviation_icao_heights_test_type @@ -24,10 +24,14 @@ module aviation_icao_heights_kernel_mod_test real(kind=r_def), allocatable :: icao_heights(:) real(kind=r_def), allocatable :: expect(:) + ! g / dry air gas constant. Can't seem to import from planet_constants. + real(kind=r_def) :: g_over_r = 0.034163195_r_def + contains procedure setUp procedure tearDown procedure test_rmdi + procedure test_low_safeguard end type aviation_icao_heights_test_type @@ -48,14 +52,6 @@ contains class(aviation_icao_heights_test_type), intent(inout) :: this -! call feign_extrusion_config( planet_radius=6000000.0_r_def, & -! domain_height=10.0_r_def, & -! method=method_uniform, & -! number_of_layers=2_i_def, & -! stretching_method=stretching_method_linear, & -! stretching_height=1.0_r_def, & -! eta_values=(/0.5_r_def/) ) - allocate(this%icao_heights(arr_len), & this%pressure(arr_len), & this%expect(arr_len) ) @@ -77,42 +73,39 @@ contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @test subroutine test_rmdi ( this ) -! use aviation_icao_heights_kernel_mod, only: aviation_icao_heights_kernel_code - implicit none - class(aviation_icao_heights_test_type), intent(inout) :: this - this%pressure = (/ rmdi, -999.0_r_def /) - this%expect = (/ rmdi, -999.0_r_def /) + this%pressure = (/ rmdi, -999.0_r_def /) + this%expect = (/ rmdi, -999.0_r_def /) this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) -! call aviation_icao_heights_code( nlayers, & -! this%icao_heights, this%pressure, & -! ndf, undf, map ) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) @assertEqual( this%icao_heights, this%expect ) end subroutine test_rmdi !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! @test -! subroutine test_low_safeguard ( this ) -! use aviation_icao_heights_kernel_mod, only : aviation_icao_heights_code -! implicit none -! -! class(aviation_icao_heights_test_type), intent(inout) :: this -! -! this%pressure = (/ 500.0_r_def, -999.0_r_def /) -! this%expect = (/ 1234.0_r_def, -999.0_r_def /) -! -! this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) -! call aviation_icao_heights_code( nlayers, & -! this%icao_heights, this%pressure, & -! ndf, undf, map ) -! -! @assertEqual( this%icao_heights, this%expect ) -! -! end subroutine test_low_safeguard + @test + subroutine test_low_safeguard ( this ) + implicit none + class(aviation_icao_heights_test_type), intent(inout) :: this + + real(kind=r_def), parameter :: tolerance = 1.0e-3_r_def + + this%pressure = (/ 500.0_r_def, -999.0_r_def /) + this%expect = (/ 101.886_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) + + @assertEqual( this%expect, this%icao_heights, tolerance ) + + end subroutine test_low_safeguard end module aviation_icao_heights_kernel_mod_test From f9f8838fa49414cc0462ab89aa273249ffe2631f Mon Sep 17 00:00:00 2001 From: bblay Date: Mon, 15 Jun 2026 10:52:28 +0100 Subject: [PATCH 25/26] yes, good unit tests! --- .../aviation_icao_heights_kernel_mod_test.pf | 114 ++++++++++++++++-- 1 file changed, 105 insertions(+), 9 deletions(-) diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf index df8514cbd..5b5419fa4 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/aviation_icao_heights_kernel_mod_test.pf @@ -11,11 +11,15 @@ module aviation_icao_heights_kernel_mod_test use funit use constants_mod, only: i_def, r_def, rmdi use aviation_icao_heights_kernel_mod, only: aviation_icao_heights_kernel_code + use science_aviation_constants_mod, only: isa_press_bot, isa_press_mid, isa_press_top implicit none private - public :: test_rmdi, test_low_safeguard + public :: test_rmdi, & + test_low_low, test_low_high, & + test_mid_low, test_mid_high, & + test_high_low, test_high_limit @TestCase type, extends(TestCase), public :: aviation_icao_heights_test_type @@ -26,12 +30,16 @@ module aviation_icao_heights_kernel_mod_test ! g / dry air gas constant. Can't seem to import from planet_constants. real(kind=r_def) :: g_over_r = 0.034163195_r_def + real(kind=r_def) :: tolerance = 1.0e-3_r_def + contains procedure setUp procedure tearDown procedure test_rmdi - procedure test_low_safeguard + procedure test_low_low, test_low_high + procedure test_mid_low, test_mid_high + procedure test_high_low, test_high_limit end type aviation_icao_heights_test_type @@ -70,7 +78,7 @@ contains end subroutine tearDown - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! mdi handling @test subroutine test_rmdi ( this ) implicit none @@ -88,15 +96,103 @@ contains end subroutine test_rmdi - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! pressure in the lower section - near surface + @test + subroutine test_low_low ( this ) + implicit none + class(aviation_icao_heights_test_type), intent(inout) :: this + + this%pressure = (/ isa_press_bot - 1.0_r_def, -999.0_r_def /) + this%expect = (/ 0.0_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) + + @assertEqual( this%expect, this%icao_heights, this%tolerance ) + + end subroutine test_low_low + + ! pressure in the lower section - near isa_press_mid + @test + subroutine test_low_high ( this ) + implicit none + class(aviation_icao_heights_test_type), intent(inout) :: this + + this%pressure = (/ isa_press_mid + 1.0_r_def, -999.0_r_def /) + this%expect = (/ 36.088_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) + + @assertEqual( this%expect, this%icao_heights, this%tolerance ) + + end subroutine test_low_high + + ! pressure in the mid section - near isa_press_mid + @test + subroutine test_mid_low ( this ) + implicit none + class(aviation_icao_heights_test_type), intent(inout) :: this + + this%pressure = (/ isa_press_mid - 1.0_r_def, -999.0_r_def /) + this%expect = (/ 36.090_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) + + @assertEqual( this%expect, this%icao_heights, this%tolerance ) + + end subroutine test_mid_low + + ! pressure in the mid section - near isa_press_top @test - subroutine test_low_safeguard ( this ) + subroutine test_mid_high ( this ) implicit none class(aviation_icao_heights_test_type), intent(inout) :: this - real(kind=r_def), parameter :: tolerance = 1.0e-3_r_def + this%pressure = (/ isa_press_top + 1.0_r_def, -999.0_r_def /) + this%expect = (/ 65.613_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) + + @assertEqual( this%expect, this%icao_heights, this%tolerance ) + + end subroutine test_mid_high + + ! pressure in the top section - near isa_press_top + @test + subroutine test_high_low ( this ) + implicit none + class(aviation_icao_heights_test_type), intent(inout) :: this + + this%pressure = (/ isa_press_top - 1.0_r_def, -999.0_r_def /) + this%expect = (/ 65.620_r_def, -999.0_r_def /) + + this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) + call aviation_icao_heights_kernel_code( nlayers, & + this%icao_heights, this%pressure, this%g_over_r, & + ndf, undf, map ) + + @assertEqual( this%expect, this%icao_heights, this%tolerance ) + + end subroutine test_high_low + + ! very low pressure changed to 1000 + @test + subroutine test_high_limit ( this ) + implicit none + class(aviation_icao_heights_test_type), intent(inout) :: this - this%pressure = (/ 500.0_r_def, -999.0_r_def /) + this%pressure = (/ 500.0_r_def, -999.0_r_def /) this%expect = (/ 101.886_r_def, -999.0_r_def /) this%icao_heights = (/ -999.0_r_def, -999.0_r_def /) @@ -104,8 +200,8 @@ contains this%icao_heights, this%pressure, this%g_over_r, & ndf, undf, map ) - @assertEqual( this%expect, this%icao_heights, tolerance ) + @assertEqual( this%expect, this%icao_heights, this%tolerance ) - end subroutine test_low_safeguard + end subroutine test_high_limit end module aviation_icao_heights_kernel_mod_test From be51855de092a3bd2e849387478638102045bb64 Mon Sep 17 00:00:00 2001 From: Byron Blay Date: Mon, 15 Jun 2026 11:09:48 +0100 Subject: [PATCH 26/26] Apply suggestions from code review Co-authored-by: iboutle <135141261+iboutle@users.noreply.github.com> --- .../source/kernel/aviation_icao_heights_kernel_mod.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 index c05d3d120..6c9fae93a 100644 --- a/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_icao_heights_kernel_mod.F90 @@ -14,11 +14,10 @@ module aviation_icao_heights_kernel_mod use argument_mod, only: arg_type, & gh_field, gh_scalar, gh_logical, & - gh_read, gh_write, gh_integer, & + gh_read, gh_write, & gh_real, cell_column, & - any_discontinuous_space_1, & - any_discontinuous_space_2 - use constants_mod, only: r_def, i_def, l_def + any_discontinuous_space_1 + use constants_mod, only: r_def, i_def use kernel_mod, only: kernel_type implicit none @@ -80,7 +79,6 @@ subroutine aviation_icao_heights_kernel_code(nlayers, & ! Local variables - integer(kind=i_def) :: df real(kind=r_def) :: zp1 real(kind=r_def) :: zp2