Stable diags #11, #21, #24 #483
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,279 @@ | ||
| !----------------------------------------------------------------------------- | ||
|
Contributor
|
||
| ! (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, 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_level_always, log_level_error, & | ||
| log_level_info, log_level_debug | ||
|
|
||
| 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 | ||
|
|
||
| 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) :: thick_850_flag, thick_500_flag, snow_probability_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 output fields we produce. | ||
| type( field_type ) :: thick_850, thick_500, snow_probability | ||
|
|
||
| integer(i_def) :: i | ||
|
|
||
| character(len=str_long) :: message | ||
|
|
||
|
|
||
| ! Check the request flags. | ||
| 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. (thick_850_flag .or. & | ||
| thick_500_flag .or. & | ||
| snow_probability_flag) ) then | ||
| return | ||
| end if | ||
|
|
||
| if ( .not. use_xios_io ) then | ||
|
Contributor
There was a problem hiding this comment. I would suggest removing this - xios is really the only option, so it's doing lots of if-tests for no real benefit! |
||
| return | ||
| endif | ||
|
|
||
| ! Initialise required fields to 1. | ||
| if ( snow_probability_flag ) then | ||
| message = 'Section 20: snow probability requested' | ||
|
Contributor
There was a problem hiding this comment. I guess these are useful for debugging, but are they needed on main? I believe init_diag prints this information under some log level (perhaps trace instead of debug) |
||
| call log_event( message, log_level_debug ) | ||
| call invoke( setval_c(snow_probability, 0.0_r_def)) | ||
|
Contributor
There was a problem hiding this comment. I don't think these pre-initialisations to 0 or 1 are needed, because the field is entirely written to inside the kernel, so can be removed |
||
| end if | ||
|
|
||
| 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(thick_850, 1.0_r_def)) | ||
| end if | ||
|
|
||
| if ( thick_500_flag ) then | ||
| message = 'Section 20: geopot thick 500 requested' | ||
| call log_event( message, log_level_debug ) | ||
| call invoke( setval_c(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_error ) | ||
| return | ||
|
Contributor
There was a problem hiding this comment. log_level_error should terminate the model, so the returns are redundant |
||
| 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_error ) | ||
| 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_error ) | ||
| return | ||
| end if | ||
|
|
||
| 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 (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 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( & | ||
| thick_850, thick_500, snow_probability, & | ||
| plev_geopot, & | ||
| thick_850_flag, thick_500_flag, snow_probability_flag, & | ||
| i1000, i850, i500)) | ||
|
|
||
| ! Write the fields | ||
| 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 | ||
| deallocate(plevs) | ||
| end if | ||
|
|
||
| end subroutine from_plev_geopot | ||
|
|
||
|
|
||
| subroutine icao_heights(convection_fields) | ||
|
Contributor
There was a problem hiding this comment. Sorry, I don't think this is the right place to be calling this code - I think it should be called directly from the convection diagnostics routines, the same as happens in the UM. Currently I don't believe this will be working correctly, specifically if the ICAO height has been requested but the pressure diagnostic hasn't. I would suggest removing this. In conv_diags_mod (and the Comorph equivalent) you will need to add the init_diag calls just before the pres_cv_base/top is initialised, and use the output of this to force the pres_cv_base/top to be initialised if the ICAO equivalent has been requested. Then in the output section where the pres_cv_base/top are written, you can call the icao_kernel on these fields and then output there (if the field is requested). You shouldn't need the pre-initialisation to rmdi after making my other suggested changes to the kernel. |
||
| ! 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 | ||
|
|
||
| ! Locals. | ||
| character(len=str_long) :: message | ||
|
|
||
| ! 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 | ||
| message = 'section 20: cloud icao heights on' | ||
| call log_event( message, 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), & | ||
| aviation_icao_heights_kernel_type(conv_cloud_base_icao_height, & | ||
| conv_cloud_base_pressure, real(g_over_r, kind=r_def))) | ||
|
|
||
| call conv_cloud_base_icao_height%write_field() | ||
|
|
||
| ! Is it necessary/safe to remove this field now? | ||
|
Author
There was a problem hiding this comment. Need steering on this question please. |
||
| ! Did the get trigger an unpack which we need to now remove? | ||
| ! call convection_fields%remove_field('convection__pres_cv_base') | ||
|
|
||
| end if | ||
|
|
||
| ! 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), & | ||
| aviation_icao_heights_kernel_type(conv_cloud_top_icao_height, & | ||
| conv_cloud_top_pressure, real(g_over_r, kind=r_def))) | ||
|
|
||
| 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 | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Couldn't find documentation for the
<file>and<field_group>tags.