From 39db4989ea984cf3270635cf5e764f586ecf7239 Mon Sep 17 00:00:00 2001 From: bblay Date: Fri, 12 Dec 2025 16:33:20 +0000 Subject: [PATCH] s20 diag thickness --- CONTRIBUTORS.md | 3 +- applications/lfric_atm/example/iodef.xml | 5 + .../lfric_atm/metadata/field_def_diags.xml | 3 + .../algorithm/aviation_diags_alg_mod.x90 | 162 ++++++++++++++++++ .../algorithm/pres_lev_diags_alg_mod.x90 | 23 ++- .../kernel/aviation_diags_kernel_mod.F90 | 133 ++++++++++++++ .../file/file_def_diags_oper_nwp_gl.xml | 9 +- .../driver/gungho_diagnostics_driver_mod.F90 | 10 +- 8 files changed, 340 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 10e9ced23..8c95cbb81 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -1,6 +1,6 @@ # Contributors | GitHub user | Real Name | Affiliation | Date | -| --------------- | ----------------- | ----------- | ---------- | +|-----------------|-------------------| ----------- |------------| | james-bruten-mo | James Bruten | Met Office | 2025-12-09 | | jedbakerMO | Jed Baker | Met Office | 2025-12-29 | | jennyhickson | Jenny Hickson | Met Office | 2025-12-10 | @@ -8,6 +8,7 @@ | mo-marqh | mark Hedley | Met Office | 2025-12-11 | | yaswant | Yaswant Pradhan | Met Office | 2025-12-16 | | oakleybrunt | Oakley Brunt | Met Office | 2025-12-19 | +| bblay-mo | Byron Blay | Met Office | 2026-01-07 | | harry-shepherd | Harry Shepherd | Met Office | 2026-01-08 | | DrTVockerodtMO | Terence Vockerodt | Met Office | 2026-01-08 | | MetBenjaminWent | Benjamin Went | Met Office | 2026-01-15 | 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 007ee86df..ff1cb913b 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -989,6 +989,9 @@ + + + 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..2783df502 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -0,0 +1,162 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2025 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- +! +! The main entry point for calculating section 20 aviation diagnostics. +! +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file currently belongs in section: physics_schemes_interface +! whilst discussions are ongoing about its final location. +! +MODULE aviation_diags_alg_mod + + USE constants_mod, ONLY: r_def, i_def, l_def + USE field_mod, ONLY: field_type + USE io_config_mod, ONLY: USE_xios_io + USE initialise_diagnostics_mod, ONLY: init_diag => init_diagnostic_field + + USE lfric_xios_diag_mod, ONLY: get_axis_dimension, get_axis_values + USE log_mod, ONLY: log_event, log_scratch_space, & + LOG_LEVEL_ALWAYS, LOG_LEVEL_ERROR, & + LOG_LEVEL_INFO, LOG_LEVEL_DEBUG + + USE aviation_diags_kernel_mod, ONLY: aviation_diags_kernel_type + + IMPLICIT NONE + +PRIVATE +PUBLIC :: aviation_diags_alg + +CONTAINS + + ! Algorithm to calculate the section 20 aviation diagnostics. + SUBROUTINE aviation_diags_alg(plev_geopot) + + IMPLICIT NONE + + ! Arguments + + ! We calculate the result from this input field. + TYPE(field_type), INTENT(IN) :: plev_geopot + + + ! Local variables + + ! These flags tell us which results are requested at this time step. + LOGICAL(l_def) :: aviation_thick_850_flag, aviation_thick_500_flag + + ! The array of pressure levels. + INTEGER(I_DEF) :: nplev + REAL(R_DEF), ALLOCATABLE :: plevs(:) + INTEGER(I_DEF) :: plevs_alloc_stat + + ! Level indices for the three pressure levels we want to read. + INTEGER(I_DEF) :: i1000, i850, i500 + + ! Approximate equality tolerance - is there a common variable somewhere? + REAL(R_DEF), PARAMETER :: plev_tol = 0.1_r_def + + ! The two output fields we produce. + TYPE( field_type ) :: aviation_thick_850, aviation_thick_500 + + INTEGER(I_DEF) :: i + + + ! Check the request flags. + aviation_thick_850_flag = & + init_diag(aviation_thick_850, 'aviation__geopot_thickness_850') + aviation_thick_500_flag = & + init_diag(aviation_thick_500, 'aviation__geopot_thickness_500') + IF ( aviation_thick_850_flag .OR. aviation_thick_500_flag ) THEN + WRITE(log_scratch_space, '(A)') 'Section 20 thickness is on' + CALL log_event(log_scratch_space, LOG_LEVEL_DEBUG) + ELSE + ! Nothing requested at this time step. + RETURN + END IF + + ! Get the array of pressure levels. + nplev = get_axis_dimension('pressure_levels') + IF (nplev <= 0) THEN + WRITE(log_scratch_space, '(A, I0)') 'No pressure levels, nplev=', nplev + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + ALLOCATE(plevs(nplev), stat=plevs_alloc_stat) + IF (plevs_alloc_stat /= 0) THEN + WRITE(log_scratch_space, '(A, I0)') 'allocate(plevs) failed, stat=', & + plevs_alloc_stat + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + plevs = get_axis_values('pressure_levels',nplev) + + ! Find the level indices for our three pressures of interest. + ! Leave as -1 if not requested, telling the kernel not to calculate it. + ! Assumes hPa. Should we worry about other units? + i1000 = -1 + i850 = -1 + i500 = -1 + DO i = 1, nplev + + ! 1000? + IF ( abs(plevs(i) - 100000.0_r_def) < plev_tol ) THEN + i1000 = i + + ! 850? + ELSE IF ( abs(plevs(i) - 85000.0_r_def) < plev_tol ) THEN + i850 = i + + ! 500? + ELSE IF ( abs(plevs(i) - 50000.0_r_def) < plev_tol ) THEN + i500 = i + + END IF + END DO + + ! Check we found the required levels. + IF (i1000 == -1) THEN + WRITE(log_scratch_space, '(A)') 'could not find 1000hPa' + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + IF (aviation_thick_850_flag .AND. i850 == -1) THEN + WRITE(log_scratch_space, '(A)') 'could not find 850hPa' + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + IF (aviation_thick_500_flag .AND. i500 == -1) THEN + WRITE(log_scratch_space, '(A)') 'could not find 500hPa' + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + ! Call the kernel to subtract the levels. + ! Todo: Comment why we initialise to 1.0 and not 0.0. I can't remember. + CALL invoke( & + setval_c(aviation_thick_850, 1.0_r_def), & + setval_c(aviation_thick_500, 1.0_r_def), & + aviation_diags_kernel_type(& + aviation_thick_850, aviation_thick_500, & + plev_geopot, & + aviation_thick_850_flag, aviation_thick_500_flag, & + i1000, i850, i500) & + ) + + ! Write the fields + IF (aviation_thick_850_flag) CALL aviation_thick_850%write_field() + IF (aviation_thick_500_flag) CALL aviation_thick_500%write_field() + + ! Clean up. + IF (allocated(plevs)) THEN + DEALLOCATE(plevs) + END IF + + END SUBROUTINE aviation_diags_alg + +END MODULE aviation_diags_alg_mod 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 883293b37..3d51fa0b1 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, 'pres_lev_diags_alg' ) nplev = get_axis_dimension('pressure_levels') @@ -257,9 +263,16 @@ contains end if ! Geopotential height on pressure levels + aviation_thickness_flag = diag_samp('aviation__geopot_thickness_850') .or. & + diag_samp('aviation__geopot_thickness_500') plev_geopot_clim_flag = diag_samp('plev__geopot_clim') - plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', activate=plev_geopot_clim_flag) - if ((plev_geopot_flag .or. plev_geopot_clim_flag) .and. use_xios_io) then + plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', & + activate=plev_geopot_clim_flag .or. aviation_thickness_flag) + + plev_gepot_required = plev_geopot_flag .or. & + plev_geopot_clim_flag .or. & + aviation_thickness_flag + if (plev_gepot_required .and. use_xios_io) then height_w3 => get_height_fv(W3, exner_w3%get_mesh_id()) call invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, & 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..9e0e3558c --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/aviation_diags_kernel_mod.F90 @@ -0,0 +1,133 @@ +!------------------------------------------------------------------------------- +! (c) Crown copyright 2021 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!------------------------------------------------------------------------------- +! +! Section 20 aviation diagnostics kernel. +! +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file currently belongs in section: physics_schemes_interface +! whilst discussions are ongoing about its final location. +! +MODULE aviation_diags_kernel_mod + + USE argument_mod, ONLY: arg_type, & + GH_FIELD, GH_SCALAR, GH_LOGICAL, & + GH_READ, GH_WRITE, GH_INTEGER, & + GH_REAL, CELL_COLUMN, & + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + USE kernel_mod, ONLY: kernel_type + USE constants_mod, ONLY: r_def, i_def, l_def + + IMPLICIT NONE + + ! The aviation diagnostics kernel type. + TYPE, EXTENDS(kernel_type) :: aviation_diags_kernel_type + TYPE(arg_type), DIMENSION(8) :: meta_args = (/ & + + ! Output fields. + arg_type(gh_field, gh_real, gh_write, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(gh_field, gh_real, gh_write, ANY_DISCONTINUOUS_SPACE_1), & + + ! Source field. + arg_type(gh_field, gh_real, gh_read, ANY_DISCONTINUOUS_SPACE_2), & + + ! Request flags. + arg_type(gh_scalar, gh_logical, gh_read), & + arg_type(gh_scalar, gh_logical, gh_read), & + + ! Level indices. + arg_type(gh_scalar, gh_integer, gh_read), & + arg_type(gh_scalar, gh_integer, gh_read), & + arg_type(gh_scalar, gh_integer, gh_read) & + /) + + INTEGER :: operates_on = cell_column + + CONTAINS + PROCEDURE, NOPASS :: code => aviation_diags_kernel_code + END TYPE aviation_diags_kernel_type + +CONTAINS + + SUBROUTINE aviation_diags_kernel_code(nlayers, & + ! Output fields. + thickness_850, thickness_500, & + + ! Source field. + plev_geopot, & + + ! Request flags. + thickness_850_flag, thickness_500_flag, & + + ! Level incides. + i1000, i850, i500, & + + ! Kernel stuff. + result_ndf, result_undf, result_map, & + source_ndf, source_undf, source_map) + + ! Subtract geopotential heights at 850 and 500 hPa from that at 1000 hPa + ! to calculate thickness fields. + + IMPLICIT NONE + + ! Arguments (kernel) + + ! The number of layers in a column. + INTEGER(KIND=i_def), INTENT(IN) :: nlayers + + ! Number of degrees of freedom (columns) in the cell we're processing. + INTEGER(KIND=i_def), INTENT(IN) :: result_ndf, source_ndf + + ! Number of unique degrees of freedom in the fields. + INTEGER(KIND=i_def), INTENT(IN) :: result_undf, source_undf + + ! Degrees of freedom maps. Offsets to the bottom of each column. + INTEGER(KIND=i_def), INTENT(IN), DIMENSION(result_ndf) :: result_map + INTEGER(KIND=i_def), INTENT(IN), DIMENSION(source_ndf) :: source_map + + + ! Arguments (algorithm) + + ! Output thickness fields. + REAL(KIND=r_def), INTENT(OUT), DIMENSION(result_undf) :: thickness_850 + REAL(KIND=r_def), INTENT(OUT), DIMENSION(result_undf) :: thickness_500 + + ! Geopotential height at pressure levels. + REAL(KIND=r_def), INTENT(IN), DIMENSION(source_undf) :: plev_geopot + + ! Request flags. + LOGICAL(KIND=l_def), INTENT(IN) :: thickness_850_flag, thickness_500_flag + + ! Level indices. For i850 and i500, -1 means "not requested". + INTEGER(KIND=i_def), INTENT(IN) :: i1000, i850, i500 + + + ! Local variables + INTEGER(KIND=i_def) :: df + REAL(KIND=r_def) :: gph_1000 + + + ! Process every DOF in this cell. + DO df = 1, result_ndf + + gph_1000 = plev_geopot(source_map(df) + i1000-1) + + IF (i850 /= -1) THEN + thickness_850(result_map(df)) = & + plev_geopot(source_map(df)+i850-1) - gph_1000 + END IF + + IF (i500 /= -1) THEN + thickness_500(result_map(df)) = & + plev_geopot(source_map(df)+i500-1) - gph_1000 + END IF + + END DO + + END SUBROUTINE aviation_diags_kernel_code + +END MODULE aviation_diags_kernel_mod 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 661368f12..5764cf07a 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 aff34838f..8bd4003a0 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(:) => null() 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 => null() type( field_type), pointer :: u => null() type( field_type), pointer :: h_u => null() @@ -315,9 +321,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')