From 9548db13402f30372844aaa9e0e4b76be7e2c998 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Fri, 24 Jul 2026 15:23:30 +0100 Subject: [PATCH 01/11] 10m t and q --- .../lfric_atm/metadata/field_def_diags.xml | 2 + .../source/algorithm/jules_exp_alg_mod.x90 | 11 +++-- .../source/algorithm/jules_imp_alg_mod.x90 | 13 ++++- .../diagnostics/jules_exp_diags_mod.f90 | 29 +++++++++++- .../diagnostics/jules_imp_diags_mod.x90 | 20 +++++++- .../source/kernel/jules_exp_kernel_mod.F90 | 15 +++++- .../source/kernel/jules_imp_kernel_mod.F90 | 47 ++++++++++++++++++- .../source/psy/psykal_lite_phys_mod.F90 | 10 ++-- 8 files changed, 132 insertions(+), 15 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index a17d884f90..f0399a66eb 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -533,9 +533,11 @@ + + diff --git a/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 b/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 index ee0bb5f8d6..a741650eff 100644 --- a/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 +++ b/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 @@ -90,7 +90,7 @@ contains 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(inout) :: surface_fields type( field_collection_type ), intent(in) :: soil_fields type( field_collection_type ), intent(in) :: snow_fields type( field_collection_type ), intent(in) :: aerosol_fields @@ -219,7 +219,7 @@ contains ! local variables type(mesh_type), pointer :: mesh type( field_type ) :: gross_prim_prod, & - z0h_eff, soil_respiration + z0h_eff, soil_respiration, chr10m type( field_type ) :: mr_ice integer(i_def) :: stencil_depth, ncells, ncells_halo @@ -371,7 +371,7 @@ contains height_w3 => get_height_fv( W3, mesh%get_id() ) call initialise_diags_for_jules_exp(z0h_eff, gross_prim_prod, & - soil_respiration) + soil_respiration, chr10m, mesh) call zh%copy_field_properties(recip_l_mo_sea) call zh%copy_field_properties(rhostar) @@ -447,11 +447,14 @@ contains urbhwr, urbhgt, urbztm, urbdisp, & rhostar, recip_l_mo_sea, & t1_sd_2d, q1_sd_2d, & - gross_prim_prod, z0h_eff, ocn_cpl_point, & + gross_prim_prod, z0h_eff, chr10m, ocn_cpl_point,& stencil_depth) ! Switch UM back to columns call um_sizes_init(1_i_def) + ! Add to surface fields collection for passing to implicit jules + call surface_fields%add_field(chr10m) + if ( LPROF ) call stop_timing( id, 'jules.explicit' ) diff --git a/interfaces/jules_interface/source/algorithm/jules_imp_alg_mod.x90 b/interfaces/jules_interface/source/algorithm/jules_imp_alg_mod.x90 index 355b37f748..77639a16ba 100644 --- a/interfaces/jules_interface/source/algorithm/jules_imp_alg_mod.x90 +++ b/interfaces/jules_interface/source/algorithm/jules_imp_alg_mod.x90 @@ -143,6 +143,7 @@ contains type( field_type ), pointer :: canhc_tile => null() type( field_type ), pointer :: ustar => null() type( field_type ), pointer :: tile_water_extract => null() + type( field_type ), pointer :: chr10m type( integer_field_type ), pointer :: ocn_cpl_point => null() type( field_type ), pointer :: soil_temperature => null() @@ -157,6 +158,7 @@ contains type( field_type ) :: rh1p5m, rh1p5m_ssi, rh1p5m_land type( field_type ) :: t1p5m_surft, q1p5m_surft + type( field_type ) :: t10m_ssi, q10m_ssi type( field_type ) :: snomlt_surf_htf type( field_type ) :: soil_evap type( field_type ) :: soil_surf_ht_flux @@ -221,6 +223,7 @@ contains call surface_fields%get_field('canhc_tile', canhc_tile) call surface_fields%get_field('ustar', ustar) call surface_fields%get_field('tile_water_extract', tile_water_extract) + call surface_fields%get_field('chr10m', chr10m) call surface_fields%get_field('ocn_cpl_point', ocn_cpl_point) ! Soil fields @@ -255,6 +258,7 @@ contains t1p5m, q1p5m, qcl1p5m, rh1p5m, & t1p5m_ssi, q1p5m_ssi, & qcl1p5m_ssi, rh1p5m_ssi, & + t10m_ssi, q10m_ssi, & t1p5m_land, q1p5m_land, & qcl1p5m_land, rh1p5m_land, & latent_heat, & @@ -298,6 +302,7 @@ contains t1p5m_surft, q1p5m_surft, t1p5m, & q1p5m, qcl1p5m, rh1p5m, t1p5m_ssi, & q1p5m_ssi, qcl1p5m_ssi, rh1p5m_ssi, & + chr10m, t10m_ssi, q10m_ssi, & t1p5m_land, q1p5m_land, qcl1p5m_land, & rh1p5m_land, latent_heat, & snomlt_surf_htf, soil_evap, & @@ -305,6 +310,11 @@ contains surf_radnet, surf_lw_up, surf_lw_down,& ocn_cpl_point) ) + ! Remove from field collection as we're done for this timestep + if (loop == 2 .and. outer == outer_iterations) then + call surface_fields%remove_field('chr10m') + end if + if ( LPROF ) call stop_timing( id, 'jules.implicit' ) ! Output the BL diagnostics @@ -317,7 +327,8 @@ contains t1p5m_surft, q1p5m_surft, & t1p5m, q1p5m, qcl1p5m, rh1p5m, & t1p5m_ssi, q1p5m_ssi, qcl1p5m_ssi, & - rh1p5m_ssi, t1p5m_land, q1p5m_land, & + rh1p5m_ssi, t10m_ssi, q10m_ssi, & + t1p5m_land, q1p5m_land, & qcl1p5m_land, rh1p5m_land, latent_heat, & snomlt_surf_htf, & soil_evap, & diff --git a/interfaces/jules_interface/source/diagnostics/jules_exp_diags_mod.f90 b/interfaces/jules_interface/source/diagnostics/jules_exp_diags_mod.f90 index 4bac420e35..83a1e17d40 100644 --- a/interfaces/jules_interface/source/diagnostics/jules_exp_diags_mod.f90 +++ b/interfaces/jules_interface/source/diagnostics/jules_exp_diags_mod.f90 @@ -12,7 +12,15 @@ module jules_exp_diags_mod use integer_field_mod, only: integer_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, & + samp_diag => diagnostic_to_be_sampled + use empty_data_mod, only : empty_real_data + use extrusion_mod, only: TWOD + use mesh_collection_mod, only: mesh_collection + use mesh_mod, only: mesh_type + use function_space_collection_mod, only: function_space_collection + use function_space_mod, only: function_space_type + use fs_continuity_mod, only: W3 implicit none @@ -31,14 +39,21 @@ module jules_exp_diags_mod !> @param[inout] z0h_eff Gridbox mean effective roughness length for scalars !> @param[inout] gross_prim_prod Gross Primary Productivity !> @param[inout] soil_respiration Soil heterotrophic respiration + !> @param[inout] chr10m 10m transfer coefficient subroutine initialise_diags_for_jules_exp(z0h_eff, gross_prim_prod, & - soil_respiration) + soil_respiration, chr10m, mesh) implicit none type( field_type ), intent(inout) :: z0h_eff type( field_type ), intent(inout) :: gross_prim_prod type( field_type ), intent(inout) :: soil_respiration + type( field_type ), intent(inout) :: chr10m + + type( mesh_type ), intent(in), pointer :: mesh + type(mesh_type), pointer :: twod_mesh + type(function_space_type), pointer :: vector_space + integer( tik ) :: id if ( LPROF ) call start_timing( id, 'diags.jules_exp' ) @@ -46,6 +61,16 @@ subroutine initialise_diags_for_jules_exp(z0h_eff, gross_prim_prod, & z0h_eff_flag = init_diag(z0h_eff, 'surface__z0h_eff') gross_prim_prod_flag = init_diag(gross_prim_prod, 'surface__gross_prim_prod') soil_respiration_flag = init_diag(soil_respiration, 'surface__soil_respiration') + ! chr10m is required if either t10m_ssi or q10m_ssi are requested + ! but needed as an empty field otherwise + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + vector_space => function_space_collection%get_fs(twod_mesh, 0, 0, W3) + if (samp_diag('surface__t10m_ssi') .or. samp_diag('surface__q10m_ssi')) then + call chr10m%initialise(vector_space, name='chr10m') + else + call chr10m%initialise(vector_space, name='chr10m', & + override_data = empty_real_data) + end if if ( LPROF ) call stop_timing( id, 'diags.jules_exp' ) diff --git a/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 b/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 index 13a36d4bf6..8a98a2c164 100644 --- a/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 +++ b/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 @@ -32,6 +32,8 @@ module jules_imp_diags_mod logical( l_def ) :: q1p5m_ssi_flag logical( l_def ) :: qcl1p5m_ssi_flag logical( l_def ) :: rh1p5m_ssi_flag + logical( l_def ) :: t10m_ssi_flag + logical( l_def ) :: q10m_ssi_flag logical( l_def ) :: t1p5m_land_flag logical( l_def ) :: q1p5m_land_flag logical( l_def ) :: qcl1p5m_land_flag @@ -62,6 +64,8 @@ contains !> @param[in,out] q1p5m_ssi 1.5m specific humidity over sea and sea-ice !> @param[in,out] qcl1p5m_ssi 1.5m specific cloud water over sea and sea-ice !> @param[in,out] rh1p5m_ssi 1.5m relative humidity over sea and sea-ice + !> @param[in,out] t10m_ssi 10m temperature over sea and sea-ice + !> @param[in,out] q10m_ssi 10m specific humidity over sea and sea-ice !> @param[in,out] t1p5m_land 1.5m temperature over land !> @param[in,out] q1p5m_land 1.5m specific humidity over land !> @param[in,out] qcl1p5m_land 1.5m specific cloud water over land @@ -78,6 +82,7 @@ contains t1p5m, q1p5m, qcl1p5m, rh1p5m, & t1p5m_ssi, q1p5m_ssi, & qcl1p5m_ssi, rh1p5m_ssi, & + t10m_ssi, q10m_ssi, & t1p5m_land, q1p5m_land, & qcl1p5m_land, rh1p5m_land, & latent_heat, & @@ -92,6 +97,7 @@ contains type( field_type ), intent(inout) :: t1p5m_surft, q1p5m_surft type( field_type ), intent(inout) :: t1p5m, q1p5m, qcl1p5m, rh1p5m type( field_type ), intent(inout) :: t1p5m_ssi, q1p5m_ssi, qcl1p5m_ssi, rh1p5m_ssi + type( field_type ), intent(inout) :: t10m_ssi, q10m_ssi type( field_type ), intent(inout) :: t1p5m_land, q1p5m_land, qcl1p5m_land, rh1p5m_land type( field_type ), intent(inout) :: latent_heat type( field_type ), intent(inout) :: snomlt_surf_htf @@ -109,6 +115,13 @@ contains ! 2D fields rh1p5m_flag = init_diag(rh1p5m, 'surface__rh1p5m') rh1p5m_ssi_flag = init_diag(rh1p5m_ssi, 'surface__rh1p5m_ssi') + t10m_ssi_flag = init_diag(t10m_ssi, 'surface__t10m_ssi') + q10m_ssi_flag = init_diag(q10m_ssi, 'surface__q10m_ssi', activate=t10m_ssi_flag) + ! both diagnostics are required even if not output. + ! If q has been requested and t has not, ensure t is allocated + if (q10m_ssi_flag .and. .not. t10m_ssi_flag) then + t10m_ssi_flag = init_diag(t10m_ssi, 'surface__t10m_ssi', activate=q10m_ssi_flag) + end if rh1p5m_land_flag = init_diag(rh1p5m_land, 'surface__rh1p5m_land') snomlt_surf_htf_flag = init_diag(snomlt_surf_htf, 'surface__snomlt_surf_htf') soil_evap_flag = init_diag(soil_evap, 'surface__soil_evap') @@ -205,6 +218,8 @@ contains !> @param[in] q1p5m_ssi 1.5m specific humidity over sea and sea-ice !> @param[in] qcl1p5m_ssi 1.5m specific cloud water over sea and sea-ice !> @param[in] rh1p5m_ssi 1.5m relative humidity over sea and sea-ice + !> @param[in] t10m_ssi 10m temperature over sea and sea-ice + !> @param[in] q10m_ssi 10m specific humidity over sea and sea-ice !> @param[in] t1p5m_land 1.5m temperature over land !> @param[in] q1p5m_land 1.5m specific humidity over land !> @param[in] qcl1p5m_land 1.5m specific cloud water over land @@ -225,7 +240,7 @@ contains t1p5m_surft, q1p5m_surft, & t1p5m, q1p5m, qcl1p5m, rh1p5m, & t1p5m_ssi, q1p5m_ssi, qcl1p5m_ssi, & - rh1p5m_ssi, & + rh1p5m_ssi, t10m_ssi, q10m_ssi, & t1p5m_land, q1p5m_land, qcl1p5m_land, & rh1p5m_land, & latent_heat, & @@ -248,6 +263,7 @@ contains type( field_type ), intent(in) :: t1p5m_surft, q1p5m_surft type( field_type ), intent(in) :: t1p5m, q1p5m, qcl1p5m, rh1p5m type( field_type ), intent(in) :: t1p5m_ssi, q1p5m_ssi, qcl1p5m_ssi, rh1p5m_ssi + type( field_type ), intent(in) :: t10m_ssi, q10m_ssi type( field_type ), intent(in) :: t1p5m_land, q1p5m_land, qcl1p5m_land, rh1p5m_land type( field_type ), intent(in) :: latent_heat, snomlt_surf_htf type( field_type ), intent(in) :: soil_evap @@ -328,6 +344,8 @@ contains if (q1p5m_ssi_flag) call q1p5m_ssi%write_field() if (qcl1p5m_ssi_flag) call qcl1p5m_ssi%write_field() if (rh1p5m_ssi_flag) call rh1p5m_ssi%write_field() + if (t10m_ssi_flag) call t10m_ssi%write_field() + if (q10m_ssi_flag) call q10m_ssi%write_field() if (t1p5m_land_flag) call t1p5m_land%write_field() if (q1p5m_land_flag) call q1p5m_land%write_field() if (qcl1p5m_land_flag) call qcl1p5m_land%write_field() diff --git a/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 b/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 index 24e2aceb5e..e031c19e80 100644 --- a/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 +++ b/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 @@ -46,7 +46,7 @@ module jules_exp_kernel_mod !> type, public, extends(kernel_type) :: jules_exp_kernel_type private - type(arg_type) :: meta_args(108) = (/ & + type(arg_type) :: meta_args(109) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! theta_in_wth arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! exner_in_wth arg_type(GH_FIELD, GH_REAL, GH_READ, W3, STENCIL(REGION)), &! u_in_w3 @@ -154,6 +154,7 @@ module jules_exp_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! q1_sd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! diag__gross_prim_prod arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! diag__z0h_eff + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! diag__chr10m arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_1) &! ocn_cpl_point /) integer :: operates_on = DOMAIN @@ -276,7 +277,8 @@ module jules_exp_kernel_mod !> @param[in,out] q1_sd_2d StDev of level 1 humidity !> @param[in,out] gross_prim_prod Diagnostic: Gross Primary Productivity !> @param[in,out] z0h_eff Diagnostic: Gridbox mean effective roughness length for scalars - !> @param[in,out] ocn_cpl_point Diagnostic: Coupling point mask + !> @param[in,out] chr10m Diagnostic: 10m transfer coefficient + !> @param[in,out] ocn_cpl_point Coupling point mask !> @param[in] ndf_wth Number of DOFs per cell for potential temperature space !> @param[in] undf_wth Number of unique DOFs for potential temperature space !> @param[in] map_wth Dofmap for the cell at the base of the column for potential temperature space @@ -426,6 +428,7 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & q1_sd_2d, & gross_prim_prod, & z0h_eff, & + chr10m, & ocn_cpl_point, & ndf_wth, undf_wth, map_wth, & ndf_w3, undf_w3, map_w3, & @@ -696,6 +699,7 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & real(kind=r_def), dimension(undf_dust), intent(inout) :: dust_div_flux real(kind=r_def), pointer, intent(inout) :: z0h_eff(:), gross_prim_prod(:) + real(kind=r_def), pointer, intent(inout) :: chr10m(:) real(kind=r_def), intent(in) :: flux_h real(kind=r_def), intent(in) :: flux_e @@ -1393,6 +1397,8 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & ! needed to ensure z0h_eff is saved if wanted sf_diag%l_z0h_eff_gb = .not. associated(z0h_eff, empty_real_data) sf_diag%l_z0m_gb = .true. + sf_diag%l_t10m = .not. associated(chr10m, empty_real_data) + sf_diag%l_q10m = .not. associated(chr10m, empty_real_data) call alloc_sf_expl(sf_diag, .true., land_field) !----------------------------------------------------------------------- @@ -1855,6 +1861,11 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & z0h_eff(map_2d(1,i)) = sf_diag%z0h_eff_gb(i,1) end do end if + if (.not. associated(chr10m, empty_real_data) ) then + do i = 1, seg_len + chr10m(map_2d(1,i)) = sf_diag%chr10m(i,1) + end do + end if ! deallocate diagnostics call dealloc_sf_expl(sf_diag) diff --git a/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 b/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 index 618b6f307a..6e9f5e8186 100644 --- a/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 +++ b/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 @@ -38,7 +38,7 @@ module jules_imp_kernel_mod !> type, public, extends(kernel_type) :: jules_imp_kernel_type private - type(arg_type) :: meta_args(85) = (/ & + type(arg_type) :: meta_args(88) = (/ & arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! outer arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! loop arg_type(GH_FIELD, GH_REAL, GH_READ, W3), &! wetrho_in_w3 @@ -111,6 +111,9 @@ module jules_imp_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! q1p5m_ssi arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! qcl1p5m_ssi arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! rh1p5m_ssi + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1),&! chr10m + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! t10m_ssi + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! q10m_ssi arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! t1p5m_land arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! q1p5m_land arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! qcl1p5m_land @@ -209,6 +212,9 @@ module jules_imp_kernel_mod !> @param[in,out] q1p5m_ssi Diagnostic: 1.5m specific humidity over sea and sea-ice !> @param[in,out] qcl1p5m_ssi Diagnostic: 1.5m specific cloud water over sea and sea-ice !> @param[in,out] rh1p5m_ssi Diagnostic: 1.5m relative humidity over sea and sea-ice + !> @param[in] chr10m 10m transfer coefficient + !> @param[in,out] t10m_ssi Diagnostic: 10m temperature over sea and sea-ice + !> @param[in,out] q10m_ssi Diagnostic: 10m specific humidity over sea and sea-ice !> @param[in,out] t1p5m_land Diagnostic: 1.5m temperature over land !> @param[in,out] q1p5m_land Diagnostic: 1.5m specific humidity over land !> @param[in,out] qcl1p5m_land Diagnostic: 1.5m specific cloud water over land @@ -310,6 +316,7 @@ subroutine jules_imp_code(nlayers, seg_len, & t1p5m, q1p5m, qcl1p5m, rh1p5m, & t1p5m_ssi, q1p5m_ssi, & qcl1p5m_ssi, rh1p5m_ssi, & + chr10m, t10m_ssi, q10m_ssi, & t1p5m_land, q1p5m_land, & qcl1p5m_land, rh1p5m_land, & latent_heat, snomlt_surf_htf, & @@ -494,6 +501,7 @@ subroutine jules_imp_code(nlayers, seg_len, & real(kind=r_def), intent(in) :: dqw1_2d(undf_2d) real(kind=r_def), intent(in) :: dtl1_2d(undf_2d) real(kind=r_def), intent(in) :: ct_ctq1_2d(undf_2d) + real(kind=r_def), intent(in) :: chr10m(undf_2d) real(kind=r_def), intent(inout) :: surf_ht_flux(undf_tile) real(kind=r_def), pointer, intent(inout) :: t1p5m_surft(:) real(kind=r_def), pointer, intent(inout) :: q1p5m_surft(:) @@ -501,6 +509,7 @@ subroutine jules_imp_code(nlayers, seg_len, & real(kind=r_def), pointer, intent(inout) :: qcl1p5m(:), rh1p5m(:) real(kind=r_def), pointer, intent(inout) :: t1p5m_ssi(:), q1p5m_ssi(:) real(kind=r_def), pointer, intent(inout) :: qcl1p5m_ssi(:), rh1p5m_ssi(:) + real(kind=r_def), pointer, intent(inout) :: t10m_ssi(:), q10m_ssi(:) real(kind=r_def), pointer, intent(inout) :: t1p5m_land(:), q1p5m_land(:) real(kind=r_def), pointer, intent(inout) :: qcl1p5m_land(:), rh1p5m_land(:) real(kind=r_def), pointer, intent(inout) :: latent_heat(:) @@ -714,6 +723,8 @@ subroutine jules_imp_code(nlayers, seg_len, & .or. .not. associated(surf_lw_down, empty_real_data) sf_diag%l_lw_up_sice_weighted_cat = .not. associated(surf_lw_up, empty_real_data) sf_diag%sq1p5 = .true. + sf_diag%l_t10m = .not. associated(t10m_ssi, empty_real_data) + sf_diag%l_q10m = .not. associated(q10m_ssi, empty_real_data) call alloc_sf_imp(sf_diag, outer == outer_iterations, land_field) if (sf_diag%simlt) then do n = 1, nice @@ -722,6 +733,12 @@ subroutine jules_imp_code(nlayers, seg_len, & end do end do end if + if (sf_diag%l_t10m .or. sf_diag%l_q10m) then + allocate(sf_diag%chr10m(seg_len,1)) + do i = 1, seg_len + sf_diag%chr10m(i,1) = chr10m(map_2d(1,i)) + end do + end if allocate(epot_surft(land_field,ntiles)) allocate(tscrndcl_surft(land_field,ntiles)) @@ -1408,6 +1425,13 @@ subroutine jules_imp_code(nlayers, seg_len, & end do end do end if + if (sf_diag%l_t10m .or. sf_diag%l_q10m) then + ! Convert to mixing ratios + do i = 1, seg_len + sf_diag%q10m(i,1) = sf_diag%q10m(i,1) / & + (1.0_r_um+sf_diag%q10m(i,1)+qcf_latest(i,1)) + end do + end if do i = 1, seg_len ! Dummy values as unused in ls_cld for lowest level @@ -1489,6 +1513,27 @@ subroutine jules_imp_code(nlayers, seg_len, & end do end if + ! Sea and sea-ice 10m diagnostics + if (.not. associated(t10m_ssi, empty_real_data) .or. & + .not. associated(q10m_ssi, empty_real_data) ) then + + call ls_cld( & + forcing%pstar_ij, rhcpt, 1, 1, seg_len, 1, ntml, cumulus, & + .false., sf_diag%t10m, work_2d_1, sf_diag%q10m, & + qcf_latest, qcl1p5m_loc, work_2d_2, work_2d_3, error_code ) + end if + + if (.not. associated(t10m_ssi, empty_real_data) ) then + do i = 1, seg_len + t10m_ssi(map_2d(1,i)) = sf_diag%t10m(i,1) + end do + end if + if (.not. associated(q10m_ssi, empty_real_data) ) then + do i = 1, seg_len + q10m_ssi(map_2d(1,i)) = sf_diag%q10m(i,1) + end do + end if + ! Land screen level diagnostics if (.not. associated(t1p5m_land, empty_real_data) .or. & .not. associated(q1p5m_land, empty_real_data) .or. & diff --git a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 index a064523c8d..c900df9998 100644 --- a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 +++ b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 @@ -338,7 +338,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth &soil_moist_avail, snow_unload_rate, albedo_obs_scaling, soil_clay, soil_sand, dust_mrel, dust_flux, day_of_year, second_of_day, & flux_e, flux_h, urbwrr, urbhwr, urbhgt, urbztm, urbdisp, & &rhostar, recip_l_mo_sea, & -&t1_sd_2d, q1_sd_2d, gross_prim_prod, z0h_eff, ocn_cpl_point, stencil_depth) +&t1_sd_2d, q1_sd_2d, gross_prim_prod, z0h_eff, chr10m, ocn_cpl_point, stencil_depth) USE jules_exp_kernel_mod, ONLY: jules_exp_code USE mesh_mod, ONLY: mesh_type USE stencil_dofmap_mod, ONLY: STENCIL_REGION @@ -359,7 +359,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth &albedo_obs_scaling, soil_clay, soil_sand, dust_mrel, dust_flux, & urbwrr, urbhwr, urbhgt, urbztm, urbdisp, & rhostar, recip_l_mo_sea, t1_sd_2d, q1_sd_2d, & -&gross_prim_prod, z0h_eff +&gross_prim_prod, z0h_eff, chr10m TYPE(integer_field_type), intent(in) :: n_snow_layers, blend_height_tq, ocn_cpl_point INTEGER(KIND=i_def), intent(in) :: stencil_depth, ncells, ncells_halo, day_of_year, second_of_day REAL(KIND=r_def), intent(in) :: flux_e, flux_h @@ -386,7 +386,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth &soil_clay_proxy, soil_sand_proxy, dust_mrel_proxy, dust_flux_proxy, & urbwrr_proxy, urbhwr_proxy, urbhgt_proxy, urbztm_proxy, urbdisp_proxy, & rhostar_proxy, recip_l_mo_sea_proxy, & -&t1_sd_2d_proxy, q1_sd_2d_proxy, gross_prim_prod_proxy, z0h_eff_proxy +&t1_sd_2d_proxy, q1_sd_2d_proxy, gross_prim_prod_proxy, z0h_eff_proxy, chr10m_proxy INTEGER(KIND=i_def), pointer :: map_adspc10_dust_mrel(:,:) => null(), map_adspc1_zh(:,:) => null(), & &map_adspc2_tile_fraction(:,:) => null(), map_adspc3_leaf_area_index(:,:) => null(), & &map_adspc4_sea_ice_temperature(:,:) => null(), map_adspc5_snow_layer_thickness(:,:) => null(), & @@ -519,6 +519,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth q1_sd_2d_proxy = q1_sd_2d%get_proxy() gross_prim_prod_proxy = gross_prim_prod%get_proxy() z0h_eff_proxy = z0h_eff%get_proxy() + chr10m_proxy = chr10m%get_proxy() ocn_cpl_point_proxy = ocn_cpl_point%get_proxy() ! ! Initialise number of layers @@ -673,7 +674,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth urbdisp_proxy%data, & rhostar_proxy%data, recip_l_mo_sea_proxy%data, & &t1_sd_2d_proxy%data, q1_sd_2d_proxy%data, gross_prim_prod_proxy%data, & -z0h_eff_proxy%data, ocn_cpl_point_proxy%data, ndf_wtheta, & +z0h_eff_proxy%data, chr10m_proxy%data, ocn_cpl_point_proxy%data, ndf_wtheta, & &undf_wtheta, map_wtheta, ndf_w3, undf_w3, map_w3, ndf_adspc1_zh, undf_adspc1_zh, map_adspc1_zh, & &ndf_adspc2_tile_fraction, undf_adspc2_tile_fraction, map_adspc2_tile_fraction, ndf_adspc3_leaf_area_index, & &undf_adspc3_leaf_area_index, map_adspc3_leaf_area_index, ndf_adspc4_sea_ice_temperature, undf_adspc4_sea_ice_temperature, & @@ -726,6 +727,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth CALL q1_sd_2d_proxy%set_dirty() CALL gross_prim_prod_proxy%set_dirty() CALL z0h_eff_proxy%set_dirty() + CALL chr10m_proxy%set_dirty() ! ! END SUBROUTINE invoke_jules_exp_kernel_type From b240db0521475673b1834bbe243b39c6f082346a Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Mon, 27 Jul 2026 09:21:02 +0100 Subject: [PATCH 02/11] thermal speed diagnostic --- .../lfric_atm/metadata/field_def_diags.xml | 1 + .../diagnostics/jules_imp_diags_mod.x90 | 3 +- .../source/diagnostics/bl_imp_diags_mod.x90 | 8 ++++- .../kernel/bl_extra_diags_kernel_mod.F90 | 35 ++++++++++++++++--- .../file/file_def_diags_oper_nwp_gl.xml | 1 + 5 files changed, 42 insertions(+), 6 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index f0399a66eb..f7db93e4b7 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -659,6 +659,7 @@ + diff --git a/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 b/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 index 8a98a2c164..f23cc4df65 100644 --- a/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 +++ b/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 @@ -170,7 +170,8 @@ contains qcl1p5m_land_flag = init_diag(qcl1p5m_land, 'surface__qcl1p5m_land') ! Initialise fields where another diagnostic depends on them - if (samp_diag('surface__ustar_implicit') .or. & + if (samp_diag('turbulence__thermal_speed') .or. & + samp_diag('surface__ustar_implicit') .or. & samp_diag('surface__wind_gust') .or. & samp_diag('surface__scale_dep_wind_gust') .or. & samp_diag('surface__visibility_no_precip') .or. & diff --git a/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 index bf02ad998e..5e5b59331d 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 @@ -34,6 +34,7 @@ module bl_imp_diags_mod logical( l_def ) :: dw_bl_flag logical( l_def ) :: u10m_flag logical( l_def ) :: v10m_flag + logical( l_def ) :: thermal_speed_flag logical( l_def ) :: wind_gust_flag logical( l_def ) :: scale_dep_wind_gust_flag logical( l_def ) :: wspd10m_neut_flag @@ -265,6 +266,7 @@ contains type( field_type ), intent(inout) :: dcff_bl ! Local variables + type( field_type ) :: thermal_speed type( field_type ) :: wind_gust type( field_type ) :: scale_dep_wind_gust type( field_type ) :: visibility_with_precip @@ -391,6 +393,7 @@ contains if (ustar_implicit_flag) call ustar_implicit%write_field() ! Local diagnostic calculated here + thermal_speed_flag = init_diag(thermal_speed, 'turbulence__thermal_speed') wind_gust_flag = init_diag(wind_gust, 'surface__wind_gust') scale_dep_wind_gust_flag = init_diag(scale_dep_wind_gust, 'surface__scale_dep_wind_gust') vis_no_precip_flag = init_diag(visibility_no_precip, 'surface__visibility_no_precip') @@ -408,7 +411,8 @@ contains ignore = init_diag(visibility_no_precip, 'surface__visibility_no_precip', activate=.true.) end if - if ( wind_gust_flag .or. scale_dep_wind_gust_flag .or. & + if ( thermal_speed_flag .or. & + wind_gust_flag .or. scale_dep_wind_gust_flag .or. & vis_no_precip_flag .or. vis_with_precip_flag .or. & fog_fraction_flag .or. fog_fraction_ssi_flag .or. & fog_fraction_land_flag .or. vis_prob_5km_flag .or. dew_point_flag & @@ -429,6 +433,7 @@ contains ls_rain, ls_snow, lsca_2d, & conv_rain, conv_snow, cca_2d, & ustar_implicit, & + thermal_speed, & wind_gust, scale_dep_wind_gust, & fog_fraction, fog_fraction_ssi, & fog_fraction_land, vis_prob_5km,& @@ -437,6 +442,7 @@ contains visibility_with_precip, & visibility_no_precip ) ) + if (thermal_speed_flag) call thermal_speed%write_field() if (wind_gust_flag) call wind_gust%write_field() if (scale_dep_wind_gust_flag) call scale_dep_wind_gust%write_field() if (vis_no_precip_flag) call visibility_no_precip%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 index 323d262e84..d001ac1bcd 100644 --- a/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 @@ -26,7 +26,7 @@ module bl_extra_diags_kernel_mod !> type, public, extends(kernel_type) :: bl_extra_diags_kernel_type private - type(arg_type) :: meta_args(41) = (/ & + type(arg_type) :: meta_args(42) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! rho_in_w3 arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! wetrho_in_w3 arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! heat_flux_bl @@ -57,6 +57,7 @@ module bl_extra_diags_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & ! conv_snow_2d arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & ! cca_2d_in arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & ! ustar_implicit + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! thermal_speed arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! wind_gust arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! scale_dep_wind_gust arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! fog_fraction @@ -111,6 +112,7 @@ module bl_extra_diags_kernel_mod !> @param[in] conv_snow_2d Surface convective snowfall rate !> @param[in] cca_2d_in 2D convective cloud fraction !> @param[in] ustar_implicit Implicit friction velocity + !> @param[in,out] thermal_speed Thermal speed !> @param[in,out] wind_gust Wind gust !> @param[in,out] scale_dep_wind_gust Scale dependent wind gust !> @param[in,out] fog_fraction Fog_fraction @@ -151,7 +153,9 @@ subroutine bl_extra_diags_code( nlayers, & lsca_2d, & conv_rain_2d, & conv_snow_2d, cca_2d_in, & - ustar_implicit, wind_gust,& + ustar_implicit, & + thermal_speed, & + wind_gust, & scale_dep_wind_gust, & fog_fraction, & fog_fraction_ssi, & @@ -177,6 +181,7 @@ subroutine bl_extra_diags_code( nlayers, & use cloud_inputs_mod, only : rhcrit use dewpnt_mod, only : dewpnt use fog_fr_mod, only : fog_fr + use jules_surface_mod, only : min_ustar use mphys_constants_mod, only : mprog_min use nlsizes_namelist_mod, only : row_length, rows use planet_config_mod, only : p_zero, kappa, gravity, cp @@ -221,6 +226,7 @@ subroutine bl_extra_diags_code( nlayers, & real(kind=r_def), intent(in), pointer :: t1p5m_land(:), q1p5m_land(:), qcl1p5m_land(:) real(kind=r_def), intent(in), pointer :: wspd10m(:), z0m_eff(:) real(kind=r_def), intent(inout), pointer :: ustar_implicit(:) + real(kind=r_def), intent(inout), pointer :: thermal_speed(:) real(kind=r_def), intent(inout), pointer :: wind_gust(:), scale_dep_wind_gust(:) real(kind=r_def), intent(inout), pointer :: fog_fraction(:), vis_prob_5km(:) real(kind=r_def), intent(inout), pointer :: fog_fraction_ssi(:), fog_fraction_land(:) @@ -234,6 +240,8 @@ subroutine bl_extra_diags_code( nlayers, & ! Tunable parameters used in the calculation of the wind gust real(kind=r_def), parameter :: c_ws = 1.0_r_def/24.0_r_def real(kind=r_def), parameter :: gust_const = 2.29_r_def + ! Parameter used in the thermal speed calculation + real(kind=r_def), parameter :: a_stab = 1.5_r_def ! Switches needed for visibility calculations logical(l_def), parameter :: pct = .false. ! Cloud amounts are in % @@ -260,17 +268,36 @@ subroutine bl_extra_diags_code( nlayers, & ! Local scalars real(kind=r_def) :: ftl_surf, fqw_surf, & - wstar3_imp, std_dev, gust_contribution + wstar3_imp, std_dev, gust_contribution, z_on_l, f_stab integer(kind=i_def) :: k, icode, i,j if ( .not. associated(wind_gust, empty_real_data) .or. & - .not. associated(scale_dep_wind_gust, empty_real_data) ) then + .not. associated(scale_dep_wind_gust, empty_real_data) .or. & + .not. associated(thermal_speed, empty_real_data) ) then ftl_surf = heat_flux_bl(map_w3(1)) / cp fqw_surf = moist_flux_bl(map_w3(1)) wstar3_imp = zh(map_2d(1)) * gravity * ( ftl_surf/t1p5m(map_2d(1)) + & fqw_surf*c_virtual ) / & rho_in_w3(map_w3(1)) + end if + + if (.not. associated(thermal_speed, empty_real_data) ) then + if (wstar3_imp > tiny(1.0_r_def)) then + if (ustar_implicit(map_2d(1)) > min_ustar) then + z_on_l = vkman * wstar3_imp / ustar_implicit(map_2d(1))**3.0_r_def + f_stab = a_stab * z_on_l / (1.0_r_def + a_stab * z_on_l) + else + f_stab = 1.0_r_def + end if + thermal_speed(map_2d(1)) = f_stab * wstar3_imp**one_third + else + thermal_speed(map_2d(1)) = 0.0_r_def + end if + end if + + if ( .not. associated(wind_gust, empty_real_data) .or. & + .not. associated(scale_dep_wind_gust, empty_real_data) ) then if ( wstar3_imp > 0.0_r_def ) then ! Include the stability dependence std_dev = gust_const * ( ustar_implicit(map_2d(1))**3.0_r_def + & 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 fb5ee25716..bf3e5cf44b 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 @@ -141,6 +141,7 @@ + From 75d9ed6db222ff1aceba0d4720eca39691978924 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Mon, 27 Jul 2026 15:52:32 +0100 Subject: [PATCH 03/11] vis diags --- .../lfric_atm/metadata/field_def_diags.xml | 2 + .../diagnostics/jules_imp_diags_mod.x90 | 6 +- .../source/diagnostics/bl_imp_diags_mod.x90 | 13 +- .../kernel/bl_extra_diags_kernel_mod.F90 | 145 ++++++++++++++---- 4 files changed, 133 insertions(+), 33 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index f7db93e4b7..668032b683 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -555,6 +555,8 @@ + + diff --git a/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 b/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 index f23cc4df65..356453b881 100644 --- a/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 +++ b/interfaces/jules_interface/source/diagnostics/jules_imp_diags_mod.x90 @@ -184,14 +184,16 @@ contains if (.not. qcl1p5m_flag) ignore = init_diag(qcl1p5m, 'surface__qcl1p5m', activate=.true.) end if ! sea and sea-ice derived diags - if (samp_diag('surface__fog_fraction_ssi') .or. & + if (samp_diag('surface__visibility_with_precip_ssi') .or. & + samp_diag('surface__fog_fraction_ssi') .or. & samp_diag('surface__dew_point_ssi')) then if (.not. t1p5m_ssi_flag) ignore = init_diag(t1p5m_ssi, 'surface__t1p5m_ssi', activate=.true.) if (.not. q1p5m_ssi_flag) ignore = init_diag(q1p5m_ssi, 'surface__q1p5m_ssi', activate=.true.) if (.not. qcl1p5m_ssi_flag) ignore = init_diag(qcl1p5m_ssi, 'surface__qcl1p5m_ssi', activate=.true.) end if ! land derived diags - if (samp_diag('surface__fog_fraction_land') .or. & + if (samp_diag('surface__visibility_with_precip_land') .or. & + samp_diag('surface__fog_fraction_land') .or. & samp_diag('surface__dew_point_land')) then if (.not. t1p5m_land_flag) ignore = init_diag(t1p5m_land, 'surface__t1p5m_land', activate=.true.) if (.not. q1p5m_land_flag) ignore = init_diag(q1p5m_land, 'surface__q1p5m_land', activate=.true.) diff --git a/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 index 5e5b59331d..6be6501d94 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/bl_imp_diags_mod.x90 @@ -47,6 +47,8 @@ module bl_imp_diags_mod logical( l_def ) :: pseudotauy_flag logical( l_def ) :: vis_no_precip_flag logical( l_def ) :: vis_with_precip_flag + logical( l_def ) :: vis_with_precip_land_flag + logical( l_def ) :: vis_with_precip_ssi_flag logical( l_def ) :: fog_fraction_flag logical( l_def ) :: fog_fraction_ssi_flag logical( l_def ) :: fog_fraction_land_flag @@ -271,6 +273,8 @@ contains type( field_type ) :: scale_dep_wind_gust type( field_type ) :: visibility_with_precip type( field_type ) :: visibility_no_precip + type( field_type ) :: vis_with_precip_land + type( field_type ) :: vis_with_precip_ssi type( field_type ) :: fog_fraction type( field_type ) :: fog_fraction_ssi type( field_type ) :: fog_fraction_land @@ -398,6 +402,8 @@ contains scale_dep_wind_gust_flag = init_diag(scale_dep_wind_gust, 'surface__scale_dep_wind_gust') vis_no_precip_flag = init_diag(visibility_no_precip, 'surface__visibility_no_precip') vis_with_precip_flag = init_diag(visibility_with_precip, 'surface__visibility_with_precip') + vis_with_precip_land_flag = init_diag(vis_with_precip_land, 'surface__visibility_with_precip_land') + vis_with_precip_ssi_flag = init_diag(vis_with_precip_ssi, 'surface__visibility_with_precip_ssi') fog_fraction_flag = init_diag(fog_fraction, 'surface__fog_fraction') fog_fraction_ssi_flag = init_diag(fog_fraction_ssi, 'surface__fog_fraction_ssi') fog_fraction_land_flag = init_diag(fog_fraction_land, 'surface__fog_fraction_land') @@ -414,6 +420,7 @@ contains if ( thermal_speed_flag .or. & wind_gust_flag .or. scale_dep_wind_gust_flag .or. & vis_no_precip_flag .or. vis_with_precip_flag .or. & + vis_with_precip_land_flag .or. vis_with_precip_ssi_flag .or. & fog_fraction_flag .or. fog_fraction_ssi_flag .or. & fog_fraction_land_flag .or. vis_prob_5km_flag .or. dew_point_flag & .or. dew_point_ssi_flag .or. dew_point_land_flag ) then @@ -440,13 +447,17 @@ contains dew_point, dew_point_ssi, & dew_point_land, & visibility_with_precip, & - visibility_no_precip ) ) + visibility_no_precip, & + vis_with_precip_land, & + vis_with_precip_ssi ) ) if (thermal_speed_flag) call thermal_speed%write_field() if (wind_gust_flag) call wind_gust%write_field() if (scale_dep_wind_gust_flag) call scale_dep_wind_gust%write_field() if (vis_no_precip_flag) call visibility_no_precip%write_field() if (vis_with_precip_flag) call visibility_with_precip%write_field() + if (vis_with_precip_land_flag) call vis_with_precip_land%write_field() + if (vis_with_precip_ssi_flag) call vis_with_precip_ssi%write_field() if (fog_fraction_flag) call fog_fraction%write_field() if (fog_fraction_ssi_flag) call fog_fraction_ssi%write_field() if (fog_fraction_land_flag) call fog_fraction_land%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 index d001ac1bcd..745a854584 100644 --- a/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/bl_extra_diags_kernel_mod.F90 @@ -12,7 +12,7 @@ module bl_extra_diags_kernel_mod GH_READ, GH_WRITE, & CELL_COLUMN, & ANY_DISCONTINUOUS_SPACE_1 - use constants_mod, only : r_def, i_def, i_um, r_um, l_def + use constants_mod, only : r_def, i_def, i_um, r_um, l_def, rmdi use empty_data_mod, only : empty_real_data use fs_continuity_mod, only : Wtheta, W3 use kernel_mod, only : kernel_type @@ -26,7 +26,7 @@ module bl_extra_diags_kernel_mod !> type, public, extends(kernel_type) :: bl_extra_diags_kernel_type private - type(arg_type) :: meta_args(42) = (/ & + type(arg_type) :: meta_args(44) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! rho_in_w3 arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! wetrho_in_w3 arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! heat_flux_bl @@ -68,7 +68,9 @@ module bl_extra_diags_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! dew_point_ssi arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! dew_point_land arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! visibility_with_precip - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1) & ! visibility_no_precip + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! visibility_no_precip + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! vis_with_precip_land + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1) & ! vis_with_precip_ssi /) integer :: operates_on = CELL_COLUMN contains @@ -124,6 +126,8 @@ module bl_extra_diags_kernel_mod !> @param[in,out] dew_point_land Dew point temperature over land !> @param[in,out] visibility_with_precip Visibility with precip included !> @param[in,out] visibility_no_precip Visibility without including precip + !> @param[in,out] vis_with_precip_land Visibility with precip included over land + !> @param[in,out] vis_with_precip_ssi Visibility with precip included over sea/sea-ice !> @param[in] ndf_w3 Number of degrees of freedom per cell for density space !> @param[in] undf_w3 Number unique of degrees of freedom for density space !> @param[in] map_w3 Dofmap for the cell at the base of the column for density space @@ -165,6 +169,8 @@ subroutine bl_extra_diags_code( nlayers, & dew_point_land, & visibility_with_precip, & visibility_no_precip, & + vis_with_precip_land, & + vis_with_precip_ssi, & ndf_w3, & undf_w3, & map_w3, & @@ -234,6 +240,8 @@ subroutine bl_extra_diags_code( nlayers, & real(kind=r_def), intent(inout), pointer :: dew_point_ssi(:), dew_point_land(:) real(kind=r_def), intent(inout), pointer :: visibility_with_precip(:) real(kind=r_def), intent(inout), pointer :: visibility_no_precip(:) + real(kind=r_def), intent(inout), pointer :: vis_with_precip_land(:) + real(kind=r_def), intent(inout), pointer :: vis_with_precip_ssi(:) real(kind=r_def), parameter :: one_third = 1.0_r_def/3.0_r_def @@ -327,6 +335,8 @@ subroutine bl_extra_diags_code( nlayers, & ! map main input fields if (.not. associated(visibility_no_precip, empty_real_data) .or. & .not. associated(visibility_with_precip, empty_real_data) .or. & + .not. associated(vis_with_precip_land, empty_real_data) .or. & + .not. associated(vis_with_precip_ssi, empty_real_data) .or. & .not. associated(fog_fraction, empty_real_data) .or. & .not. associated(fog_fraction_ssi, empty_real_data) .or. & .not. associated(fog_fraction_land, empty_real_data) .or. & @@ -344,6 +354,38 @@ subroutine bl_extra_diags_code( nlayers, & qcl1p5m_loc(1,1) = qcl1p5m(map_2d(1)) end if + if (.not. associated(visibility_with_precip, empty_real_data) .or. & + .not. associated(vis_with_precip_land, empty_real_data) .or. & + .not. associated(vis_with_precip_ssi, empty_real_data)) then + + ! map additional input fields required for visibility with precipitation + ! level 1 rho + rho1(1,1) = wetrho_in_w3(map_w3(1)) + ! level 1 cloud ice mixing ratio + qcf1(1,1) = mci(map_wth(1) + 1) + ! level 1 rain mixing ratio + qrain1(1,1) = mr(map_wth(1) + 1) + ! surface rain and snow rates from large-scale microphysics + ls_rain(1,1) = ls_rain_2d(map_2d(1)) + ls_snow(1,1) = ls_snow_2d(map_2d(1)) + ! surface rain and snow rates from convection + conv_rain(1,1) = conv_rain_2d(map_2d(1)) + conv_snow(1,1) = conv_snow_2d(map_2d(1)) + ! cca_2d + cca_2d(1,1) = cca_2d_in(map_2d(1)) + ! prob of ls precip - just use existing rain area fraction + plsp(1,1) = lsca_2d(map_2d(1)) + + ! number prognostics used in the visibility calculation + ! We only copy these if casim is enabled. Otherwise they will + ! not be used. + if (microphysics_casim) then + rainnumber(1,1,1) = nr_mphys(map_wth(1) + 1) + snownumber(1,1,1) = ns_mphys(map_wth(1) + 1) + end if + + end if + ! Visibility if ( .not. associated(visibility_no_precip, empty_real_data) .or. & .not. associated(visibility_with_precip, empty_real_data) ) then @@ -357,31 +399,6 @@ subroutine bl_extra_diags_code( nlayers, & ! Visibility at 1.5 m including precipitation if ( .not. associated(visibility_with_precip, empty_real_data) ) then - ! map additional input fields - ! level 1 rho - rho1(1,1) = wetrho_in_w3(map_w3(1)) - ! level 1 cloud ice mixing ratio - qcf1(1,1) = mci(map_wth(1) + 1) - ! level 1 rain mixing ratio - qrain1(1,1) = mr(map_wth(1) + 1) - ! surface rain and snow rates from large-scale microphysics - ls_rain(1,1) = ls_rain_2d(map_2d(1)) - ls_snow(1,1) = ls_snow_2d(map_2d(1)) - ! surface rain and snow rates from convection - conv_rain(1,1) = conv_rain_2d(map_2d(1)) - conv_snow(1,1) = conv_snow_2d(map_2d(1)) - ! cca_2d - cca_2d(1,1) = cca_2d_in(map_2d(1)) - ! prob of ls precip - just use existing rain area fraction - plsp(1,1) = lsca_2d(map_2d(1)) - - ! number prognostics used in the visibility calculation - ! We only copy these if casim is enabled. Otherwise they will - ! not be used. - if (microphysics_casim) then - rainnumber(1,1,1) = nr_mphys(map_wth(1) + 1) - snownumber(1,1,1) = ns_mphys(map_wth(1) + 1) - end if call beta_precip( ls_rain, ls_snow, & conv_rain, conv_snow, qcf1, qrain1, & @@ -429,7 +446,8 @@ subroutine bl_extra_diags_code( nlayers, & end if ! sea and sea-ice diagnostics - if (.not. associated(fog_fraction_ssi, empty_real_data) .or. & + if (.not. associated(vis_with_precip_ssi) .or. & + .not. associated(fog_fraction_ssi, empty_real_data) .or. & .not. associated(dew_point_ssi, empty_real_data) ) then ! copy of screen variables t1p5m_loc(1,1) = t1p5m_ssi(map_2d(1)) @@ -437,6 +455,39 @@ subroutine bl_extra_diags_code( nlayers, & qcl1p5m_loc(1,1) = qcl1p5m_ssi(map_2d(1)) end if + ! Visibility + if ( .not. associated(vis_with_precip_ssi, empty_real_data) ) then + if (t1p5m_loc(1,1) == 0.0_r_def) then + ! If the sea/sea-ice temperature is zero, then we are not on sea/sea-ice, so set + ! the visibility to rmdi. + vis_with_precip_ssi(map_2d(1)) = rmdi + else + call visbty( & + ! inputs + p_star, t1p5m_loc, q1p5m_loc, qcl1p5m_loc, aerosol1, & + calc_prob_of_vis, rhcrit(1), murk_visibility, 1, & + ! output + vis_no_precip ) + + call beta_precip( ls_rain, ls_snow, & + conv_rain, conv_snow, qcf1, qrain1, & + rho1, t1p5m_loc, p_star, snownumber, rainnumber, & + plsp,cca_2d,pct,avg, & + 1, 1, 1, & + beta_ls_rain, beta_ls_snow, & + beta_c_rain, beta_c_snow ) + call vis_precip( vis_no_precip, & + plsp,cca_2d,pct, & + beta_ls_rain, beta_ls_snow, & + beta_c_rain, beta_c_snow, & + 1, 1, 1, & + vis,vis_ls_precip,vis_c_precip, & + icode ) + vis_with_precip_ssi(map_2d(1)) = vis(1,1) + end if + + end if ! any vis + if ( .not. associated(fog_fraction_ssi, empty_real_data) ) then do k = 1, 1 vis_threshold(1,1,1,k)=vis_thresh(k) @@ -458,7 +509,8 @@ subroutine bl_extra_diags_code( nlayers, & end if ! land diagnostics - if (.not. associated(fog_fraction_land, empty_real_data) .or. & + if (.not. associated(vis_with_precip_land, empty_real_data) .or. & + .not. associated(fog_fraction_land, empty_real_data) .or. & .not. associated(dew_point_land, empty_real_data) ) then ! copy of screen variables t1p5m_loc(1,1) = t1p5m_land(map_2d(1)) @@ -466,6 +518,39 @@ subroutine bl_extra_diags_code( nlayers, & qcl1p5m_loc(1,1) = qcl1p5m_land(map_2d(1)) end if + ! Visibility + if ( .not. associated(vis_with_precip_land, empty_real_data) ) then + if (t1p5m_loc(1,1) == 0.0_r_def) then + ! If the land temperature is zero, then we are not on land, so set + ! the visibility to rmdi. + vis_with_precip_land(map_2d(1)) = rmdi + else + call visbty( & + ! inputs + p_star, t1p5m_loc, q1p5m_loc, qcl1p5m_loc, aerosol1, & + calc_prob_of_vis, rhcrit(1), murk_visibility, 1, & + ! output + vis_no_precip ) + + call beta_precip( ls_rain, ls_snow, & + conv_rain, conv_snow, qcf1, qrain1, & + rho1, t1p5m_loc, p_star, snownumber, rainnumber, & + plsp,cca_2d,pct,avg, & + 1, 1, 1, & + beta_ls_rain, beta_ls_snow, & + beta_c_rain, beta_c_snow ) + call vis_precip( vis_no_precip, & + plsp,cca_2d,pct, & + beta_ls_rain, beta_ls_snow, & + beta_c_rain, beta_c_snow, & + 1, 1, 1, & + vis,vis_ls_precip,vis_c_precip, & + icode ) + vis_with_precip_land(map_2d(1)) = vis(1,1) + end if + + end if ! any vis + if ( .not. associated(fog_fraction_land, empty_real_data) ) then do k = 1, 1 vis_threshold(1,1,1,k)=vis_thresh(k) From 18c0bd502b389e46b472967e07f8f347bba499ca Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Mon, 27 Jul 2026 15:55:02 +0100 Subject: [PATCH 04/11] add to op diags file --- rose-stem/app/lfric_atm/file/file_def_diags_oper_nwp_gl.xml | 2 ++ 1 file changed, 2 insertions(+) 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 bf3e5cf44b..55c497cb90 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 @@ -26,6 +26,8 @@ + + From c6e9812b90afd0fa2184650ee0958f69b6368ab4 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Tue, 28 Jul 2026 10:07:03 +0100 Subject: [PATCH 05/11] cloud base 4p5 okta --- .../lfric_atm/metadata/field_def_diags.xml | 1 + .../source/diagnostics/cld_diags_mod.x90 | 8 +++++++ .../source/kernel/cld_diags_kernel_mod.F90 | 23 ++++++++++++++++++- .../file/file_def_diags_oper_nwp_gl.xml | 1 + 4 files changed, 32 insertions(+), 1 deletion(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 668032b683..be1a17405e 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -331,6 +331,7 @@ + diff --git a/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 index 6b6eb9a400..f1005cdb52 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 @@ -39,6 +39,7 @@ module cld_diags_mod cld_amount_maxrnd_flag, & ceil_cld_amount_maxrnd_flag, & cld_base_altitude_flag, & + cld_base_4p5_okta_flag, & low_cld_base_altitude_flag, & very_low_cld_amount_flag, & low_cld_amount_flag, & @@ -124,6 +125,7 @@ contains type( field_type ) :: cld_amount_maxrnd type( field_type ) :: ceil_cld_amount_maxrnd type( field_type ) :: cld_base_altitude + type( field_type ) :: cld_base_4p5_okta type( field_type ) :: low_cld_base_altitude type( field_type ) :: very_low_cld_amount type( field_type ) :: low_cld_amount @@ -156,6 +158,7 @@ contains cld_amount_maxrnd_flag = init_diag(cld_amount_maxrnd, 'cloud__cloud_amount_maxrnd') ceil_cld_amount_maxrnd_flag = init_diag(ceil_cld_amount_maxrnd, 'cloud__ceilometer_cloud_amount_maxrnd') cld_base_altitude_flag = init_diag(cld_base_altitude, 'cloud__cloud_base_altitude') + cld_base_4p5_okta_flag = init_diag(cld_base_4p5_okta, 'cloud__cloud_base_4p5_okta') low_cld_base_altitude_flag = init_diag(low_cld_base_altitude, 'cloud__low_cloud_base_altitude') very_low_cld_amount_flag = init_diag(very_low_cld_amount, 'cloud__very_low_type_cloud_amount') low_cld_amount_flag = init_diag(low_cld_amount, 'cloud__low_type_cloud_amount') @@ -176,6 +179,7 @@ contains cld_amount_maxrnd_flag .or. & ceil_cld_amount_maxrnd_flag .or. & cld_base_altitude_flag .or. & + cld_base_4p5_okta_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & low_cld_amount_flag .or. & @@ -196,6 +200,7 @@ contains cld_amount_maxrnd_flag .or. & ceil_cld_amount_maxrnd_flag .or. & cld_base_altitude_flag .or. & + cld_base_4p5_okta_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & low_cld_amount_flag .or. & @@ -235,6 +240,7 @@ contains cld_amount_maxrnd_flag .or. & ceil_cld_amount_maxrnd_flag .or. & cld_base_altitude_flag .or. & + cld_base_4p5_okta_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & low_cld_amount_flag .or. & @@ -256,6 +262,7 @@ contains cld_amount_maxrnd, & ceil_cld_amount_maxrnd, & cld_base_altitude, & + cld_base_4p5_okta, & low_cld_base_altitude, & very_low_cld_amount, & low_cld_amount, & @@ -355,6 +362,7 @@ contains if (cld_amount_maxrnd_flag) call cld_amount_maxrnd%write_field() if (ceil_cld_amount_maxrnd_flag) call ceil_cld_amount_maxrnd%write_field() if (cld_base_altitude_flag) call cld_base_altitude%write_field() + if (cld_base_4p5_okta_flag) call cld_base_4p5_okta%write_field() if (low_cld_base_altitude_flag) call low_cld_base_altitude%write_field() if (very_low_cld_amount_flag) call very_low_cld_amount%write_field() if (low_cld_amount_flag) call low_cld_amount%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 index 9083258ecf..4431211bea 100644 --- a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 @@ -26,13 +26,14 @@ module cld_diags_kernel_mod !> type, public, extends(kernel_type) :: cld_diags_kernel_type private - type(arg_type) :: meta_args(20) = (/ & + type(arg_type) :: meta_args(21) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! combined_cld_amount_wth arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_max arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_rnd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_maxrnd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! ceil_cld_amount_maxrnd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_altitude + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_4p5_okta arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! low_cld_base_altitude arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! very_low_cld_amount arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! low_cld_amount @@ -67,6 +68,7 @@ module cld_diags_kernel_mod !> @param[in,out] cld_amount_maxrnd Cloud amount maximum-random overlap !> @param[in,out] ceil_cld_amount_maxrnd Ceilometer filtered cloud amount maximum-random overlap !> @param[in,out] cld_base_altitude Cloud base altitude wrt sea level + !> @param[in,out] cld_base_4p5_okta Cloud base altitude for 4.5 okta cloud amount !> @param[in,out] low_cld_base_altitude Cloud base altitude wrt sea level for very low amount of cloud !> @param[in,out] very_low_cld_amount Maximum cloud amount below 111m !> @param[in,out] low_cld_amount Maximum cloud amount between 111 and 1949m above sea level @@ -102,6 +104,7 @@ subroutine cld_diags_code( nlayers, & cld_amount_maxrnd, & ceil_cld_amount_maxrnd, & cld_base_altitude, & + cld_base_4p5_okta, & low_cld_base_altitude, & very_low_cld_amount, & low_cld_amount, & @@ -145,6 +148,7 @@ subroutine cld_diags_code( nlayers, & real(kind=r_def), pointer, intent(inout) :: cld_amount_maxrnd(:) real(kind=r_def), pointer, intent(inout) :: ceil_cld_amount_maxrnd(:) real(kind=r_def), pointer, intent(inout) :: cld_base_altitude(:) + real(kind=r_def), pointer, intent(inout) :: cld_base_4p5_okta(:) real(kind=r_def), pointer, intent(inout) :: low_cld_base_altitude(:) real(kind=r_def), pointer, intent(inout) :: very_low_cld_amount(:) real(kind=r_def), pointer, intent(inout) :: low_cld_amount(:) @@ -183,6 +187,7 @@ subroutine cld_diags_code( nlayers, & real(kind=r_def), parameter :: high_to_very_high = 13608.0_r_def ! metres ! When looking for cloud base: how much cloud cover defines cloud base. real(kind=r_def), parameter :: cld_cover_for_cld_base = 2.5_r_def/8.0_r_def + real(kind=r_def), parameter :: cld_cover_for_4p5_okta = 4.5_r_def/8.0_r_def ! When looking for cloud base sometimes want much smaller amount of cloud. real(kind=r_def), parameter :: low_cld_cover_for_cld_base = 0.05_r_def ! For tracking whether cloud has been found. @@ -306,6 +311,7 @@ subroutine cld_diags_code( nlayers, & ! Find heights above sea level (asl) if required. if (.not. associated(cld_base_altitude, empty_real_data) .or. & + .not. associated(cld_base_4p5_okta, empty_real_data) .or. & .not. associated(low_cld_base_altitude, empty_real_data) .or. & .not. associated(very_low_cld_amount, empty_real_data) .or. & .not. associated(low_cld_amount, empty_real_data) .or. & @@ -338,6 +344,21 @@ subroutine cld_diags_code( nlayers, & end do end if + ! cld_base_4p5_okta (in kilofeet) + if (.not. associated(cld_base_4p5_okta, empty_real_data) ) then + ! As a default, set cloud-base to beyond top of model + cld_base_4p5_okta(map_2d(1)) = 1.1_r_def * & + z_asl_centre_of_levels(nlayers) * & + m_to_kfeet + + do k = 1, nlayers + if ( combined_cld_amount(k) >= cld_cover_for_4p5_okta ) then + cld_base_4p5_okta(map_2d(1)) = z_asl_base_of_levels(k) * m_to_kfeet + exit + end if + end do + end if + ! low_cld_base_altitude (in feet) for low threshold cloud amounts if (.not. associated(low_cld_base_altitude, empty_real_data) ) then ! As a default, set cloud-base to beyond top of model 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 55c497cb90..ef3b43125e 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 @@ -58,6 +58,7 @@ + From 1d6e6854c525a5aba79742d8a82945e3cbb6cf39 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Tue, 28 Jul 2026 12:21:06 +0100 Subject: [PATCH 06/11] bugfixes --- .../source/diagnostics/bl_extra_diags_mod.x90 | 12 ++++++------ .../source/kernel/cld_diags_kernel_mod.F90 | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/diagnostics/bl_extra_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/bl_extra_diags_mod.x90 index 0ce369f485..e66226d8a4 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/bl_extra_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/bl_extra_diags_mod.x90 @@ -91,9 +91,9 @@ contains type( field_type ), pointer :: murk ! Local variables + type( field_type ) :: thermal_speed type( field_type ) :: wind_gust type( field_type ) :: scale_dep_wind_gust - type( field_type ) :: thermal_speed type( field_type ) :: visibility_with_precip type( field_type ) :: visibility_no_precip type( field_type ) :: vis_with_precip_land @@ -107,9 +107,9 @@ contains type( field_type ) :: dew_point_land ! Logical indicating whether diagnostics are requested + logical( l_def ) :: thermal_speed_flag logical( l_def ) :: wind_gust_flag logical( l_def ) :: scale_dep_wind_gust_flag - logical( l_def ) :: thermal_speed_flag logical( l_def ) :: vis_no_precip_flag logical( l_def ) :: vis_with_precip_flag logical( l_def ) :: vis_with_precip_land_flag @@ -148,9 +148,9 @@ contains if (samp_diag('aerosol__murk')) call murk%write_field('aerosol__murk') ! Local diagnostic calculated here + thermal_speed_flag = init_diag(thermal_speed, 'turbulence__thermal_speed') wind_gust_flag = init_diag(wind_gust, 'surface__wind_gust') scale_dep_wind_gust_flag = init_diag(scale_dep_wind_gust, 'surface__scale_dep_wind_gust') - thermal_speed_flag = init_diag(thermal_speed, 'surface__thermal_speed') vis_no_precip_flag = init_diag(visibility_no_precip, 'surface__visibility_no_precip') vis_with_precip_flag = init_diag(visibility_with_precip, 'surface__visibility_with_precip') vis_with_precip_land_flag = init_diag(vis_with_precip_land, 'surface__visibility_with_precip_land') @@ -168,8 +168,8 @@ contains ignore = init_diag(visibility_no_precip, 'surface__visibility_no_precip', activate=.true.) end if - if ( wind_gust_flag .or. scale_dep_wind_gust_flag .or. & - thermal_speed_flag .or. & + if ( thermal_speed_flag .or. & + wind_gust_flag .or. scale_dep_wind_gust_flag .or. & vis_no_precip_flag .or. vis_with_precip_flag .or. & vis_with_precip_land_flag .or. vis_with_precip_ssi_flag .or. & fog_fraction_flag .or. fog_fraction_ssi_flag .or. & @@ -202,9 +202,9 @@ contains vis_with_precip_land, & vis_with_precip_ssi ) ) + if (thermal_speed_flag) call thermal_speed%write_field() if (wind_gust_flag) call wind_gust%write_field() if (scale_dep_wind_gust_flag) call scale_dep_wind_gust%write_field() - if (thermal_speed_flag) call thermal_speed%write_field() if (vis_no_precip_flag) call visibility_no_precip%write_field() if (vis_with_precip_flag) call visibility_with_precip%write_field() if (vis_with_precip_land_flag) call vis_with_precip_land%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 index ea48c80289..4bacc3e4a4 100644 --- a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 @@ -26,7 +26,7 @@ module cld_diags_kernel_mod !> type, public, extends(kernel_type) :: cld_diags_kernel_type private - type(arg_type) :: meta_args(21) = (/ & + type(arg_type) :: meta_args(22) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! combined_cld_amount_wth arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_max arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_rnd From e57a6d28e829ff05a5ecfcf745b88c829463e8b3 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Thu, 6 Aug 2026 09:32:41 +0100 Subject: [PATCH 07/11] fix comment --- .../jules_interface/source/kernel/jules_imp_kernel_mod.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 b/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 index 6e9f5e8186..228548a0fa 100644 --- a/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 +++ b/interfaces/jules_interface/source/kernel/jules_imp_kernel_mod.F90 @@ -1426,7 +1426,6 @@ subroutine jules_imp_code(nlayers, seg_len, & end do end if if (sf_diag%l_t10m .or. sf_diag%l_q10m) then - ! Convert to mixing ratios do i = 1, seg_len sf_diag%q10m(i,1) = sf_diag%q10m(i,1) / & (1.0_r_um+sf_diag%q10m(i,1)+qcf_latest(i,1)) From c01dff0ddd5e4de5195658477a6197b81d411e0f Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Wed, 12 Aug 2026 10:11:52 +0000 Subject: [PATCH 08/11] sd_orog diagnostic Output the standard deviation of orography ancillary (UM STASH 6-203) from the orographic drag diagnostics. Some of the content of this change has been produced with the assistance of Anthropic Claude Opus 5 (Claude Code). --- applications/lfric_atm/metadata/field_def_diags.xml | 1 + .../source/algorithm/orographic_drag_alg_mod.x90 | 6 +++++- .../diagnostics/orographic_drag_diags_mod.x90 | 13 +++++++++++-- 3 files changed, 17 insertions(+), 3 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 0e3cf00d9f..882af5bc94 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -498,6 +498,7 @@ orographic_drag__tauy_orog_blk + orographic_drag__tauy_orog_gwd + diff --git a/interfaces/physics_schemes_interface/source/algorithm/orographic_drag_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/orographic_drag_alg_mod.x90 index 8964530ad5..511b3d4963 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/orographic_drag_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/orographic_drag_alg_mod.x90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- !> @brief Interface to the orographic drag parametrization scheme module orographic_drag_alg_mod @@ -191,7 +194,8 @@ contains du_orographic_drag, dv_orographic_drag, dtemp_orographic_drag,& dtheta_orographic_drag, du_orog_blk, dv_orog_blk, & dtemp_orog_blk, du_orog_gwd, dv_orog_gwd, dtemp_orog_gwd, & - taux_orog_blk, tauy_orog_blk, taux_orog_gwd, tauy_orog_gwd) + taux_orog_blk, tauy_orog_blk, taux_orog_gwd, tauy_orog_gwd, & + sd_orog) end if diff --git a/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 index 5bb96cac45..7970ad795f 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- !> @brief Processes diagnostics for orographic_drag_alg module orographic_drag_diags_mod @@ -70,7 +73,8 @@ contains du_orographic_drag, dv_orographic_drag, dtemp_orographic_drag, & dtheta_orographic_drag, du_orog_blk, dv_orog_blk, dtemp_orog_blk, & du_orog_gwd, dv_orog_gwd, dtemp_orog_gwd, & - taux_orog_blk, tauy_orog_blk, taux_orog_gwd, tauy_orog_gwd) + taux_orog_blk, tauy_orog_blk, taux_orog_gwd, tauy_orog_gwd, & + sd_orog) implicit none @@ -78,7 +82,8 @@ contains type( field_type ), intent(in) :: du_orographic_drag, dv_orographic_drag, dtemp_orographic_drag, & dtheta_orographic_drag, du_orog_blk, dv_orog_blk, dtemp_orog_blk, & du_orog_gwd, dv_orog_gwd, dtemp_orog_gwd, & - taux_orog_blk, tauy_orog_blk, taux_orog_gwd, tauy_orog_gwd + taux_orog_blk, tauy_orog_blk, taux_orog_gwd, tauy_orog_gwd, & + sd_orog integer( tik ) :: id if ( LPROF ) call start_timing( id, 'diags.orog_gwd' ) @@ -95,6 +100,10 @@ contains call dv_orog_gwd%write_field('orographic_drag__dv_orog_gwd') call dtemp_orog_gwd%write_field('orographic_drag__dtemp_orog_gwd') + ! Ancillary from the orography collection. It is time-constant, so an + ! output stream only needs it once. + call sd_orog%write_field('orographic_drag__sd_orog') + ! Diagnostics computed within the kernel if (taux_orog_blk_flag) call taux_orog_blk%write_field() if (tauy_orog_blk_flag) call tauy_orog_blk%write_field() From 40c02716c423556a1e593618272d4fa055575e79 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Wed, 12 Aug 2026 11:29:13 +0100 Subject: [PATCH 09/11] fix comment --- .../source/diagnostics/orographic_drag_diags_mod.x90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 index 7970ad795f..f24ed3514c 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/orographic_drag_diags_mod.x90 @@ -100,8 +100,7 @@ contains call dv_orog_gwd%write_field('orographic_drag__dv_orog_gwd') call dtemp_orog_gwd%write_field('orographic_drag__dtemp_orog_gwd') - ! Ancillary from the orography collection. It is time-constant, so an - ! output stream only needs it once. + ! Ancillary from the orography collection. call sd_orog%write_field('orographic_drag__sd_orog') ! Diagnostics computed within the kernel From bc0762c0d4486b8bb077a0be1c7c0aae10e20a3f Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Fri, 14 Aug 2026 09:33:47 +0100 Subject: [PATCH 10/11] 7p9 okta cloud base --- .../lfric_atm/metadata/field_def_diags.xml | 1 + .../source/diagnostics/cld_diags_mod.x90 | 8 +++++++ .../source/kernel/cld_diags_kernel_mod.F90 | 23 ++++++++++++++++++- 3 files changed, 31 insertions(+), 1 deletion(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 882af5bc94..7c2d9910c7 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -335,6 +335,7 @@ + diff --git a/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 index 7f9399d182..475596801e 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 @@ -42,6 +42,7 @@ module cld_diags_mod ceil_cld_amount_maxrnd_flag, & cld_base_altitude_flag, & cld_base_4p5_okta_flag, & + cld_base_7p9_okta_flag, & cld_top_altitude_flag, & low_cld_base_altitude_flag, & very_low_cld_amount_flag, & @@ -132,6 +133,7 @@ contains type( field_type ) :: ceil_cld_amount_maxrnd type( field_type ) :: cld_base_altitude type( field_type ) :: cld_base_4p5_okta + type( field_type ) :: cld_base_7p9_okta type( field_type ) :: cld_top_altitude type( field_type ) :: low_cld_base_altitude type( field_type ) :: very_low_cld_amount @@ -166,6 +168,7 @@ contains ceil_cld_amount_maxrnd_flag = init_diag(ceil_cld_amount_maxrnd, 'cloud__ceilometer_cloud_amount_maxrnd') cld_base_altitude_flag = init_diag(cld_base_altitude, 'cloud__cloud_base_altitude') cld_base_4p5_okta_flag = init_diag(cld_base_4p5_okta, 'cloud__cloud_base_4p5_okta') + cld_base_7p9_okta_flag = init_diag(cld_base_7p9_okta, 'cloud__cloud_base_7p9_okta') cld_top_altitude_flag = init_diag(cld_top_altitude, 'cloud__cloud_top_altitude') low_cld_base_altitude_flag = init_diag(low_cld_base_altitude, 'cloud__low_cloud_base_altitude') very_low_cld_amount_flag = init_diag(very_low_cld_amount, 'cloud__very_low_type_cloud_amount') @@ -188,6 +191,7 @@ contains ceil_cld_amount_maxrnd_flag .or. & cld_base_altitude_flag .or. & cld_base_4p5_okta_flag .or. & + cld_base_7p9_okta_flag .or. & cld_top_altitude_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & @@ -210,6 +214,7 @@ contains ceil_cld_amount_maxrnd_flag .or. & cld_base_altitude_flag .or. & cld_base_4p5_okta_flag .or. & + cld_base_7p9_okta_flag .or. & cld_top_altitude_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & @@ -251,6 +256,7 @@ contains ceil_cld_amount_maxrnd_flag .or. & cld_base_altitude_flag .or. & cld_base_4p5_okta_flag .or. & + cld_base_7p9_okta_flag .or. & cld_top_altitude_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & @@ -275,6 +281,7 @@ contains ceil_cld_amount_maxrnd, & cld_base_altitude, & cld_base_4p5_okta, & + cld_base_7p9_okta, & cld_top_altitude, & low_cld_base_altitude, & very_low_cld_amount, & @@ -376,6 +383,7 @@ contains if (ceil_cld_amount_maxrnd_flag) call ceil_cld_amount_maxrnd%write_field() if (cld_base_altitude_flag) call cld_base_altitude%write_field() if (cld_base_4p5_okta_flag) call cld_base_4p5_okta%write_field() + if (cld_base_7p9_okta_flag) call cld_base_7p9_okta%write_field() if (cld_top_altitude_flag) call cld_top_altitude%write_field() if (low_cld_base_altitude_flag) call low_cld_base_altitude%write_field() if (very_low_cld_amount_flag) call very_low_cld_amount%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 index 4bacc3e4a4..d8353d332d 100644 --- a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 @@ -26,7 +26,7 @@ module cld_diags_kernel_mod !> type, public, extends(kernel_type) :: cld_diags_kernel_type private - type(arg_type) :: meta_args(22) = (/ & + type(arg_type) :: meta_args(23) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! combined_cld_amount_wth arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_max arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_rnd @@ -34,6 +34,7 @@ module cld_diags_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! ceil_cld_amount_maxrnd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_altitude arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_4p5_okta + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_7p9_okta arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_top_altitude arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! low_cld_base_altitude arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! very_low_cld_amount @@ -70,6 +71,7 @@ module cld_diags_kernel_mod !> @param[in,out] ceil_cld_amount_maxrnd Ceilometer filtered cloud amount maximum-random overlap !> @param[in,out] cld_base_altitude Cloud base altitude wrt sea level !> @param[in,out] cld_base_4p5_okta Cloud base altitude for 4.5 okta cloud amount + !> @param[in,out] cld_base_7p9_okta Cloud base altitude for 7.9 okta cloud amount !> @param[in,out] cld_top_altitude Cloud top altitude wrt sea level !> @param[in,out] low_cld_base_altitude Cloud base altitude wrt sea level for very low amount of cloud !> @param[in,out] very_low_cld_amount Maximum cloud amount below 111m @@ -107,6 +109,7 @@ subroutine cld_diags_code( nlayers, & ceil_cld_amount_maxrnd, & cld_base_altitude, & cld_base_4p5_okta, & + cld_base_7p9_okta, & cld_top_altitude, & low_cld_base_altitude, & very_low_cld_amount, & @@ -153,6 +156,7 @@ subroutine cld_diags_code( nlayers, & real(kind=r_def), pointer, intent(inout) :: ceil_cld_amount_maxrnd(:) real(kind=r_def), pointer, intent(inout) :: cld_base_altitude(:) real(kind=r_def), pointer, intent(inout) :: cld_base_4p5_okta(:) + real(kind=r_def), pointer, intent(inout) :: cld_base_7p9_okta(:) real(kind=r_def), pointer, intent(inout) :: cld_top_altitude(:) real(kind=r_def), pointer, intent(inout) :: low_cld_base_altitude(:) real(kind=r_def), pointer, intent(inout) :: very_low_cld_amount(:) @@ -193,6 +197,7 @@ subroutine cld_diags_code( nlayers, & ! When looking for cloud base or top: how much cloud cover defines cloud boundary. real(kind=r_def), parameter :: cld_cover_for_cld_bdry = 2.5_r_def/8.0_r_def real(kind=r_def), parameter :: cld_cover_for_4p5_okta = 4.5_r_def/8.0_r_def + real(kind=r_def), parameter :: cld_cover_for_7p9_okta = 7.9_r_def/8.0_r_def ! When looking for cloud base sometimes want much smaller amount of cloud. real(kind=r_def), parameter :: low_cld_cover_for_cld_base = 0.05_r_def ! For tracking whether cloud has been found. @@ -317,6 +322,7 @@ subroutine cld_diags_code( nlayers, & ! Find heights above sea level (asl) if required. if (.not. associated(cld_base_altitude, empty_real_data) .or. & .not. associated(cld_base_4p5_okta, empty_real_data) .or. & + .not. associated(cld_base_7p9_okta, empty_real_data) .or. & .not. associated(cld_top_altitude, empty_real_data) .or. & .not. associated(low_cld_base_altitude, empty_real_data) .or. & .not. associated(very_low_cld_amount, empty_real_data) .or. & @@ -365,6 +371,21 @@ subroutine cld_diags_code( nlayers, & end do end if + ! cld_base_7p9_okta (in kilofeet) + if (.not. associated(cld_base_7p9_okta, empty_real_data) ) then + ! As a default, set cloud-base to beyond top of model + cld_base_7p9_okta(map_2d(1)) = 1.1_r_def * & + z_asl_centre_of_levels(nlayers) * & + m_to_kfeet + + do k = 1, nlayers + if ( combined_cld_amount(k) >= cld_cover_for_7p9_okta ) then + cld_base_7p9_okta(map_2d(1)) = z_asl_base_of_levels(k) * m_to_kfeet + exit + end if + end do + end if + ! cld_top_altitude (in kilofeet) if (.not. associated(cld_top_altitude, empty_real_data) ) then ! As a default, set cloud-top to missing data From 6fae3e0e6fd4ee360d93c3ec3b9b306ffe09f849 Mon Sep 17 00:00:00 2001 From: iboutle <135141261+iboutle@users.noreply.github.com> Date: Mon, 24 Aug 2026 10:50:11 +0100 Subject: [PATCH 11/11] revised cloud base diags --- .../lfric_atm/metadata/field_def_diags.xml | 5 +- .../source/diagnostics/cld_diags_mod.x90 | 40 +++++++----- .../source/kernel/cld_diags_kernel_mod.F90 | 65 ++++++++++++------- .../file/file_def_diags_idealised.xml | 2 +- .../file/file_def_diags_idealised1.xml | 2 +- .../file/file_def_diags_ls_and_jedi.xml | 2 +- .../app/lfric_atm/file/file_def_diags_ral.xml | 2 +- .../app/lfric_atm/file/file_def_diags_ver.xml | 4 +- .../file/iodef_basic_gal.xml | 2 +- 9 files changed, 77 insertions(+), 47 deletions(-) diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index 7c2d9910c7..a410f74caa 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -333,10 +333,11 @@ - + + - + diff --git a/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 b/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 index 475596801e..14e4ae55a6 100644 --- a/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 +++ b/interfaces/physics_schemes_interface/source/diagnostics/cld_diags_mod.x90 @@ -40,10 +40,11 @@ module cld_diags_mod cld_amount_rnd_flag, & cld_amount_maxrnd_flag, & ceil_cld_amount_maxrnd_flag, & - cld_base_altitude_flag, & + cld_base_0p1_okta_flag, & + cld_base_2p5_okta_flag, & cld_base_4p5_okta_flag, & cld_base_7p9_okta_flag, & - cld_top_altitude_flag, & + cld_top_2p5_okta_flag, & low_cld_base_altitude_flag, & very_low_cld_amount_flag, & low_cld_amount_flag, & @@ -131,10 +132,11 @@ contains type( field_type ) :: cld_amount_rnd type( field_type ) :: cld_amount_maxrnd type( field_type ) :: ceil_cld_amount_maxrnd - type( field_type ) :: cld_base_altitude + type( field_type ) :: cld_base_0p1_okta + type( field_type ) :: cld_base_2p5_okta type( field_type ) :: cld_base_4p5_okta type( field_type ) :: cld_base_7p9_okta - type( field_type ) :: cld_top_altitude + type( field_type ) :: cld_top_2p5_okta type( field_type ) :: low_cld_base_altitude type( field_type ) :: very_low_cld_amount type( field_type ) :: low_cld_amount @@ -166,10 +168,11 @@ contains cld_amount_rnd_flag = init_diag(cld_amount_rnd, 'cloud__cloud_amount_rnd') cld_amount_maxrnd_flag = init_diag(cld_amount_maxrnd, 'cloud__cloud_amount_maxrnd') ceil_cld_amount_maxrnd_flag = init_diag(ceil_cld_amount_maxrnd, 'cloud__ceilometer_cloud_amount_maxrnd') - cld_base_altitude_flag = init_diag(cld_base_altitude, 'cloud__cloud_base_altitude') + cld_base_0p1_okta_flag = init_diag(cld_base_0p1_okta, 'cloud__cloud_base_0p1_okta') + cld_base_2p5_okta_flag = init_diag(cld_base_2p5_okta, 'cloud__cloud_base_2p5_okta') cld_base_4p5_okta_flag = init_diag(cld_base_4p5_okta, 'cloud__cloud_base_4p5_okta') cld_base_7p9_okta_flag = init_diag(cld_base_7p9_okta, 'cloud__cloud_base_7p9_okta') - cld_top_altitude_flag = init_diag(cld_top_altitude, 'cloud__cloud_top_altitude') + cld_top_2p5_okta_flag = init_diag(cld_top_2p5_okta, 'cloud__cloud_top_2p5_okta') low_cld_base_altitude_flag = init_diag(low_cld_base_altitude, 'cloud__low_cloud_base_altitude') very_low_cld_amount_flag = init_diag(very_low_cld_amount, 'cloud__very_low_type_cloud_amount') low_cld_amount_flag = init_diag(low_cld_amount, 'cloud__low_type_cloud_amount') @@ -189,10 +192,11 @@ contains cld_amount_rnd_flag .or. & cld_amount_maxrnd_flag .or. & ceil_cld_amount_maxrnd_flag .or. & - cld_base_altitude_flag .or. & + cld_base_0p1_okta_flag .or. & + cld_base_2p5_okta_flag .or. & cld_base_4p5_okta_flag .or. & cld_base_7p9_okta_flag .or. & - cld_top_altitude_flag .or. & + cld_top_2p5_okta_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & low_cld_amount_flag .or. & @@ -212,10 +216,11 @@ contains cld_amount_rnd_flag .or. & cld_amount_maxrnd_flag .or. & ceil_cld_amount_maxrnd_flag .or. & - cld_base_altitude_flag .or. & + cld_base_0p1_okta_flag .or. & + cld_base_2p5_okta_flag .or. & cld_base_4p5_okta_flag .or. & cld_base_7p9_okta_flag .or. & - cld_top_altitude_flag .or. & + cld_top_2p5_okta_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & low_cld_amount_flag .or. & @@ -254,10 +259,11 @@ contains cld_amount_rnd_flag .or. & cld_amount_maxrnd_flag .or. & ceil_cld_amount_maxrnd_flag .or. & - cld_base_altitude_flag .or. & + cld_base_0p1_okta_flag .or. & + cld_base_2p5_okta_flag .or. & cld_base_4p5_okta_flag .or. & cld_base_7p9_okta_flag .or. & - cld_top_altitude_flag .or. & + cld_top_2p5_okta_flag .or. & low_cld_base_altitude_flag .or. & very_low_cld_amount_flag .or. & low_cld_amount_flag .or. & @@ -279,10 +285,11 @@ contains cld_amount_rnd, & cld_amount_maxrnd, & ceil_cld_amount_maxrnd, & - cld_base_altitude, & + cld_base_0p1_okta, & + cld_base_2p5_okta, & cld_base_4p5_okta, & cld_base_7p9_okta, & - cld_top_altitude, & + cld_top_2p5_okta, & low_cld_base_altitude, & very_low_cld_amount, & low_cld_amount, & @@ -381,10 +388,11 @@ contains if (cld_amount_rnd_flag) call cld_amount_rnd%write_field() if (cld_amount_maxrnd_flag) call cld_amount_maxrnd%write_field() if (ceil_cld_amount_maxrnd_flag) call ceil_cld_amount_maxrnd%write_field() - if (cld_base_altitude_flag) call cld_base_altitude%write_field() + if (cld_base_0p1_okta_flag) call cld_base_0p1_okta%write_field() + if (cld_base_2p5_okta_flag) call cld_base_2p5_okta%write_field() if (cld_base_4p5_okta_flag) call cld_base_4p5_okta%write_field() if (cld_base_7p9_okta_flag) call cld_base_7p9_okta%write_field() - if (cld_top_altitude_flag) call cld_top_altitude%write_field() + if (cld_top_2p5_okta_flag) call cld_top_2p5_okta%write_field() if (low_cld_base_altitude_flag) call low_cld_base_altitude%write_field() if (very_low_cld_amount_flag) call very_low_cld_amount%write_field() if (low_cld_amount_flag) call low_cld_amount%write_field() diff --git a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 index d8353d332d..6d25ffc7b9 100644 --- a/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/cld_diags_kernel_mod.F90 @@ -26,16 +26,17 @@ module cld_diags_kernel_mod !> type, public, extends(kernel_type) :: cld_diags_kernel_type private - type(arg_type) :: meta_args(23) = (/ & + type(arg_type) :: meta_args(24) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! combined_cld_amount_wth arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_max arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_rnd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_amount_maxrnd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! ceil_cld_amount_maxrnd - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_altitude + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_0p1_okta + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_2p5_okta arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_4p5_okta arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_base_7p9_okta - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_top_altitude + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! cld_top_2p5_okta arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! low_cld_base_altitude arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! very_low_cld_amount arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! low_cld_amount @@ -69,10 +70,11 @@ module cld_diags_kernel_mod !> @param[in,out] cld_amount_rnd Cloud amount random overlap !> @param[in,out] cld_amount_maxrnd Cloud amount maximum-random overlap !> @param[in,out] ceil_cld_amount_maxrnd Ceilometer filtered cloud amount maximum-random overlap - !> @param[in,out] cld_base_altitude Cloud base altitude wrt sea level + !> @param[in,out] cld_base_0p1_okta Cloud base altitude for 0.1 okta cloud amount + !> @param[in,out] cld_base_2p5_okta Cloud base altitude wrt sea level !> @param[in,out] cld_base_4p5_okta Cloud base altitude for 4.5 okta cloud amount !> @param[in,out] cld_base_7p9_okta Cloud base altitude for 7.9 okta cloud amount - !> @param[in,out] cld_top_altitude Cloud top altitude wrt sea level + !> @param[in,out] cld_top_2p5_okta Cloud top altitude wrt sea level !> @param[in,out] low_cld_base_altitude Cloud base altitude wrt sea level for very low amount of cloud !> @param[in,out] very_low_cld_amount Maximum cloud amount below 111m !> @param[in,out] low_cld_amount Maximum cloud amount between 111 and 1949m above sea level @@ -107,10 +109,11 @@ subroutine cld_diags_code( nlayers, & cld_amount_rnd, & cld_amount_maxrnd, & ceil_cld_amount_maxrnd, & - cld_base_altitude, & + cld_base_0p1_okta, & + cld_base_2p5_okta, & cld_base_4p5_okta, & cld_base_7p9_okta, & - cld_top_altitude, & + cld_top_2p5_okta, & low_cld_base_altitude, & very_low_cld_amount, & low_cld_amount, & @@ -154,10 +157,11 @@ subroutine cld_diags_code( nlayers, & real(kind=r_def), pointer, intent(inout) :: cld_amount_rnd(:) real(kind=r_def), pointer, intent(inout) :: cld_amount_maxrnd(:) real(kind=r_def), pointer, intent(inout) :: ceil_cld_amount_maxrnd(:) - real(kind=r_def), pointer, intent(inout) :: cld_base_altitude(:) + real(kind=r_def), pointer, intent(inout) :: cld_base_0p1_okta(:) + real(kind=r_def), pointer, intent(inout) :: cld_base_2p5_okta(:) real(kind=r_def), pointer, intent(inout) :: cld_base_4p5_okta(:) real(kind=r_def), pointer, intent(inout) :: cld_base_7p9_okta(:) - real(kind=r_def), pointer, intent(inout) :: cld_top_altitude(:) + real(kind=r_def), pointer, intent(inout) :: cld_top_2p5_okta(:) real(kind=r_def), pointer, intent(inout) :: low_cld_base_altitude(:) real(kind=r_def), pointer, intent(inout) :: very_low_cld_amount(:) real(kind=r_def), pointer, intent(inout) :: low_cld_amount(:) @@ -195,7 +199,8 @@ subroutine cld_diags_code( nlayers, & ! ... 150 hPa: real(kind=r_def), parameter :: high_to_very_high = 13608.0_r_def ! metres ! When looking for cloud base or top: how much cloud cover defines cloud boundary. - real(kind=r_def), parameter :: cld_cover_for_cld_bdry = 2.5_r_def/8.0_r_def + real(kind=r_def), parameter :: cld_cover_for_0p1_okta = 0.1_r_def/8.0_r_def + real(kind=r_def), parameter :: cld_cover_for_2p5_okta = 2.5_r_def/8.0_r_def real(kind=r_def), parameter :: cld_cover_for_4p5_okta = 4.5_r_def/8.0_r_def real(kind=r_def), parameter :: cld_cover_for_7p9_okta = 7.9_r_def/8.0_r_def ! When looking for cloud base sometimes want much smaller amount of cloud. @@ -320,10 +325,11 @@ subroutine cld_diags_code( nlayers, & end if ! Find heights above sea level (asl) if required. - if (.not. associated(cld_base_altitude, empty_real_data) .or. & + if (.not. associated(cld_base_0p1_okta, empty_real_data) .or. & + .not. associated(cld_base_2p5_okta, empty_real_data) .or. & .not. associated(cld_base_4p5_okta, empty_real_data) .or. & .not. associated(cld_base_7p9_okta, empty_real_data) .or. & - .not. associated(cld_top_altitude, empty_real_data) .or. & + .not. associated(cld_top_2p5_okta, empty_real_data) .or. & .not. associated(low_cld_base_altitude, empty_real_data) .or. & .not. associated(very_low_cld_amount, empty_real_data) .or. & .not. associated(low_cld_amount, empty_real_data) .or. & @@ -341,16 +347,31 @@ subroutine cld_diags_code( nlayers, & end do end if - ! cld_base_altitude (in kilofeet) - if (.not. associated(cld_base_altitude, empty_real_data) ) then + ! cld_base_0p1_okta (in kilofeet) + if (.not. associated(cld_base_0p1_okta, empty_real_data) ) then ! As a default, set cloud-base to beyond top of model - cld_base_altitude(map_2d(1)) = 1.1_r_def * & + cld_base_0p1_okta(map_2d(1)) = 1.1_r_def * & z_asl_centre_of_levels(nlayers) * & m_to_kfeet do k = 1, nlayers - if ( combined_cld_amount(k) >= cld_cover_for_cld_bdry ) then - cld_base_altitude(map_2d(1)) = z_asl_base_of_levels(k) * m_to_kfeet + if ( combined_cld_amount(k) >= cld_cover_for_0p1_okta ) then + cld_base_0p1_okta(map_2d(1)) = z_asl_base_of_levels(k) * m_to_kfeet + exit + end if + end do + end if + + ! cld_base_2p5_okta (in kilofeet) + if (.not. associated(cld_base_2p5_okta, empty_real_data) ) then + ! As a default, set cloud-base to beyond top of model + cld_base_2p5_okta(map_2d(1)) = 1.1_r_def * & + z_asl_centre_of_levels(nlayers) * & + m_to_kfeet + + do k = 1, nlayers + if ( combined_cld_amount(k) >= cld_cover_for_2p5_okta ) then + cld_base_2p5_okta(map_2d(1)) = z_asl_base_of_levels(k) * m_to_kfeet exit end if end do @@ -386,14 +407,14 @@ subroutine cld_diags_code( nlayers, & end do end if - ! cld_top_altitude (in kilofeet) - if (.not. associated(cld_top_altitude, empty_real_data) ) then + ! cld_top_2p5_okta (in kilofeet) + if (.not. associated(cld_top_2p5_okta, empty_real_data) ) then ! As a default, set cloud-top to missing data - cld_top_altitude(map_2d(1)) = rmdi + cld_top_2p5_okta(map_2d(1)) = rmdi do k = nlayers-1, 1, -1 - if ( combined_cld_amount(k) >= cld_cover_for_cld_bdry ) then - cld_top_altitude(map_2d(1)) = z_asl_base_of_levels(k+1) * m_to_kfeet + if ( combined_cld_amount(k) >= cld_cover_for_2p5_okta ) then + cld_top_2p5_okta(map_2d(1)) = z_asl_base_of_levels(k+1) * m_to_kfeet exit end if end do diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_idealised.xml b/rose-stem/app/lfric_atm/file/file_def_diags_idealised.xml index 76e92d2188..2c03fa4f41 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_idealised.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_idealised.xml @@ -34,7 +34,7 @@ - + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_idealised1.xml b/rose-stem/app/lfric_atm/file/file_def_diags_idealised1.xml index 089f2a3007..5aa802ad2b 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_idealised1.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_idealised1.xml @@ -34,7 +34,7 @@ - + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_ls_and_jedi.xml b/rose-stem/app/lfric_atm/file/file_def_diags_ls_and_jedi.xml index e45ce213db..7e0feedcac 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_ls_and_jedi.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_ls_and_jedi.xml @@ -10,7 +10,7 @@ - + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_ral.xml b/rose-stem/app/lfric_atm/file/file_def_diags_ral.xml index baf6aa3ca0..b80952620b 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_ral.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_ral.xml @@ -19,7 +19,7 @@ - + diff --git a/rose-stem/app/lfric_atm/file/file_def_diags_ver.xml b/rose-stem/app/lfric_atm/file/file_def_diags_ver.xml index 5e28ad32eb..fde265167f 100644 --- a/rose-stem/app/lfric_atm/file/file_def_diags_ver.xml +++ b/rose-stem/app/lfric_atm/file/file_def_diags_ver.xml @@ -22,7 +22,7 @@ - + @@ -56,7 +56,7 @@ - + diff --git a/rose-stem/app/lfric_coupled_atmosphere/file/iodef_basic_gal.xml b/rose-stem/app/lfric_coupled_atmosphere/file/iodef_basic_gal.xml index beca0e13e8..dfd26a3d79 100644 --- a/rose-stem/app/lfric_coupled_atmosphere/file/iodef_basic_gal.xml +++ b/rose-stem/app/lfric_coupled_atmosphere/file/iodef_basic_gal.xml @@ -271,7 +271,7 @@ - +