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..e472a05d4 100644 --- a/applications/lfric_atm/example/iodef.xml +++ b/applications/lfric_atm/example/iodef.xml @@ -94,6 +94,13 @@ + + + + + + + diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 007ee86df..94d045e8b 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -989,6 +989,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..d56908851 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/algorithm/aviation_diags_alg_mod.x90 @@ -0,0 +1,253 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2025 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- +! +! The main entry point for calculating section 20 aviation diagnostics. +! +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file currently belongs in section: physics_schemes_interface +! whilst discussions are ongoing about its final location. +! +MODULE aviation_diags_alg_mod + + USE constants_mod, ONLY: r_def, i_def, l_def, RMDI + USE field_mod, ONLY: field_type + USE field_collection_mod, ONLY: field_collection_type + USE io_config_mod, ONLY: USE_xios_io + USE initialise_diagnostics_mod, ONLY: init_diag => init_diagnostic_field + + USE lfric_xios_diag_mod, ONLY: get_axis_dimension, get_axis_values + USE log_mod, ONLY: log_event, log_scratch_space, & + LOG_LEVEL_ALWAYS, LOG_LEVEL_ERROR, & + LOG_LEVEL_INFO, LOG_LEVEL_DEBUG + + USE diags_geopot_kernel_mod, ONLY: diags_geopot_kernel_type + USE diags_icao_heights_kernel_mod, ONLY: diags_icao_heights_kernel_type + + IMPLICIT NONE + +PRIVATE +PUBLIC :: aviation_diags_alg + +CONTAINS + + ! Algorithm to calculate the section 20 aviation diagnostics. + SUBROUTINE aviation_diags_alg(plev_geopot, convection_fields) + + IMPLICIT NONE + + ! Arguments + TYPE(field_type), INTENT(IN) :: plev_geopot + TYPE(field_collection_type), INTENT(INOUT) :: convection_fields + + + ! Calculate thickness from geopotential height at pressure levels. + CALL from_plev_geopot(plev_geopot) + + + ! Calculate icao heights of top and base of cloud from their pressures. + CALL icao_heights(convection_fields) + + + END SUBROUTINE aviation_diags_alg + + SUBROUTINE from_plev_geopot(plev_geopot) + ! Calculate diagnostics from plvel_geopot. + ! - Thickness from geopotential height at pressure levels. + + ! Arguments + TYPE(field_type), INTENT(IN) :: plev_geopot + + ! Local variables + + ! These flags tell us which results are requested at this time step. + LOGICAL(l_def) :: aviation_thick_850_flag, aviation_thick_500_flag + + ! The array of pressure levels. + INTEGER(I_DEF) :: nplev + REAL(R_DEF), ALLOCATABLE :: plevs(:) + INTEGER(I_DEF) :: plevs_alloc_stat + + ! Level indices for the three pressure levels we want to read. + INTEGER(I_DEF) :: i1000, i850, i500 + + ! Approximate equality tolerance - is there a common variable somewhere? + REAL(R_DEF), PARAMETER :: plev_tol = 0.1_r_def + + ! The two output fields we produce. + TYPE( field_type ) :: aviation_thick_850, aviation_thick_500 + + INTEGER(I_DEF) :: i + + + ! Check the request flags. + aviation_thick_850_flag = & + init_diag(aviation_thick_850, 'aviation__geopot_thickness_850') + aviation_thick_500_flag = & + init_diag(aviation_thick_500, 'aviation__geopot_thickness_500') + IF ( aviation_thick_850_flag .OR. aviation_thick_500_flag ) THEN + WRITE(log_scratch_space, '(A)') 'Section 20 thickness is on' + CALL log_event(log_scratch_space, LOG_LEVEL_DEBUG) + ELSE + ! Nothing requested at this time step. + RETURN + END IF + + ! Get the array of pressure levels. + nplev = get_axis_dimension('pressure_levels') + IF (nplev <= 0) THEN + WRITE(log_scratch_space, '(A, I0)') 'No pressure levels, nplev=', nplev + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + ALLOCATE(plevs(nplev), stat=plevs_alloc_stat) + IF (plevs_alloc_stat /= 0) THEN + WRITE(log_scratch_space, '(A, I0)') 'allocate(plevs) failed, stat=', & + plevs_alloc_stat + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + plevs = get_axis_values('pressure_levels',nplev) + + ! Find the level indices for our three pressures of interest. + ! Leave as -1 if not requested, telling the kernel not to calculate it. + ! Assumes hPa. Should we worry about other units? + i1000 = -1 + i850 = -1 + i500 = -1 + DO i = 1, nplev + + ! 1000? + IF ( abs(plevs(i) - 100000.0_r_def) < plev_tol ) THEN + i1000 = i + + ! 850? + ELSE IF ( abs(plevs(i) - 85000.0_r_def) < plev_tol ) THEN + i850 = i + + ! 500? + ELSE IF ( abs(plevs(i) - 50000.0_r_def) < plev_tol ) THEN + i500 = i + + END IF + END DO + + ! Check we found the required levels. + IF (i1000 == -1) THEN + WRITE(log_scratch_space, '(A)') 'could not find 1000hPa' + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + IF (aviation_thick_850_flag .AND. i850 == -1) THEN + WRITE(log_scratch_space, '(A)') 'could not find 850hPa' + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + IF (aviation_thick_500_flag .AND. i500 == -1) THEN + WRITE(log_scratch_space, '(A)') 'could not find 500hPa' + CALL log_event(log_scratch_space, LOG_LEVEL_ERROR) + RETURN + END IF + + ! Call the kernel to subtract the levels. + ! Todo: Comment why we initialise to 1.0 and not 0.0. I can't remember. + CALL invoke( & + setval_c(aviation_thick_850, 1.0_r_def), & + setval_c(aviation_thick_500, 1.0_r_def), & + diags_geopot_kernel_type( & + aviation_thick_850, aviation_thick_500, & + plev_geopot, & + aviation_thick_850_flag, aviation_thick_500_flag, & + i1000, i850, i500) & + ) + + ! Write the fields + IF (aviation_thick_850_flag) CALL aviation_thick_850%write_field() + IF (aviation_thick_500_flag) CALL aviation_thick_500%write_field() + + ! Clean up. + IF (allocated(plevs)) THEN + DEALLOCATE(plevs) + END IF + + END SUBROUTINE from_plev_geopot + + + SUBROUTINE icao_heights(convection_fields) + ! Calculate iaco heights of top and base of cloud from their pressures. + + ! Arguments + TYPE(field_collection_type), INTENT(INOUT) :: convection_fields + + ! Locals + + ! Pressure fields to unpack. + type(field_type), pointer :: conv_cloud_base_pressure + type(field_type), pointer :: conv_cloud_top_pressure + + ! Height fields to create and write. + TYPE(field_type) :: conv_cloud_base_icao_height + TYPE(field_type) :: conv_cloud_top_icao_height + + ! Request flags. + LOGICAL(l_def) :: base_height_flag + LOGICAL(l_def) :: top_height_flag + + ! which heights are requested this time step? + top_height_flag = init_diag( & + conv_cloud_top_icao_height, 'aviation__conv_cloud_top_icao_height') + + base_height_flag = init_diag( & + conv_cloud_base_icao_height, 'aviation__conv_cloud_base_icao_height') + + if (base_height_flag .or. top_height_flag) then + WRITE(log_scratch_space, '(A)') 'Section 20 cloud ICAO heights on' + CALL log_event(log_scratch_space, LOG_LEVEL_DEBUG) + else + ! nothing requested at this time step + return + end if + + + ! cloud base + IF (base_height_flag) THEN + call convection_fields%get_field( & + 'convection__pres_cv_base', conv_cloud_base_pressure) + call invoke( & + setval_c(conv_cloud_base_icao_height, RMDI), & + diags_icao_heights_kernel_type(conv_cloud_base_icao_height, & + conv_cloud_base_pressure)) + + call conv_cloud_base_icao_height%write_field() + + ! Is it necessary/safe to remove this field now? + ! Did the get trigger an unpack which we need to now remove? + ! Based on https://code.metoffice.gov.uk/trac/lfric_apps/ticket/189 + ! call convection_fields%remove_field('convection__pres_cv_base') + + END IF + + ! could top + IF (top_height_flag) THEN + call convection_fields%get_field( & + 'convection__pres_cv_top', conv_cloud_top_pressure) + call invoke( & + setval_c(conv_cloud_top_icao_height, RMDI), & + diags_icao_heights_kernel_type(conv_cloud_top_icao_height, & + conv_cloud_top_pressure)) + + call conv_cloud_top_icao_height%write_field() + + ! Is it necessary/safe to remove this field now? + ! call convection_fields%remove_field('convection__pres_cv_top') + + END IF + + END SUBROUTINE icao_heights + +END MODULE aviation_diags_alg_mod 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 489742c65..df47e248d 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 @@ -27,6 +27,8 @@ module conv_gr_alg_mod use log_mod, only: log_event, LOG_LEVEL_DEBUG use mesh_mod, only: mesh_type + use initialise_diagnostics_mod, only : diag_samp => diagnostic_to_be_sampled + implicit none private @@ -62,14 +64,14 @@ contains type( field_type ), intent( inout ) :: mr(nummr) type( field_type ), intent( in ) :: rho, exner - type( field_collection_type ), intent(in) :: derived_fields - type( field_collection_type ), intent(in) :: turbulence_fields - type( field_collection_type ), intent(in) :: convection_fields - type( field_collection_type ), intent(in) :: cloud_fields - type( field_collection_type ), intent(in) :: surface_fields - type( field_collection_type ), intent(in) :: chemistry_fields - type( field_collection_type ), intent(in) :: aerosol_fields - type( field_collection_type ), intent(in) :: microphysics_fields + type( field_collection_type ), intent(in) :: derived_fields + type( field_collection_type ), intent(in) :: turbulence_fields + type( field_collection_type ), intent(inout) :: convection_fields + type( field_collection_type ), intent(in) :: cloud_fields + type( field_collection_type ), intent(in) :: surface_fields + type( field_collection_type ), intent(in) :: chemistry_fields + type( field_collection_type ), intent(in) :: aerosol_fields + type( field_collection_type ), intent(in) :: microphysics_fields integer(kind=i_def), intent(in) :: outer @@ -767,6 +769,20 @@ contains dth_conv_noshal, & dmv_conv_noshal) + + ! Aviation dependencies to activate pres_cv_base and pres_cv_top + if (diag_samp('aviation__conv_cloud_base_icao_height')) then + if (.not. convection_fields%field_exists('convection__pres_cv_base')) then + call convection_fields%add_field(pres_cv_base) + end if + end if + + if (diag_samp('aviation__conv_cloud_top_icao_height')) then + if (.not. convection_fields%field_exists('convection__pres_cv_top')) then + call convection_fields%add_field(pres_cv_top) + end if + end if + end if nullify( mesh ) 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/diagnostics/conv_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/conv_diags_mod.x90 index 1d06d96e8..f3e6455b1 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/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/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 661368f12..287e7f3c2 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,16 @@ - + + + + + + + + + + diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 index aff34838f..ebe8d8af8 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, 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 +