diff --git a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 deleted file mode 100644 index eb7dc2d09..000000000 --- a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 +++ /dev/null @@ -1,1265 +0,0 @@ -!----------------------------------------------------------------------------- -! (c) Crown copyright 2020 Met Office. All rights reserved. -! The file LICENCE, distributed with this code, contains details of the terms -! under which the code may be used. -!----------------------------------------------------------------------------- - -!> @brief A two time-level iterative time-discretisation of the nonlinear -!! 3D equations. -module semi_implicit_timestep_alg_mod - - use constants_mod, only: i_def, r_def, l_def, str_def - use fs_continuity_mod, only: Wtheta, W2 - use log_mod, only: log_event, & - log_scratch_space, & - LOG_LEVEL_INFO - use extrusion_mod, only: TWOD - use namelist_mod, only: namelist_type - use sci_fem_constants_mod, only: get_mass_matrix_fe, & - get_mass_matrix_fv - use sci_field_bundle_builtins_mod, & - only: clone_bundle, & - bundle_axpy, & - add_bundle, & - copy_bundle, & - set_bundle_scalar - - ! Parent of this module's semi-implicit timestep type - use timestep_method_mod, only: timestep_method_type - - ! Configuration options - use section_choice_config_mod, only: cloud, cloud_um, & - aerosol, aerosol_um - use physics_config_mod, only: blayer_placement, & - blayer_placement_fast, & - convection_placement, & - convection_placement_fast, & - stochastic_physics_placement, & - stochastic_physics_placement_fast, & - smagorinsky_placement, & - smagorinsky_placement_outer - - use aerosol_config_mod, only: glomap_mode, & - glomap_mode_dust_and_clim, & - glomap_mode_ukca - - use formulation_config_mod, only: use_physics, dlayer_on, & - use_wavedynamics, & - moisture_formulation, & - moisture_formulation_dry, & - exner_from_eos - use io_config_mod, only: write_conservation_diag, write_diag, & - use_xios_io, diagnostic_frequency, & - checkpoint_read - use initialization_config_mod, only: init_option, & - init_option_checkpoint_dump, & - lbc_option_um2lfric_file - use mixed_solver_config_mod, only: guess_np1, reference_reset_time - use timestepping_config_mod, only: alpha, & - outer_iterations, inner_iterations, & - spinup_alpha - use transport_config_mod, only: cheap_update, & - transport_ageofair - use derived_config_mod, only: bundle_size - use boundaries_config_mod, only: limited_area, blend_frequency, & - blend_frequency_inner, & - blend_frequency_outer, & - blend_frequency_final - use finite_element_config_mod, only: element_order_h, & - element_order_v - - ! PSyKAl PSyclone kernels - use matrix_vector_kernel_mod, only: matrix_vector_kernel_type - use dg_inc_matrix_vector_kernel_mod, & - only: dg_inc_matrix_vector_kernel_type - use sci_enforce_bc_kernel_mod, only: enforce_bc_kernel_type - - ! Derived Types - use field_array_mod, only: field_array_type - use field_mod, only: field_type - use field_parent_mod, only: field_parent_type - use field_collection_mod, only: field_collection_type - use r_tran_field_mod, only: r_tran_field_type - use io_value_mod, only: io_value_type, & - get_io_value - use driver_modeldb_mod, only: modeldb_type - use mesh_mod, only: mesh_type - use mesh_collection_mod, only: mesh_collection - use model_clock_mod, only: model_clock_type - use operator_mod, only: operator_type - - ! Algorithms - use rhs_alg_mod, only: rhs_alg - use gungho_transport_control_alg_mod, & - only: gungho_transport_control_alg_init, & - gungho_transport_control_alg - - use si_operators_alg_mod, only: create_si_operators, & - compute_si_operators, & - final_si_operators - use fast_physics_alg_mod, only: fast_physics - use slow_physics_alg_mod, only: slow_physics - use checks_and_balances_alg_mod, & - only: check_fields - - use semi_implicit_solver_alg_mod, & - only: semi_implicit_solver_alg_init, & - semi_implicit_solver_alg_step, & - semi_implicit_solver_alg_final - use derive_exner_from_eos_alg_mod, & - only: derive_exner_from_eos - use sci_mass_matrix_solver_alg_mod, & - only: mass_matrix_solver_alg - use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg - use update_prognostic_scalars_alg_mod, & - only: update_prognostic_scalars_alg - use mixing_alg_mod, only: mixing_alg - use si_diagnostics_mod, only: output_diags_for_si - use predictors_alg_mod, only: predictors_alg - - ! LAM - use limited_area_lbc_alg_mod, only: lam_solver_lbc, & - lam_blend_lbc - use lam_rhs_alg_mod, only: calc_rhs_lbc, & - apply_mask_rhs - - ! Field mappings - use calc_phys_predictors_alg_mod, & - only: calc_phys_predictors_alg - use map_physics_fields_alg_mod, only: map_physics_fields_alg - - ! Moisture species - use mr_indices_mod, only: nummr - use moist_dyn_mod, only: num_moist_factors, gas_law - - ! Field indices - use field_indices_mod, only: igh_u, igh_t, igh_d, igh_p - - ! Mixing settings - use mixing_config_mod, only: smagorinsky - use smagorinsky_alg_mod, only: smagorinsky_alg - ! Physics routines called -#ifdef UM_PHYSICS - use cld_alg_mod, only: cld_alg - use aerosol_ukca_alg_mod, only: aerosol_ukca_alg - use casim_activate_alg_mod, only: casim_activate_alg -#endif - use cld_incs_mod, only: cld_incs_init, cld_incs_output - - use timing_mod, only: start_timing, stop_timing, tik, LPROF - use ageofair_alg_mod, only: ageofair_update - - implicit none - - private - - ! State object for the semi-implicit time-stepping method - type, extends(timestep_method_type), public :: semi_implicit_timestep_type - private - logical(kind=l_def) :: use_moisture - logical(kind=l_def) :: output_cld_incs - ! holds the latest estimate of the prognostic fields through the timestep - type(field_type), allocatable :: state(:) - ! prognostic fields at the start of timestep - type(field_type), allocatable :: state_n(:) - - ! prognostic fields after slow physics, i.e. state_n + slow incs - type( field_type ), allocatable :: state_after_slow(:) - type( field_type ), allocatable :: advected_state(:) - type( field_type ), allocatable :: mr_n(:), mr_after_adv(:), mr_after_slow(:) - type( field_type ), allocatable :: rhs_n(:), rhs_np1(:), rhs_adv(:) - type( field_type ), allocatable :: adv_inc_prev(:) - type( field_type ), allocatable :: rhs_phys(:), rhs_lbc(:) - type( field_type ) :: dtheta, dtheta_cld ! increment to theta - type( field_type ) :: du ! increment to u - type( field_type ) :: wind_prev ! u from previous iteration used for cheap transport update - type( r_tran_field_type ) :: total_dry_flux - - contains - private - - procedure, public :: step => semi_implicit_alg_step - procedure, public :: finalise => semi_implicit_alg_final - procedure, nopass :: run_init - procedure, nopass :: run_step - procedure, nopass :: conditional_collection_copy - - end type semi_implicit_timestep_type - - ! Constructor for type - interface semi_implicit_timestep_type - module procedure semi_implicit_alg_init - end interface semi_implicit_timestep_type - -contains - - !> Extracts data from modeldb to prepare for initialising object - !> @param[in] modeldb Holds the model state - function semi_implicit_alg_init(modeldb) result(self) - implicit none - - type(semi_implicit_timestep_type) :: self - - type(modeldb_type), intent(in), target :: modeldb - - type( field_collection_type ), pointer :: prognostic_fields - type( field_collection_type ), pointer :: moisture_fields - - type( field_type), pointer :: u - type( field_type), pointer :: rho - type( field_type), pointer :: theta - type( field_type), pointer :: exner - type( field_array_type ), pointer :: mr_array - type( field_type ), pointer :: mr(:) - class(model_clock_type), pointer :: model_clock - - ! Get pointer to clock to use downstream - model_clock => modeldb%clock - - ! Get pointers to field collections for use downstream - prognostic_fields => modeldb%fields%get_field_collection("prognostic_fields") - - ! Get pointers to fields in the prognostic/diagnostic field collections - ! for use downstream - call prognostic_fields%get_field('theta', theta) - call prognostic_fields%get_field('u', u) - call prognostic_fields%get_field('rho', rho) - call prognostic_fields%get_field('exner', exner) - - moisture_fields => modeldb%fields%get_field_collection("moisture_fields") - call moisture_fields%get_field("mr", mr_array) - mr => mr_array%bundle - - ! Run the initialisation process - call run_init(self, u, rho, theta, exner, mr, prognostic_fields, moisture_fields, model_clock ) - - nullify ( prognostic_fields, moisture_fields, mr_array, model_clock, & - u, rho, theta, exner, mr ) - - end function semi_implicit_alg_init - - !> Extracts data from input objects to prepare for running timestep - !> @param[in] modeldb The gungho model data object - subroutine semi_implicit_alg_step( self, modeldb ) - - implicit none - - class(semi_implicit_timestep_type), intent(inout) :: self - - type(modeldb_type), intent(in), target :: modeldb - - type(mesh_type), pointer :: mesh - type(mesh_type), pointer :: twod_mesh - - class(model_clock_type), pointer :: model_clock - - type( field_collection_type ), pointer :: prognostic_fields - type( field_collection_type ), pointer :: moisture_fields - type( field_type), pointer :: u - type( field_type), pointer :: rho - type( field_type), pointer :: theta - type( field_type), pointer :: exner - type( field_array_type ), pointer :: mr_array - type( field_type ), pointer :: mr(:) - type( field_array_type ), pointer :: moist_dyn_array - type( field_type ), pointer :: moist_dyn(:) - - type( field_collection_type ), pointer :: adv_tracer_all_outer - type( field_collection_type ), pointer :: adv_tracer_last_outer - type( field_collection_type ), pointer :: con_tracer_all_outer - type( field_collection_type ), pointer :: con_tracer_last_outer - type( field_collection_type ), pointer :: derived_fields - type( field_collection_type ), pointer :: radiation_fields - type( field_collection_type ), pointer :: microphysics_fields - type( field_collection_type ), pointer :: electric_fields - type( field_collection_type ), pointer :: orography_fields - type( field_collection_type ), pointer :: turbulence_fields - type( field_collection_type ), pointer :: convection_fields - type( field_collection_type ), pointer :: cloud_fields - type( field_collection_type ), pointer :: surface_fields - type( field_collection_type ), pointer :: soil_fields - type( field_collection_type ), pointer :: snow_fields - type( field_collection_type ), pointer :: chemistry_fields - type( field_collection_type ), pointer :: aerosol_fields - type( field_collection_type ), pointer :: stph_fields - type( field_collection_type ), pointer :: lbc_fields - - type( io_value_type ), pointer :: temp_corr_io_value - real( r_def ) :: dt - real( r_def ) :: dtemp_encorr - - ! Get pointers to field collections for use downstream - prognostic_fields => modeldb%fields%get_field_collection("prognostic_fields") - model_clock => modeldb%clock - - ! Get pointers to fields for use downstream - call prognostic_fields%get_field('theta', theta) - call prognostic_fields%get_field('u', u) - call prognostic_fields%get_field('rho', rho) - call prognostic_fields%get_field('exner', exner) - - ! Get timestep parameters from clock - dt = real(model_clock%get_seconds_per_step(), r_def) - - moisture_fields => modeldb%fields%get_field_collection("moisture_fields") - lbc_fields => modeldb%fields%get_field_collection("lbc_fields") - radiation_fields => modeldb%fields%get_field_collection("radiation_fields") - call moisture_fields%get_field("mr", mr_array) - mr => mr_array%bundle - call moisture_fields%get_field("moist_dyn", moist_dyn_array) - moist_dyn => moist_dyn_array%bundle - - ! Get the mesh from one of the fields - mesh => theta%get_mesh() - twod_mesh => mesh_collection%get_mesh(mesh, TWOD) - - adv_tracer_all_outer => modeldb%fields%get_field_collection("adv_tracer_all_outer") - adv_tracer_last_outer => modeldb%fields%get_field_collection("adv_tracer_last_outer") - con_tracer_all_outer => modeldb%fields%get_field_collection("con_tracer_all_outer") - con_tracer_last_outer => modeldb%fields%get_field_collection("con_tracer_last_outer") - derived_fields => modeldb%fields%get_field_collection("derived_fields") - microphysics_fields => modeldb%fields%get_field_collection("microphysics_fields") - turbulence_fields => modeldb%fields%get_field_collection("turbulence_fields") - convection_fields => modeldb%fields%get_field_collection("convection_fields") - cloud_fields => modeldb%fields%get_field_collection("cloud_fields") - surface_fields => modeldb%fields%get_field_collection("surface_fields") - soil_fields => modeldb%fields%get_field_collection("soil_fields") - snow_fields => modeldb%fields%get_field_collection("snow_fields") - chemistry_fields => modeldb%fields%get_field_collection("chemistry_fields") - aerosol_fields => modeldb%fields%get_field_collection("aerosol_fields") - stph_fields => modeldb%fields%get_field_collection("stph_fields") - electric_fields => modeldb%fields%get_field_collection("electric_fields") - orography_fields => modeldb%fields%get_field_collection("orography_fields") - - temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') - ! Get temperature increment for energy correction - dtemp_encorr = dt * temp_corr_io_value%data(1) - - ! Run the timestep - call run_step(self, modeldb, & - u, rho, theta, exner, mr, moist_dyn, & - adv_tracer_all_outer,adv_tracer_last_outer, & - con_tracer_all_outer,con_tracer_last_outer, & - prognostic_fields, moisture_fields, & - derived_fields, radiation_fields, & - microphysics_fields, electric_fields, & - orography_fields, & - turbulence_fields, convection_fields, & - cloud_fields, surface_fields, soil_fields, & - snow_fields, chemistry_fields, & - aerosol_fields, stph_fields, lbc_fields, & - model_clock, dtemp_encorr, & - mesh, twod_mesh) - - nullify( model_clock, mesh, & - prognostic_fields, moisture_fields, & - u, rho, theta, exner, mr, moist_dyn, & - adv_tracer_all_outer,adv_tracer_last_outer, & - con_tracer_all_outer,con_tracer_last_outer, & - derived_fields, radiation_fields, & - microphysics_fields, electric_fields, & - orography_fields, & - turbulence_fields, convection_fields, & - cloud_fields, surface_fields, soil_fields, & - snow_fields, chemistry_fields, & - aerosol_fields, stph_fields, lbc_fields, & - temp_corr_io_value, & - mr_array, moist_dyn_array) - - end subroutine semi_implicit_alg_step - - !> @details Initialisation procedure for the timestepping algorithm - !> Initialises various internal fields - !> @param[in,out] u 3D wind field - !> @param[in,out] rho Density - !> @param[in,out] theta Potential temperature - !> @param[in,out] exner Exner pressure - !> @param[in,out] mr Mixing ratios - !> @param[in,out] prognostic_fields Prognostic field collection - !> @param[in,out] moisture_fields Moisture field collection - !> @param[in] model_clock Time in the model. - subroutine run_init( self, u, rho, theta, exner, mr, prognostic_fields, moisture_fields, model_clock) - - implicit none - - type(semi_implicit_timestep_type), intent(inout) :: self - - ! Prognostic fields - type( field_type ), intent( in ) :: u, rho, theta, exner - type( field_type ), intent( in ), optional :: mr(:) - - ! field groups - type( field_collection_type ), intent( inout ) :: prognostic_fields - type( field_collection_type ), intent( inout ) :: moisture_fields - - ! Clock - class(model_clock_type), intent(in) :: model_clock - - ! Mesh - type(mesh_type), pointer :: mesh - - ! Reference fields - type( field_type ), pointer :: rho_ref, theta_ref, exner_ref - - ! Reference moist dynamics factors - type( field_array_type ), pointer :: moist_dyn_ref_array - type( field_type ), pointer :: moist_dyn_ref(:) - - - ! Timestep - real(kind=r_def) :: cast_dt - - ! Reset frequency of semi-implicit operators - integer(kind=i_def) :: reference_reset_freq - - ! Get the mesh from one of the fields - mesh => theta%get_mesh() - - !-------------------------------------------------------------------- - ! Allocate internal state field arrays - !-------------------------------------------------------------------- - allocate(self%state(bundle_size)) - allocate(self%state_n(bundle_size)) - allocate(self%state_after_slow(bundle_size)) - allocate(self%advected_state(bundle_size)) - allocate(self%rhs_n(bundle_size)) - allocate(self%rhs_np1(bundle_size)) - allocate(self%rhs_adv(bundle_size)) - allocate(self%adv_inc_prev(bundle_size)) - allocate(self%rhs_phys(bundle_size)) - allocate(self%rhs_lbc(bundle_size)) - allocate(self%mr_n(nummr)) - allocate(self%mr_after_slow(nummr)) - allocate(self%mr_after_adv(nummr)) - - self%use_moisture = ( moisture_formulation /= moisture_formulation_dry ) - self%output_cld_incs = ( self%use_moisture .and. write_diag .and. & - use_xios_io ) - - !-------------------------------------------------------------------- - ! Initialise internal state field objects - !-------------------------------------------------------------------- - - call u%copy_field_properties( self%state(igh_u) ) - call theta%copy_field_properties( self%state(igh_t) ) - call rho%copy_field_properties( self%state(igh_d) ) - call exner%copy_field_properties( self%state(igh_p) ) - - call clone_bundle(self%state, self%state_n, bundle_size) - call clone_bundle(self%state, self%state_after_slow, bundle_size) - call clone_bundle(self%state, self%advected_state, bundle_size) - call clone_bundle(self%state, self%rhs_n, bundle_size) - call clone_bundle(self%state, self%rhs_np1, bundle_size) - call clone_bundle(self%state, self%rhs_adv, bundle_size) - call clone_bundle(self%state, self%rhs_phys, bundle_size) - call clone_bundle(self%state, self%rhs_lbc, bundle_size) - call clone_bundle(self%state, self%adv_inc_prev, bundle_size) - - call theta%copy_field_properties(self%dtheta) - call theta%copy_field_properties(self%dtheta_cld) - call invoke(setval_c(self%dtheta_cld, 0.0_r_def)) - call u%copy_field_properties(self%du) - call u%copy_field_properties(self%wind_prev) - - call self%total_dry_flux%initialise( u%get_function_space() ) - - call clone_bundle(mr, self%mr_n, nummr) - call clone_bundle(mr, self%mr_after_slow, nummr) - if (self%use_moisture) then - call clone_bundle(mr, self%mr_after_adv, nummr) - else - call set_bundle_scalar(0.0_r_def, self%mr_n, nummr) - call set_bundle_scalar(0.0_r_def, self%mr_after_slow, nummr) - end if - - !-------------------------------------------------------------------- - ! Initialise the physics increments to 0 - !-------------------------------------------------------------------- - call set_bundle_scalar(0.0_r_def, self%rhs_phys, bundle_size) - - !-------------------------------------------------------------------- - ! Operators for si solves - !-------------------------------------------------------------------- - call create_si_operators( mesh ) - - ! If using checkpointed reference state, then calculate semi-implicit - ! operators using the checkpointed reference state - cast_dt = real(model_clock%get_seconds_per_step(), r_def) - reference_reset_freq = nint(reference_reset_time / cast_dt, i_def) - if (mod(model_clock%get_first_step()-1, reference_reset_freq) /= 0) then - call moisture_fields%get_field("moist_dyn_ref", moist_dyn_ref_array) - moist_dyn_ref => moist_dyn_ref_array%bundle - call prognostic_fields%get_field('theta_ref', theta_ref) - call prognostic_fields%get_field('rho_ref', rho_ref) - call prognostic_fields%get_field('exner_ref', exner_ref) - call compute_si_operators(theta_ref, rho_ref, exner_ref, model_clock, moist_dyn_ref) - nullify(theta_ref, rho_ref, exner_ref, moist_dyn_ref, moist_dyn_ref_array) - end if - - if (use_wavedynamics) then - call gungho_transport_control_alg_init( mesh ) - - ! Construct semi-implicit solver - call semi_implicit_solver_alg_init( self%state ) - end if - - nullify(mesh) - - call log_event( "semi_implicit_timestep: initialised timestepping algorithm", LOG_LEVEL_INFO ) - - end subroutine run_init - - !> @details An algorithm for timestepping the 3D nonlinear equations - !> using an iterative process of the same form as endgame. - !> The algorithm splits all processes into one of three parts. - !> Old time level forcings computed once per timestep. - !> Advection terms computed in an outer loop using time-averaged - !> fields. - !> New time level forcings computed in an inner loop. - !> If matching ENDGame, 2 outer and 2 inner loops are used per timestep - !> by default. This means that there is one evaluation of old time - !> level terms, 2 evaluation of advective terms and 4 evaluations - !> of new time level terms and increment updates per timestep - !> @param[in] modeldb Holds the model state - !> @param[in,out] u 3D wind field - !> @param[in,out] rho Density - !> @param[in,out] theta Potential temperature - !> @param[in,out] exner Exner pressure - !> @param[in,out] mr Mixing ratios - !> @param[in,out] moist_dyn Factors for moist dynamics - !> @param[in,out] adv_tracer_all_outer Group of fields to be advected - !> @param[in,out] adv_tracer_last_outer Group of fields to be advected - !> @param[in,out] con_tracer_all_outer Second group of fields to be advected - !> @param[in,out] con_tracer_last_outer Second group of fields to be advected - !> @param[in,out] prognostic_fields Prognostic field collection - !> @param[in,out] moisture_fields Moisture field collection - !> @param[in,out] derived_fields Group of derived fields - !> @param[in,out] radiation_fields Fields for radiation scheme - !> @param[in,out] microphysics_fields Fields for mphys scheme - !> @param[in,out] electric_fields Fields for electric (lighting) scheme - !> @param[in] orography_fields Fields for orog drag scheme - !> @param[in,out] turbulence_fields Fields for turbulence scheme - !> @param[in,out] convection_fields Fields for convection scheme - !> @param[in,out] cloud_fields Fields for cloud scheme - !> @param[in,out] surface_fields Fields for surface scheme - !> @param[in,out] soil_fields Fields for soil hydrology scheme - !> @param[in,out] snow_fields Fields for snow scheme - !> @param[in,out] chemistry_fields Fields for chemistry scheme - !> @param[in,out] aerosol_fields Fields for aerosol scheme - !> @param[in,out] stph_fields Fields for stohcastic physics schemes - !> @param[in,out] lbc_fields Fields for lateral boundaries - !> @param[in] model_clock Time in the model. - !> @param[in] dtemp_encorr Temperature increment for energy - !> correction - !> @param[in] mesh The current mesh - !> @param[in] twod_mesh The current 2d mesh - subroutine run_step(self, modeldb, & - u, rho, theta, exner, mr, moist_dyn, & - adv_tracer_all_outer,adv_tracer_last_outer, & - con_tracer_all_outer,con_tracer_last_outer, & - prognostic_fields, moisture_fields, & - derived_fields, radiation_fields, & - microphysics_fields, electric_fields, & - orography_fields, & - turbulence_fields, convection_fields, & - cloud_fields, surface_fields, soil_fields, & - snow_fields, chemistry_fields, & - aerosol_fields, stph_fields, lbc_fields, & - model_clock, dtemp_encorr, & - mesh, twod_mesh) - - implicit none - - type(semi_implicit_timestep_type), intent(inout), target :: self - type(modeldb_type), intent(in), target :: modeldb - - - ! Prognostic fields - type( field_type ), intent( inout ) :: u, rho, theta, exner - type( field_type ), intent( inout ) :: mr(nummr) - type( field_type ), intent( inout ) :: moist_dyn(num_moist_factors) - ! field groups - type( field_collection_type ), intent( inout ) :: adv_tracer_all_outer - type( field_collection_type ), intent( inout ) :: adv_tracer_last_outer - type( field_collection_type ), intent( inout ) :: con_tracer_all_outer - type( field_collection_type ), intent( inout ) :: con_tracer_last_outer - type( field_collection_type ), intent( inout ) :: prognostic_fields - type( field_collection_type ), intent( inout ) :: moisture_fields - type( field_collection_type ), intent( inout ) :: derived_fields - type( field_collection_type ), intent( inout ) :: radiation_fields - type( field_collection_type ), intent( inout ) :: microphysics_fields - type( field_collection_type ), intent( inout ) :: electric_fields - type( field_collection_type ), intent( in ) :: orography_fields - type( field_collection_type ), intent( inout ) :: turbulence_fields - type( field_collection_type ), intent( inout ) :: convection_fields - type( field_collection_type ), intent( inout ) :: cloud_fields - type( field_collection_type ), intent( inout ) :: surface_fields - type( field_collection_type ), intent( inout ) :: soil_fields - type( field_collection_type ), intent( inout ) :: snow_fields - type( field_collection_type ), intent( inout ) :: chemistry_fields - type( field_collection_type ), intent( inout ) :: aerosol_fields - type( field_collection_type ), intent( inout ) :: stph_fields - type( field_collection_type ), intent( inout ) :: lbc_fields - - class(model_clock_type), intent(in) :: model_clock - - real(kind=r_def), intent(in) :: dtemp_encorr - - type(mesh_type), intent(in), pointer :: mesh - type(mesh_type), intent(in), pointer :: twod_mesh - - ! Reference fields - type( field_type ), pointer :: rho_ref - type( field_type ), pointer :: theta_ref - type( field_type ), pointer :: exner_ref - type( field_type ), pointer :: moist_dyn_ref(:) - type( field_array_type ), pointer :: moist_dyn_ref_array - - ! Moisture field to transport - type( field_type ), pointer :: mr_to_adv(:) - - type( field_type ) :: dcfl_tot, dcff_tot, dbcf_tot, & - dcfl_adv, dcff_adv, dbcf_adv - character(str_def), parameter :: sec_tot='processed' - character(str_def), parameter :: suffix_tot='tot' - character(str_def), parameter :: sec_adv='advection' - character(str_def), parameter :: suffix_adv='adv' - - real(kind=r_def) :: cast_dt - - type(field_type), pointer :: ageofair - - type(operator_type), pointer :: mm_wt - type(operator_type), pointer :: mm_vel - - integer(kind=i_def) :: outer, inner, reference_reset_freq - real(kind=r_def) :: varalpha, varbeta ! alpha, beta weight to use - ! these may differ from input values - ! during the spinup period - - ! Density field used for moisture conservation diagnostics and predictor - logical(kind=l_def) :: write_moisture_diag - - ! Fields after slow physics to be advected (i.e. field_n + slow phys inc) - type( field_collection_type ) :: adv_tracer_all_outer_after_slow - type( field_collection_type ) :: adv_tracer_last_outer_after_slow - type( field_collection_type ) :: con_tracer_all_outer_after_slow - type( field_collection_type ) :: con_tracer_last_outer_after_slow - - ! Reference fields are checkpointed - logical(kind=l_def) :: checkpoint_reference_fields - - ! Namelist parameters - type(namelist_type), pointer :: base_mesh_nml - type(namelist_type), pointer :: initialization_nml - type(namelist_type), pointer :: microphysics_nml - type(namelist_type), pointer :: aerosol_nml - type(namelist_type), pointer :: timestepping_nml - - character(str_def) :: prime_mesh_name - integer(i_def) :: lbc_option - logical(l_def) :: microphysics_casim - logical(l_def) :: murk_lbc - real(r_def) :: tau_r - integer(tik) :: id - - if ( LPROF ) call start_timing( id, 'semi_implicit_timestep' ) - - cast_dt = real(model_clock%get_seconds_per_step(), r_def) - - if (limited_area .and. use_wavedynamics) then - base_mesh_nml => modeldb%configuration%get_namelist('base_mesh') - initialization_nml => modeldb%configuration%get_namelist('initialization') - timestepping_nml => modeldb%configuration%get_namelist('timestepping') - - call base_mesh_nml%get_value( 'prime_mesh_name', prime_mesh_name ) - call initialization_nml%get_value( 'lbc_option', lbc_option ) - call timestepping_nml%get_value( 'tau_r', tau_r ) - - if (lbc_option == lbc_option_um2lfric_file) then - aerosol_nml => modeldb%configuration%get_namelist('aerosol') - call aerosol_nml%get_value( 'murk_lbc', murk_lbc ) - end if - end if - if (lbc_option == lbc_option_um2lfric_file .or. & - (use_physics .and. cloud == cloud_um)) then - microphysics_nml => modeldb%configuration%get_namelist('microphysics') - call microphysics_nml%get_value( 'microphysics_casim', microphysics_casim ) - end if - - if (element_order_h == 0 .and. element_order_v == 0) then - ! Lowest order so use finite volume constants - mm_wt => get_mass_matrix_fv(Wtheta, mesh%get_id()) - mm_vel => get_mass_matrix_fv(W2, mesh%get_id()) - else - mm_wt => get_mass_matrix_fe(Wtheta, mesh%get_id()) - mm_vel => get_mass_matrix_fe(W2, mesh%get_id()) - end if - - !-------------------------------------------------------------------- - ! Copy prognostic field data to state arrays - !-------------------------------------------------------------------- - call invoke( name = "copy_init_fields_to_state", & - setval_X(self%state(igh_u), u ), & - setval_X(self%state(igh_t), theta), & - setval_X(self%state(igh_d), rho ), & - setval_X(self%state(igh_p), exner), & - setval_c(self%total_dry_flux, 0.0_r_tran) ) - - !-------------------------------------------------------------------- - ! Compute the semi-implicit operators - !-------------------------------------------------------------------- - ! Reset the reference state in the semi-implicit operators using the latest state guess. - ! This occurs every n timesteps, where n is calculated as reference_reset_time divided by dt. - ! The reference_reset_time is specified in the configuration namelist. - ! Note that this reset can only occur at most once per timestep. - reference_reset_freq = nint(reference_reset_time / cast_dt, i_def) - - if ( mod(model_clock%get_step() - 1_i_def, reference_reset_freq) == 0_i_def ) then - ! Compute semi-implicit operators with current model state - call compute_si_operators(self%state(igh_t), self%state(igh_d), self%state(igh_p), model_clock, moist_dyn) - - checkpoint_reference_fields = & - mod(model_clock%get_first_step()-1, reference_reset_freq) /= 0 .or. & - mod(model_clock%get_last_step(), reference_reset_freq) /= 0 - if (checkpoint_read .or. init_option == init_option_checkpoint_dump) then - ! If the first timestep of this run IS an operator calc timestep, but - ! the first timestep of the next run IS NOT, then checkpoint_flag - ! must be false to allow model to start running, as the operator - ! prognostics will not be in the initial dump - if (mod(model_clock%get_first_step()-1, reference_reset_freq) == 0 .and. & - mod(model_clock%get_last_step(), reference_reset_freq) /= 0) then - checkpoint_reference_fields = .false. - end if - end if - ! Copy over model state for checkpointing if required - if (checkpoint_reference_fields) then - call prognostic_fields%get_field('theta_ref', theta_ref) - call prognostic_fields%get_field('rho_ref', rho_ref) - call prognostic_fields%get_field('exner_ref', exner_ref) - call moisture_fields%get_field("moist_dyn_ref", moist_dyn_ref_array) - moist_dyn_ref => moist_dyn_ref_array%bundle - call copy_bundle(moist_dyn, moist_dyn_ref, num_moist_factors) - call invoke( setval_X(theta_ref, self%state(igh_t)), & - setval_X(rho_ref, self%state(igh_d)), & - setval_X(exner_ref, self%state(igh_p)) ) - nullify(moist_dyn_ref_array, moist_dyn_ref, theta_ref, rho_ref, exner_ref) - end if - end if - - !-------------------------------------------------------------------- - ! If off-centering is being spun up then modify the alpha value - !-------------------------------------------------------------------- - if (spinup_alpha .and. model_clock%is_spinning_up()) then - varalpha = 1.0_r_def - else - varalpha = alpha - end if - varbeta = 1.0_r_def - varalpha - - ! Perform some checking on the fields. - call check_fields(self%state, cast_dt) - - !-------------------------------------------------------------------- - ! Update state_n and mr_n with start of timestep values - !-------------------------------------------------------------------- - if (self%use_moisture) then - call copy_bundle(mr, self%mr_n, nummr) - call copy_bundle(mr, self%mr_after_slow, nummr) - end if - call copy_bundle(self%state, self%state_n, bundle_size) - call copy_bundle(self%state, self%state_after_slow, bundle_size) - - !-------------------------------------------------------------------- - ! Compute slow physic updates - !-------------------------------------------------------------------- - if (use_physics) then - if (self%output_cld_incs) then - call cld_incs_init(cloud_fields, dcfl_tot, dcff_tot, dbcf_tot, & - sec_tot, suffix_tot) - end if - - call slow_physics( modeldb, self%du, self%dtheta, self%mr_after_slow, & - self%state_n(igh_t), self%state_n(igh_u), & - self%state_n(igh_d), self%state_n(igh_p), moist_dyn, & - self%mr_n, derived_fields, radiation_fields, & - microphysics_fields, electric_fields, & - orography_fields, & - turbulence_fields, convection_fields, cloud_fields, & - surface_fields, soil_fields, snow_fields, & - chemistry_fields, aerosol_fields, model_clock, & - cast_dt, dtemp_encorr, mesh, twod_mesh ) - - call invoke(name="update_from_slow_physics", & - inc_X_plus_Y(self%state_after_slow(igh_t), self%dtheta), & - inc_X_plus_Y(self%state_after_slow(igh_u), self%du) & - ) - - if (self%output_cld_incs) then - call cld_incs_init(cloud_fields, dcfl_adv, dcff_adv, dbcf_adv, & - sec_adv, suffix_adv) - end if - end if !use_physics - - !-------------------------------------------------------------------- - ! Compute the time-level n dynamics terms - !-------------------------------------------------------------------- - call rhs_alg( self%rhs_n, varbeta*cast_dt, & - self%state_after_slow, self%state_n, moist_dyn, & - compute_eos=.false., compute_rhs_t_d=.true., & - dlayer_rhs=.false., model_clock=model_clock ) - - call copy_bundle(self%state_after_slow, self%advected_state, bundle_size) - ! Set the moisture to be transported to point to the moisture after slow physics - mr_to_adv => self%mr_after_slow - ! Set predictors for transport - call predictors_alg(self%advected_state, self%state_n(igh_u), & - self%rhs_n(igh_u),varbeta,model_clock) - - !========================================================================== - ! Start the Outer (advection) loop - !========================================================================== - call conditional_collection_copy(adv_tracer_all_outer_after_slow, & - generic_fields_to_copy=adv_tracer_all_outer, & - field_list=adv_tracer_all_outer) - call conditional_collection_copy(adv_tracer_last_outer_after_slow, & - generic_fields_to_copy=adv_tracer_last_outer, & - field_list=adv_tracer_last_outer) - call conditional_collection_copy(con_tracer_all_outer_after_slow, & - generic_fields_to_copy=con_tracer_all_outer, & - field_list=con_tracer_all_outer) - call conditional_collection_copy(con_tracer_last_outer_after_slow, & - generic_fields_to_copy=con_tracer_last_outer, & - field_list=con_tracer_last_outer) - call invoke( setval_X(self%wind_prev, self%state(igh_u)) ) - - outer_dynamics_loop: do outer = 1,outer_iterations - - if (use_wavedynamics) then - call gungho_transport_control_alg( & - self%rhs_adv, self%advected_state, self%state(igh_u), & - self%state_n(igh_u), mr, mr_to_adv, model_clock, & - outer, cheap_update, self%adv_inc_prev, self%wind_prev, & - self%state_after_slow(igh_d), self%total_dry_flux, & - adv_tracer_all_outer, adv_tracer_all_outer_after_slow, & - adv_tracer_last_outer, adv_tracer_last_outer_after_slow, & - con_tracer_all_outer, con_tracer_all_outer_after_slow, & - con_tracer_last_outer, con_tracer_last_outer_after_slow & - ) - - if ( cheap_update .AND. (outer < outer_iterations) ) then - ! Copy transport increments to use in next outer iteration - ! for the cheap transport update - call copy_bundle(self%rhs_adv, self%adv_inc_prev, bundle_size) - ! Update advected state, moisture, and tracers to be transported - ! in the next outer iteration transport call - call mass_matrix_solver_alg( self%du, self%rhs_adv(igh_u) ) - call invoke( inc_X_plus_Y(self%advected_state(igh_d), self%rhs_adv(igh_d) ), & - inc_X_plus_Y(self%advected_state(igh_t), self%rhs_adv(igh_t) ), & - inc_X_plus_Y(self%advected_state(igh_u), self%du ) ) - call adv_tracer_all_outer_after_slow%clear() - call con_tracer_all_outer_after_slow%clear() - call conditional_collection_copy(adv_tracer_all_outer_after_slow, & - generic_fields_to_copy=adv_tracer_all_outer, & - field_list=adv_tracer_all_outer) - call conditional_collection_copy(con_tracer_all_outer_after_slow, & - generic_fields_to_copy=con_tracer_all_outer, & - field_list=con_tracer_all_outer) - if (self%use_moisture) then - ! Update the moisture to be transported in the next outer iteration - ! by setting mr_to_adv to point to the moisture after the previous - ! transport, mr_after_adv (which is set later in the outer loop) - mr_to_adv => self%mr_after_adv - end if - ! Store latest estimate of the wind for cheap transport update - call invoke( setval_X(self%wind_prev, self%state(igh_u)) ) - end if - - ! Convert theta increment to weak form - call invoke( setval_X(self%dtheta, self%rhs_adv(igh_t)), & - setval_c(self%rhs_adv(igh_t), 0.0_r_def), & - dg_inc_matrix_vector_kernel_type(self%rhs_adv(igh_t), & - self%dtheta, mm_wt) ) - - ! Compute the time-level n+1 dynamics terms - call rhs_alg( self%rhs_np1, -varalpha*cast_dt, & - self%state, self%state, moist_dyn, & - compute_eos=.true., compute_rhs_t_d=.true., & - dlayer_rhs=dlayer_on, model_clock=model_clock ) - else - if (self%use_moisture) call copy_bundle(self%mr_after_slow, mr, nummr) - end if - - ! Update the moisture after transport (mr_after_adv) - if (self%use_moisture) call copy_bundle(mr, self%mr_after_adv, nummr) - - if (use_physics) then - !-------------------------------------------------------------------- - ! setting predictors for fast physics - !-------------------------------------------------------------------- - if (blayer_placement == blayer_placement_fast .or. & - convection_placement == convection_placement_fast .or. & - stochastic_physics_placement == stochastic_physics_placement_fast)& - then - call calc_phys_predictors_alg( derived_fields, self%rhs_np1, self%rhs_adv, & - self%rhs_n, self%state, self%state_after_slow, & - lbc_fields, model_clock ) - end if - - if (outer == outer_iterations .and. self%output_cld_incs) then - call cld_incs_output(cloud_fields, dcfl_adv, dcff_adv, dbcf_adv, & - sec_adv, suffix_adv) - end if - !-------------------------------------------------------------------- - ! Call the fast physics terms - !-------------------------------------------------------------------- - call fast_physics(self%du, self%dtheta, mr, & - self%state_n(igh_t), self%state_n(igh_d), & - self%state_n(igh_u), self%state_n(igh_p), & - self%mr_n, derived_fields, & - radiation_fields, microphysics_fields, & - orography_fields, turbulence_fields, & - convection_fields, cloud_fields, & - surface_fields, soil_fields, snow_fields, & - chemistry_fields, aerosol_fields, stph_fields, & - outer, model_clock, cast_dt) - - if (smagorinsky .and. smagorinsky_placement == smagorinsky_placement_outer) then - call smagorinsky_alg(self%dtheta, self%du, mr, self%state(igh_t), & - self%state(igh_u), & - derived_fields, self%state(igh_d), & - cast_dt) - end if - - if (use_wavedynamics) then - ! copy increments into rhs_phys, including premultiplication by mass matrix - ! need to reset rhs_phys to 0 because matrix_vector_kernel_type - ! increments the field rather than over-writing it - call set_bundle_scalar(0.0_r_def, self%rhs_phys, bundle_size) - call invoke(name="update_rhs_phys_from_fast_physics", & - dg_inc_matrix_vector_kernel_type(self%rhs_phys(igh_t), self%dtheta, & - mm_wt), & - matrix_vector_kernel_type(self%rhs_phys(igh_u), self%du, mm_vel ), & - enforce_bc_kernel_type(self%rhs_phys(igh_u)) & - ) - end if - - end if !use_physics - - if (use_wavedynamics) then - - ! Use advective update to guess n+1 level scalar fields. - if ( guess_np1 ) then - ! Update factors for moist dynamics - if (self%use_moisture) call moist_dyn_factors_alg(moist_dyn, mr) - call update_prognostic_scalars_alg(self%state, self%rhs_n, self%rhs_adv, & - self%rhs_phys, & - moist_dyn(gas_law)) - end if - - !============================================================================ - ! Start the Inner (nonlinear, coriolis) loop - !============================================================================ - inner_dynamics_loop: do inner = 1,inner_iterations - write( log_scratch_space, '(A,2I3)' ) 'loop indices (o, i): ', & - outer, inner - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - !-------------------------------------------------------------------- - ! Compute the time-level n+1 dynamics terms - !-------------------------------------------------------------------- - if (inner > 1) call rhs_alg( self%rhs_np1, -varalpha*cast_dt, & - self%state, self%state, moist_dyn, & - compute_eos=.true., compute_rhs_t_d=.false., & - dlayer_rhs=dlayer_on, model_clock=model_clock ) - - !-------------------------------------------------------------------- - ! Compute the LAM LBCs and RHS - !-------------------------------------------------------------------- - if ( limited_area .and. inner == 1 .and. outer == 1 ) then - call lam_solver_lbc(self%state(igh_u), lbc_fields, prime_mesh_name) - call calc_rhs_lbc(self%rhs_lbc, lbc_fields, model_clock, prime_mesh_name, & - tau_r) - end if - - !-------------------------------------------------------------------- - ! Compute the residuals - ! - ! Add on advective terms: rhs = rhs_n - rhs_np1 + rhs_adv + rhs_phys - ! (reuse rhs_np1 for rhs) - !-------------------------------------------------------------------- - call bundle_axpy(-1.0_r_def, self%rhs_np1, self%rhs_n, self%rhs_np1, bundle_size) - call add_bundle(self%rhs_np1, self%rhs_adv, self%rhs_np1, bundle_size) - call add_bundle(self%rhs_np1, self%rhs_phys, self%rhs_np1, bundle_size) - - if ( limited_area ) then - if ( inner == 1 .and. outer == 1 ) then - ! Add on the RHS for LBCs - call add_bundle(self%rhs_np1, self%rhs_lbc, self%rhs_np1, bundle_size) - end if - ! Apply masks to RHS - call apply_mask_rhs(self%rhs_np1, prime_mesh_name) - end if - - ! Accelerators for inner loop convergence - if ( inner > 1 ) then - call invoke( setval_c(self%rhs_np1(igh_d), 0.0_r_def), & - setval_c(self%rhs_np1(igh_t), 0.0_r_def) ) - end if - - - write_moisture_diag = write_conservation_diag .and. & - outer == outer_iterations .and. & - inner == inner_iterations .and. & - self%use_moisture - !-------------------------------------------------------------------- - ! Solve semi-implicit system: A*inc = rhs, and incement state by inc - !-------------------------------------------------------------------- - call semi_implicit_solver_alg_step( self%state, self%rhs_np1, & - moist_dyn(gas_law), & - mr, & - write_moisture_diag, & - first_iteration=(inner==1) ) - - ! If not already done update factors for moist dynamics - if ( .not. guess_np1 .and. self%use_moisture ) & - call moist_dyn_factors_alg(moist_dyn, mr) - - if (exner_from_eos) then - call derive_exner_from_eos( self%state, & - moist_dyn(gas_law) ) - end if - - !-------------------------------------------------------------------- - ! LAM Overwrite and Blend LBCs - !-------------------------------------------------------------------- - if ( limited_area ) then - if ( ( blend_frequency == blend_frequency_inner) & - .or. & - ( blend_frequency == blend_frequency_outer & - .and. inner == inner_iterations ) & - .or. & - ( blend_frequency == blend_frequency_final & - .and. inner == inner_iterations & - .and. outer == outer_iterations ) ) then - - call lam_blend_lbc(self%state(igh_u), self%state(igh_p), & - self%state(igh_d), self%state(igh_t), & - mr, lbc_fields, microphysics_fields, & - aerosol_fields, lbc_option, & - microphysics_casim, murk_lbc, & - moisture_formulation, prime_mesh_name) - - end if - endif - - end do inner_dynamics_loop - !-------------------------------------------------------------------- - ! End of Inner (nonlinear, coriolis) loop - !-------------------------------------------------------------------- - - else ! when use_wavedynamics=false, just add physics increments here - - call invoke(X_plus_Y(self%state(igh_t), self%state_after_slow(igh_t), self%dtheta), & - X_plus_Y(self%state(igh_u), self%state_after_slow(igh_u), self%du)) - - end if ! use_wavedynamics - - end do outer_dynamics_loop - !-------------------------------------------------------------------- - ! End of Outer (advection) loop - !-------------------------------------------------------------------- - call adv_tracer_all_outer_after_slow%clear() - call adv_tracer_last_outer_after_slow%clear() - call con_tracer_all_outer_after_slow%clear() - call con_tracer_last_outer_after_slow%clear() - - if (transport_ageofair) then - call con_tracer_last_outer%get_field('ageofair', ageofair) - call ageofair_update(ageofair, model_clock) - end if - - !-------------------------------------------------------------------- - ! Apply mixing - !-------------------------------------------------------------------- - call mixing_alg(mr, self%state(igh_t), & - self%state(igh_u), derived_fields, & - self%state(igh_d), cast_dt ) - - ! ---------------------------------------------------------------------- - ! Call cloud scheme to generate cloud and latent heating after pressure - ! changes are applied from the solver. - ! ---------------------------------------------------------------------- -#ifdef UM_PHYSICS - if (use_physics .and. cloud == cloud_um ) then - call cld_alg( self%dtheta_cld, mr, & - self%state(igh_t), self%state(igh_p), self%state(igh_d), & - derived_fields, turbulence_fields, & - cloud_fields, convection_fields, & - self%state_n(igh_t), self%mr_n, & - model_clock%get_step(), cast_dt, .false. ) - call invoke(inc_X_plus_Y(self%state(igh_t), self%dtheta_cld)) - if (microphysics_casim) then - call casim_activate_alg( self%state(igh_t), mr, derived_fields, & - cloud_fields, microphysics_fields, & - convection_fields, initialise=.false. ) - end if ! microphysics_casim - end if -#endif - - ! Update derived variables for time level n+1 - if (self%use_moisture) then - call moist_dyn_factors_alg(moist_dyn, mr) - end if - if (use_physics) then - call map_physics_fields_alg(self%state(igh_u), self%state(igh_p), & - self%state(igh_d), self%state(igh_t), & - moist_dyn, derived_fields) - end if - - !-------------------------------------------------------------------- - ! Call UKCA for GLOMAP-mode prognostic aerosol updates - !-------------------------------------------------------------------- -#ifdef UM_PHYSICS - if ( aerosol == aerosol_um .and. & - ( glomap_mode == glomap_mode_ukca .or. & - glomap_mode == glomap_mode_dust_and_clim ) ) then - call aerosol_ukca_alg ( chemistry_fields, aerosol_fields, & - radiation_fields, derived_fields, & - microphysics_fields, turbulence_fields, & - convection_fields, cloud_fields, & - surface_fields, soil_fields, & - self%state, mr, model_clock ) - end if -#endif - - ! Write diagnostic output - if (use_physics .and. self%output_cld_incs) then - call cld_incs_output(cloud_fields, dcfl_tot, dcff_tot, dbcf_tot, & - sec_tot, suffix_tot) - end if - if (write_diag .and. use_xios_io .and. & - mod(model_clock%get_step(),diagnostic_frequency) == 0 ) then - call output_diags_for_si(self%state, self%state_n, self%state_after_slow,& - mr, self%mr_n, self%mr_after_slow, & - self%mr_after_adv, derived_fields, & - self%du, self%dtheta, self%dtheta_cld) - end if - - !-------------------------------------------------------------------- - ! Update fields held in the driver layer - !-------------------------------------------------------------------- - call invoke( setval_X(u, self%state(igh_u)), & - setval_X(theta, self%state(igh_t)), & - setval_X(rho, self%state(igh_d)), & - setval_X(exner, self%state(igh_p)) ) - - nullify( mm_wt, mm_vel ) - - if ( LPROF ) call stop_timing( id, 'semi_implicit_timestep' ) - - end subroutine run_step - - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !> Release all claimed resources once completed. - !> - subroutine semi_implicit_alg_final(self) - - implicit none - - class(semi_implicit_timestep_type), intent(inout) :: self - - call semi_implicit_solver_alg_final() - call final_si_operators() - - if (allocated(self%state)) deallocate(self%state) - if (allocated(self%state_n)) deallocate(self%state_n) - if (allocated(self%state_after_slow)) deallocate(self%state_after_slow) - if (allocated(self%advected_state)) deallocate(self%advected_state) - if (allocated(self%rhs_n)) deallocate(self%rhs_n) - if (allocated(self%rhs_np1)) deallocate(self%rhs_np1) - if (allocated(self%rhs_adv)) deallocate(self%rhs_adv) - if (allocated(self%rhs_phys)) deallocate(self%rhs_phys) - if (allocated(self%rhs_lbc)) deallocate(self%rhs_lbc) - if (allocated(self%adv_inc_prev)) deallocate(self%adv_inc_prev) - if (allocated(self%mr_n)) deallocate(self%mr_n) - if (allocated(self%mr_after_slow)) deallocate(self%mr_after_slow) - if (allocated(self%mr_after_adv)) deallocate(self%mr_after_adv) - - call self%dtheta%field_final() - call self%du%field_final() - call self%total_dry_flux%field_final() - - return - end subroutine semi_implicit_alg_final - - !============================================================================! - !> @brief Make a deep copy of a subset of a field collection, based on a - !> list of fields provided by another collection - !> @param[out] generic_fields_copied New collection written with fields - !> saved at some point in time - !> @param[in] generic_fields_to_copy collection we want to save fields from - !> @param[in] field_list list of fields to be saved - subroutine conditional_collection_copy(generic_fields_copied, & - generic_fields_to_copy, & - field_list) - - use field_collection_mod, only: field_collection_type - use field_collection_iterator_mod, only: field_collection_iterator_type - - implicit none - - type(field_collection_type), intent(out) :: generic_fields_copied - type(field_collection_type), intent(in) :: generic_fields_to_copy - type(field_collection_type), intent(in) :: field_list - - ! Iterator for field collection - type(field_collection_iterator_type) :: iterator - - ! One of the single fields out of the generic_fields_to_copy collection - class( field_parent_type ), pointer :: abstract_field_ptr - type(field_type), pointer :: single_generic_field - - ! The saved version of single_generic_field - type(field_type) :: copied_generic_field - - logical(kind=l_def) :: l_copy - - nullify(abstract_field_ptr) - nullify(single_generic_field) - - call generic_fields_copied%initialise(name='fields_copied') - - if ( generic_fields_to_copy%get_length() > 0 ) then - - call iterator%initialise(generic_fields_to_copy) - - do - if ( .not.iterator%has_next() ) exit - - abstract_field_ptr => iterator%next() - select type(abstract_field_ptr) - type is (field_type) - single_generic_field => abstract_field_ptr - end select - - l_copy = field_list%field_exists(single_generic_field%get_name()) - - if ( l_copy ) then - - ! We copy the field we want to save into a new field - call single_generic_field%copy_field_properties(copied_generic_field) - call invoke( setval_X(copied_generic_field, single_generic_field) ) - - ! We add it to the field collection passed in - call generic_fields_copied%add_field(copied_generic_field) - - end if ! ( l_copy ) - - end do - - end if - - end subroutine conditional_collection_copy -end module semi_implicit_timestep_alg_mod diff --git a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.f90 b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.f90 new file mode 100644 index 000000000..445799637 --- /dev/null +++ b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.f90 @@ -0,0 +1,633 @@ +MODULE semi_implicit_timestep_alg_mod + USE constants_mod, ONLY: i_def, r_def, l_def, str_def + USE fs_continuity_mod, ONLY: Wtheta, W2 + USE log_mod, ONLY: log_event, log_scratch_space, LOG_LEVEL_INFO + USE extrusion_mod, ONLY: TWOD + USE namelist_mod, ONLY: namelist_type + USE sci_fem_constants_mod, ONLY: get_mass_matrix_fe, get_mass_matrix_fv + USE sci_field_bundle_builtins_mod, ONLY: clone_bundle, bundle_axpy, add_bundle, copy_bundle, set_bundle_scalar + USE timestep_method_mod, ONLY: timestep_method_type + USE section_choice_config_mod, ONLY: cloud, cloud_um, aerosol, aerosol_um + USE physics_config_mod, ONLY: blayer_placement, blayer_placement_fast, convection_placement, convection_placement_fast, & +&stochastic_physics_placement, stochastic_physics_placement_fast, smagorinsky_placement, smagorinsky_placement_outer + USE aerosol_config_mod, ONLY: glomap_mode, glomap_mode_dust_and_clim, glomap_mode_ukca + USE formulation_config_mod, ONLY: use_physics, dlayer_on, use_wavedynamics, moisture_formulation, moisture_formulation_dry, & +&exner_from_eos + USE io_config_mod, ONLY: write_conservation_diag, write_diag, use_xios_io, diagnostic_frequency, checkpoint_read + USE initialization_config_mod, ONLY: init_option, init_option_checkpoint_dump, lbc_option_um2lfric_file + USE mixed_solver_config_mod, ONLY: guess_np1, reference_reset_time + USE timestepping_config_mod, ONLY: alpha, outer_iterations, inner_iterations, spinup_alpha + USE transport_config_mod, ONLY: cheap_update, transport_ageofair + USE derived_config_mod, ONLY: bundle_size + USE boundaries_config_mod, ONLY: limited_area, blend_frequency, blend_frequency_inner, blend_frequency_outer, & +&blend_frequency_final + USE finite_element_config_mod, ONLY: element_order_h, element_order_v + USE field_array_mod, ONLY: field_array_type + USE field_mod, ONLY: field_type + USE field_parent_mod, ONLY: field_parent_type + USE field_collection_mod, ONLY: field_collection_type + USE r_tran_field_mod, ONLY: r_tran_field_type + USE io_value_mod, ONLY: io_value_type, get_io_value + USE driver_modeldb_mod, ONLY: modeldb_type + USE mesh_mod, ONLY: mesh_type + USE mesh_collection_mod, ONLY: mesh_collection + USE model_clock_mod, ONLY: model_clock_type + USE operator_mod, ONLY: operator_type + USE rhs_alg_mod, ONLY: rhs_alg + USE gungho_transport_control_alg_mod, ONLY: gungho_transport_control_alg_init, gungho_transport_control_alg + USE si_operators_alg_mod, ONLY: create_si_operators, compute_si_operators, final_si_operators + USE fast_physics_alg_mod, ONLY: fast_physics + USE slow_physics_alg_mod, ONLY: slow_physics + USE checks_and_balances_alg_mod, ONLY: check_fields + USE semi_implicit_solver_alg_mod, ONLY: semi_implicit_solver_alg_init, semi_implicit_solver_alg_step, & +&semi_implicit_solver_alg_final + USE derive_exner_from_eos_alg_mod, ONLY: derive_exner_from_eos + USE sci_mass_matrix_solver_alg_mod, ONLY: mass_matrix_solver_alg + USE moist_dyn_factors_alg_mod, ONLY: moist_dyn_factors_alg + USE update_prognostic_scalars_alg_mod, ONLY: update_prognostic_scalars_alg + USE mixing_alg_mod, ONLY: mixing_alg + USE si_diagnostics_mod, ONLY: output_diags_for_si + USE predictors_alg_mod, ONLY: predictors_alg + USE limited_area_lbc_alg_mod, ONLY: lam_solver_lbc, lam_blend_lbc + USE lam_rhs_alg_mod, ONLY: calc_rhs_lbc, apply_mask_rhs + USE calc_phys_predictors_alg_mod, ONLY: calc_phys_predictors_alg + USE map_physics_fields_alg_mod, ONLY: map_physics_fields_alg + USE mr_indices_mod, ONLY: nummr + USE moist_dyn_mod, ONLY: num_moist_factors, gas_law + USE field_indices_mod, ONLY: igh_u, igh_t, igh_d, igh_p + USE mixing_config_mod, ONLY: smagorinsky + USE smagorinsky_alg_mod, ONLY: smagorinsky_alg + USE cld_incs_mod, ONLY: cld_incs_init, cld_incs_output + USE timing_mod, ONLY: start_timing, stop_timing, tik, LPROF + USE ageofair_alg_mod, ONLY: ageofair_update + IMPLICIT NONE + PRIVATE + TYPE, EXTENDS(timestep_method_type), PUBLIC :: semi_implicit_timestep_type + PRIVATE + LOGICAL(KIND = l_def) :: use_moisture + LOGICAL(KIND = l_def) :: output_cld_incs + TYPE(field_type), ALLOCATABLE :: state(:) + TYPE(field_type), ALLOCATABLE :: state_n(:) + TYPE(field_type), ALLOCATABLE :: state_after_slow(:) + TYPE(field_type), ALLOCATABLE :: advected_state(:) + TYPE(field_type), ALLOCATABLE :: mr_n(:), mr_after_adv(:), mr_after_slow(:) + TYPE(field_type), ALLOCATABLE :: rhs_n(:), rhs_np1(:), rhs_adv(:) + TYPE(field_type), ALLOCATABLE :: adv_inc_prev(:) + TYPE(field_type), ALLOCATABLE :: rhs_phys(:), rhs_lbc(:) + TYPE(field_type) :: dtheta, dtheta_cld + TYPE(field_type) :: du + TYPE(field_type) :: wind_prev + TYPE(r_tran_field_type) :: total_dry_flux + CONTAINS + PRIVATE + PROCEDURE, PUBLIC :: step => semi_implicit_alg_step + PROCEDURE, PUBLIC :: finalise => semi_implicit_alg_final + PROCEDURE, NOPASS :: run_init + PROCEDURE, NOPASS :: run_step + PROCEDURE, NOPASS :: conditional_collection_copy + END TYPE semi_implicit_timestep_type + INTERFACE semi_implicit_timestep_type + MODULE PROCEDURE semi_implicit_alg_init + END INTERFACE semi_implicit_timestep_type + CONTAINS + FUNCTION semi_implicit_alg_init(modeldb) RESULT(self) + IMPLICIT NONE + TYPE(semi_implicit_timestep_type) :: self + TYPE(modeldb_type), INTENT(IN), TARGET :: modeldb + TYPE(field_collection_type), POINTER :: prognostic_fields + TYPE(field_collection_type), POINTER :: moisture_fields + TYPE(field_type), POINTER :: u + TYPE(field_type), POINTER :: rho + TYPE(field_type), POINTER :: theta + TYPE(field_type), POINTER :: exner + TYPE(field_array_type), POINTER :: mr_array + TYPE(field_type), POINTER :: mr(:) + CLASS(model_clock_type), POINTER :: model_clock + model_clock => modeldb % clock + prognostic_fields => modeldb % fields % get_field_collection("prognostic_fields") + CALL prognostic_fields % get_field('theta', theta) + CALL prognostic_fields % get_field('u', u) + CALL prognostic_fields % get_field('rho', rho) + CALL prognostic_fields % get_field('exner', exner) + moisture_fields => modeldb % fields % get_field_collection("moisture_fields") + CALL moisture_fields % get_field("mr", mr_array) + mr => mr_array % bundle + CALL run_init(self, u, rho, theta, exner, mr, prognostic_fields, moisture_fields, model_clock) + NULLIFY(prognostic_fields, moisture_fields, mr_array, model_clock, u, rho, theta, exner, mr) + END FUNCTION semi_implicit_alg_init + SUBROUTINE semi_implicit_alg_step(self, modeldb) + IMPLICIT NONE + CLASS(semi_implicit_timestep_type), INTENT(INOUT) :: self + TYPE(modeldb_type), INTENT(IN), TARGET :: modeldb + TYPE(mesh_type), POINTER :: mesh + TYPE(mesh_type), POINTER :: twod_mesh + CLASS(model_clock_type), POINTER :: model_clock + TYPE(field_collection_type), POINTER :: prognostic_fields + TYPE(field_collection_type), POINTER :: moisture_fields + TYPE(field_type), POINTER :: u + TYPE(field_type), POINTER :: rho + TYPE(field_type), POINTER :: theta + TYPE(field_type), POINTER :: exner + TYPE(field_array_type), POINTER :: mr_array + TYPE(field_type), POINTER :: mr(:) + TYPE(field_array_type), POINTER :: moist_dyn_array + TYPE(field_type), POINTER :: moist_dyn(:) + TYPE(field_collection_type), POINTER :: adv_tracer_all_outer + TYPE(field_collection_type), POINTER :: adv_tracer_last_outer + TYPE(field_collection_type), POINTER :: con_tracer_all_outer + TYPE(field_collection_type), POINTER :: con_tracer_last_outer + TYPE(field_collection_type), POINTER :: derived_fields + TYPE(field_collection_type), POINTER :: radiation_fields + TYPE(field_collection_type), POINTER :: microphysics_fields + TYPE(field_collection_type), POINTER :: electric_fields + TYPE(field_collection_type), POINTER :: orography_fields + TYPE(field_collection_type), POINTER :: turbulence_fields + TYPE(field_collection_type), POINTER :: convection_fields + TYPE(field_collection_type), POINTER :: cloud_fields + TYPE(field_collection_type), POINTER :: surface_fields + TYPE(field_collection_type), POINTER :: soil_fields + TYPE(field_collection_type), POINTER :: snow_fields + TYPE(field_collection_type), POINTER :: chemistry_fields + TYPE(field_collection_type), POINTER :: aerosol_fields + TYPE(field_collection_type), POINTER :: stph_fields + TYPE(field_collection_type), POINTER :: lbc_fields + TYPE(io_value_type), POINTER :: temp_corr_io_value + REAL(KIND = r_def) :: dt + REAL(KIND = r_def) :: dtemp_encorr + prognostic_fields => modeldb % fields % get_field_collection("prognostic_fields") + model_clock => modeldb % clock + CALL prognostic_fields % get_field('theta', theta) + CALL prognostic_fields % get_field('u', u) + CALL prognostic_fields % get_field('rho', rho) + CALL prognostic_fields % get_field('exner', exner) + dt = REAL(model_clock % get_seconds_per_step(), r_def) + moisture_fields => modeldb % fields % get_field_collection("moisture_fields") + lbc_fields => modeldb % fields % get_field_collection("lbc_fields") + radiation_fields => modeldb % fields % get_field_collection("radiation_fields") + CALL moisture_fields % get_field("mr", mr_array) + mr => mr_array % bundle + CALL moisture_fields % get_field("moist_dyn", moist_dyn_array) + moist_dyn => moist_dyn_array % bundle + mesh => theta % get_mesh() + twod_mesh => mesh_collection % get_mesh(mesh, TWOD) + adv_tracer_all_outer => modeldb % fields % get_field_collection("adv_tracer_all_outer") + adv_tracer_last_outer => modeldb % fields % get_field_collection("adv_tracer_last_outer") + con_tracer_all_outer => modeldb % fields % get_field_collection("con_tracer_all_outer") + con_tracer_last_outer => modeldb % fields % get_field_collection("con_tracer_last_outer") + derived_fields => modeldb % fields % get_field_collection("derived_fields") + microphysics_fields => modeldb % fields % get_field_collection("microphysics_fields") + turbulence_fields => modeldb % fields % get_field_collection("turbulence_fields") + convection_fields => modeldb % fields % get_field_collection("convection_fields") + cloud_fields => modeldb % fields % get_field_collection("cloud_fields") + surface_fields => modeldb % fields % get_field_collection("surface_fields") + soil_fields => modeldb % fields % get_field_collection("soil_fields") + snow_fields => modeldb % fields % get_field_collection("snow_fields") + chemistry_fields => modeldb % fields % get_field_collection("chemistry_fields") + aerosol_fields => modeldb % fields % get_field_collection("aerosol_fields") + stph_fields => modeldb % fields % get_field_collection("stph_fields") + electric_fields => modeldb % fields % get_field_collection("electric_fields") + orography_fields => modeldb % fields % get_field_collection("orography_fields") + temp_corr_io_value => get_io_value(modeldb % values, 'temperature_correction_io_value') + dtemp_encorr = dt * temp_corr_io_value % data(1) + CALL run_step(self, modeldb, u, rho, theta, exner, mr, moist_dyn, adv_tracer_all_outer, adv_tracer_last_outer, & +&con_tracer_all_outer, con_tracer_last_outer, prognostic_fields, moisture_fields, derived_fields, radiation_fields, & +µphysics_fields, electric_fields, orography_fields, turbulence_fields, convection_fields, cloud_fields, surface_fields, & +&soil_fields, snow_fields, chemistry_fields, aerosol_fields, stph_fields, lbc_fields, model_clock, dtemp_encorr, mesh, twod_mesh) + NULLIFY(model_clock, mesh, prognostic_fields, moisture_fields, u, rho, theta, exner, mr, moist_dyn, adv_tracer_all_outer, & +&adv_tracer_last_outer, con_tracer_all_outer, con_tracer_last_outer, derived_fields, radiation_fields, microphysics_fields, & +&electric_fields, orography_fields, turbulence_fields, convection_fields, cloud_fields, surface_fields, soil_fields, snow_fields, & +&chemistry_fields, aerosol_fields, stph_fields, lbc_fields, temp_corr_io_value, mr_array, moist_dyn_array) + END SUBROUTINE semi_implicit_alg_step + SUBROUTINE run_init(self, u, rho, theta, exner, mr, prognostic_fields, moisture_fields, model_clock) + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_0 + IMPLICIT NONE + TYPE(semi_implicit_timestep_type), INTENT(INOUT) :: self + TYPE(field_type), INTENT(IN) :: u, rho, theta, exner + TYPE(field_type), INTENT(IN), OPTIONAL :: mr(:) + TYPE(field_collection_type), INTENT(INOUT) :: prognostic_fields + TYPE(field_collection_type), INTENT(INOUT) :: moisture_fields + CLASS(model_clock_type), INTENT(IN) :: model_clock + TYPE(mesh_type), POINTER :: mesh + TYPE(field_type), POINTER :: rho_ref, theta_ref, exner_ref + TYPE(field_array_type), POINTER :: moist_dyn_ref_array + TYPE(field_type), POINTER :: moist_dyn_ref(:) + REAL(KIND = r_def) :: cast_dt + INTEGER(KIND = i_def) :: reference_reset_freq + mesh => theta % get_mesh() + ALLOCATE(self % state(bundle_size)) + ALLOCATE(self % state_n(bundle_size)) + ALLOCATE(self % state_after_slow(bundle_size)) + ALLOCATE(self % advected_state(bundle_size)) + ALLOCATE(self % rhs_n(bundle_size)) + ALLOCATE(self % rhs_np1(bundle_size)) + ALLOCATE(self % rhs_adv(bundle_size)) + ALLOCATE(self % adv_inc_prev(bundle_size)) + ALLOCATE(self % rhs_phys(bundle_size)) + ALLOCATE(self % rhs_lbc(bundle_size)) + ALLOCATE(self % mr_n(nummr)) + ALLOCATE(self % mr_after_slow(nummr)) + ALLOCATE(self % mr_after_adv(nummr)) + self % use_moisture = (moisture_formulation /= moisture_formulation_dry) + self % output_cld_incs = (self % use_moisture .AND. write_diag .AND. use_xios_io) + CALL u % copy_field_properties(self % state(igh_u)) + CALL theta % copy_field_properties(self % state(igh_t)) + CALL rho % copy_field_properties(self % state(igh_d)) + CALL exner % copy_field_properties(self % state(igh_p)) + CALL clone_bundle(self % state, self % state_n, bundle_size) + CALL clone_bundle(self % state, self % state_after_slow, bundle_size) + CALL clone_bundle(self % state, self % advected_state, bundle_size) + CALL clone_bundle(self % state, self % rhs_n, bundle_size) + CALL clone_bundle(self % state, self % rhs_np1, bundle_size) + CALL clone_bundle(self % state, self % rhs_adv, bundle_size) + CALL clone_bundle(self % state, self % rhs_phys, bundle_size) + CALL clone_bundle(self % state, self % rhs_lbc, bundle_size) + CALL clone_bundle(self % state, self % adv_inc_prev, bundle_size) + CALL theta % copy_field_properties(self % dtheta) + CALL theta % copy_field_properties(self % dtheta_cld) + CALL invoke_0(self % dtheta_cld) + CALL u % copy_field_properties(self % du) + CALL u % copy_field_properties(self % wind_prev) + CALL self % total_dry_flux % initialise(u % get_function_space()) + CALL clone_bundle(mr, self % mr_n, nummr) + CALL clone_bundle(mr, self % mr_after_slow, nummr) + IF (self % use_moisture) THEN + CALL clone_bundle(mr, self % mr_after_adv, nummr) + ELSE + CALL set_bundle_scalar(0.0_r_def, self % mr_n, nummr) + CALL set_bundle_scalar(0.0_r_def, self % mr_after_slow, nummr) + END IF + CALL set_bundle_scalar(0.0_r_def, self % rhs_phys, bundle_size) + CALL create_si_operators(mesh) + cast_dt = REAL(model_clock % get_seconds_per_step(), r_def) + reference_reset_freq = NINT(reference_reset_time / cast_dt, i_def) + IF (MOD(model_clock % get_first_step() - 1, reference_reset_freq) /= 0) THEN + CALL moisture_fields % get_field("moist_dyn_ref", moist_dyn_ref_array) + moist_dyn_ref => moist_dyn_ref_array % bundle + CALL prognostic_fields % get_field('theta_ref', theta_ref) + CALL prognostic_fields % get_field('rho_ref', rho_ref) + CALL prognostic_fields % get_field('exner_ref', exner_ref) + CALL compute_si_operators(theta_ref, rho_ref, exner_ref, model_clock, moist_dyn_ref) + NULLIFY(theta_ref, rho_ref, exner_ref, moist_dyn_ref, moist_dyn_ref_array) + END IF + IF (use_wavedynamics) THEN + CALL gungho_transport_control_alg_init(mesh) + CALL semi_implicit_solver_alg_init(self % state) + END IF + NULLIFY(mesh) + CALL log_event("semi_implicit_timestep: initialised timestepping algorithm", LOG_LEVEL_INFO) + END SUBROUTINE run_init + SUBROUTINE run_step(self, modeldb, u, rho, theta, exner, mr, moist_dyn, adv_tracer_all_outer, adv_tracer_last_outer, & +&con_tracer_all_outer, con_tracer_last_outer, prognostic_fields, moisture_fields, derived_fields, radiation_fields, & +µphysics_fields, electric_fields, orography_fields, turbulence_fields, convection_fields, cloud_fields, surface_fields, & +&soil_fields, snow_fields, chemistry_fields, aerosol_fields, stph_fields, lbc_fields, model_clock, dtemp_encorr, mesh, twod_mesh) + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_11 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_10 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_9 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_update_rhs_phys_from_fast_physics + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_7 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_6 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_5 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_4 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_update_from_slow_physics + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_2 + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_copy_init_fields_to_state + IMPLICIT NONE + TYPE(semi_implicit_timestep_type), INTENT(INOUT), TARGET :: self + TYPE(modeldb_type), INTENT(IN), TARGET :: modeldb + TYPE(field_type), INTENT(INOUT) :: u, rho, theta, exner + TYPE(field_type), INTENT(INOUT) :: mr(nummr) + TYPE(field_type), INTENT(INOUT) :: moist_dyn(num_moist_factors) + TYPE(field_collection_type), INTENT(INOUT) :: adv_tracer_all_outer + TYPE(field_collection_type), INTENT(INOUT) :: adv_tracer_last_outer + TYPE(field_collection_type), INTENT(INOUT) :: con_tracer_all_outer + TYPE(field_collection_type), INTENT(INOUT) :: con_tracer_last_outer + TYPE(field_collection_type), INTENT(INOUT) :: prognostic_fields + TYPE(field_collection_type), INTENT(INOUT) :: moisture_fields + TYPE(field_collection_type), INTENT(INOUT) :: derived_fields + TYPE(field_collection_type), INTENT(INOUT) :: radiation_fields + TYPE(field_collection_type), INTENT(INOUT) :: microphysics_fields + TYPE(field_collection_type), INTENT(INOUT) :: electric_fields + TYPE(field_collection_type), INTENT(IN) :: orography_fields + TYPE(field_collection_type), INTENT(INOUT) :: turbulence_fields + TYPE(field_collection_type), INTENT(INOUT) :: convection_fields + TYPE(field_collection_type), INTENT(INOUT) :: cloud_fields + TYPE(field_collection_type), INTENT(INOUT) :: surface_fields + TYPE(field_collection_type), INTENT(INOUT) :: soil_fields + TYPE(field_collection_type), INTENT(INOUT) :: snow_fields + TYPE(field_collection_type), INTENT(INOUT) :: chemistry_fields + TYPE(field_collection_type), INTENT(INOUT) :: aerosol_fields + TYPE(field_collection_type), INTENT(INOUT) :: stph_fields + TYPE(field_collection_type), INTENT(INOUT) :: lbc_fields + CLASS(model_clock_type), INTENT(IN) :: model_clock + REAL(KIND = r_def), INTENT(IN) :: dtemp_encorr + TYPE(mesh_type), INTENT(IN), POINTER :: mesh + TYPE(mesh_type), INTENT(IN), POINTER :: twod_mesh + TYPE(field_type), POINTER :: rho_ref + TYPE(field_type), POINTER :: theta_ref + TYPE(field_type), POINTER :: exner_ref + TYPE(field_type), POINTER :: moist_dyn_ref(:) + TYPE(field_array_type), POINTER :: moist_dyn_ref_array + TYPE(field_type), POINTER :: mr_to_adv(:) + TYPE(field_type) :: dcfl_tot, dcff_tot, dbcf_tot, dcfl_adv, dcff_adv, dbcf_adv + CHARACTER(LEN = str_def), PARAMETER :: sec_tot = 'processed' + CHARACTER(LEN = str_def), PARAMETER :: suffix_tot = 'tot' + CHARACTER(LEN = str_def), PARAMETER :: sec_adv = 'advection' + CHARACTER(LEN = str_def), PARAMETER :: suffix_adv = 'adv' + REAL(KIND = r_def) :: cast_dt + TYPE(field_type), POINTER :: ageofair + TYPE(operator_type), POINTER :: mm_wt + TYPE(operator_type), POINTER :: mm_vel + INTEGER(KIND = i_def) :: outer, inner, reference_reset_freq + REAL(KIND = r_def) :: varalpha, varbeta + LOGICAL(KIND = l_def) :: write_moisture_diag + TYPE(field_collection_type) :: adv_tracer_all_outer_after_slow + TYPE(field_collection_type) :: adv_tracer_last_outer_after_slow + TYPE(field_collection_type) :: con_tracer_all_outer_after_slow + TYPE(field_collection_type) :: con_tracer_last_outer_after_slow + LOGICAL(KIND = l_def) :: checkpoint_reference_fields + TYPE(namelist_type), POINTER :: base_mesh_nml + TYPE(namelist_type), POINTER :: initialization_nml + TYPE(namelist_type), POINTER :: microphysics_nml + TYPE(namelist_type), POINTER :: aerosol_nml + TYPE(namelist_type), POINTER :: timestepping_nml + CHARACTER(LEN = str_def) :: prime_mesh_name + INTEGER(KIND = i_def) :: lbc_option + LOGICAL(KIND = l_def) :: microphysics_casim + LOGICAL(KIND = l_def) :: murk_lbc + REAL(KIND = r_def) :: tau_r + INTEGER(KIND = tik) :: id + IF (LPROF) CALL start_timing(id, 'semi_implicit_timestep') + cast_dt = REAL(model_clock % get_seconds_per_step(), r_def) + IF (limited_area .AND. use_wavedynamics) THEN + base_mesh_nml => modeldb % configuration % get_namelist('base_mesh') + initialization_nml => modeldb % configuration % get_namelist('initialization') + timestepping_nml => modeldb % configuration % get_namelist('timestepping') + CALL base_mesh_nml % get_value('prime_mesh_name', prime_mesh_name) + CALL initialization_nml % get_value('lbc_option', lbc_option) + CALL timestepping_nml % get_value('tau_r', tau_r) + IF (lbc_option == lbc_option_um2lfric_file) THEN + aerosol_nml => modeldb % configuration % get_namelist('aerosol') + CALL aerosol_nml % get_value('murk_lbc', murk_lbc) + END IF + END IF + IF (lbc_option == lbc_option_um2lfric_file .OR. (use_physics .AND. cloud == cloud_um)) THEN + microphysics_nml => modeldb % configuration % get_namelist('microphysics') + CALL microphysics_nml % get_value('microphysics_casim', microphysics_casim) + END IF + IF (element_order_h == 0 .AND. element_order_v == 0) THEN + mm_wt => get_mass_matrix_fv(Wtheta, mesh % get_id()) + mm_vel => get_mass_matrix_fv(W2, mesh % get_id()) + ELSE + mm_wt => get_mass_matrix_fe(Wtheta, mesh % get_id()) + mm_vel => get_mass_matrix_fe(W2, mesh % get_id()) + END IF + CALL invoke_copy_init_fields_to_state(self % state(igh_u), u, self % state(igh_t), theta, self % state(igh_d), rho, & +&self % state(igh_p), exner, self % total_dry_flux) + reference_reset_freq = NINT(reference_reset_time / cast_dt, i_def) + IF (MOD(model_clock % get_step() - 1_i_def, reference_reset_freq) == 0_i_def) THEN + CALL compute_si_operators(self % state(igh_t), self % state(igh_d), self % state(igh_p), model_clock, moist_dyn) + checkpoint_reference_fields = MOD(model_clock % get_first_step() - 1, reference_reset_freq) /= 0 .OR. MOD(model_clock % & +&get_last_step(), reference_reset_freq) /= 0 + IF (checkpoint_read .OR. init_option == init_option_checkpoint_dump) THEN + IF (MOD(model_clock % get_first_step() - 1, reference_reset_freq) == 0 .AND. MOD(model_clock % get_last_step(), & +&reference_reset_freq) /= 0) THEN + checkpoint_reference_fields = .FALSE. + END IF + END IF + IF (checkpoint_reference_fields) THEN + CALL prognostic_fields % get_field('theta_ref', theta_ref) + CALL prognostic_fields % get_field('rho_ref', rho_ref) + CALL prognostic_fields % get_field('exner_ref', exner_ref) + CALL moisture_fields % get_field("moist_dyn_ref", moist_dyn_ref_array) + moist_dyn_ref => moist_dyn_ref_array % bundle + CALL copy_bundle(moist_dyn, moist_dyn_ref, num_moist_factors) + CALL invoke_2(theta_ref, self % state(igh_t), rho_ref, self % state(igh_d), exner_ref, self % state(igh_p)) + NULLIFY(moist_dyn_ref_array, moist_dyn_ref, theta_ref, rho_ref, exner_ref) + END IF + END IF + IF (spinup_alpha .AND. model_clock % is_spinning_up()) THEN + varalpha = 1.0_r_def + ELSE + varalpha = alpha + END IF + varbeta = 1.0_r_def - varalpha + CALL check_fields(self % state, cast_dt) + IF (self % use_moisture) THEN + CALL copy_bundle(mr, self % mr_n, nummr) + CALL copy_bundle(mr, self % mr_after_slow, nummr) + END IF + CALL copy_bundle(self % state, self % state_n, bundle_size) + CALL copy_bundle(self % state, self % state_after_slow, bundle_size) + IF (use_physics) THEN + IF (self % output_cld_incs) THEN + CALL cld_incs_init(cloud_fields, dcfl_tot, dcff_tot, dbcf_tot, sec_tot, suffix_tot) + END IF + CALL slow_physics(modeldb, self % du, self % dtheta, self % mr_after_slow, self % state_n(igh_t), self % state_n(igh_u), & +&self % state_n(igh_d), self % state_n(igh_p), moist_dyn, self % mr_n, derived_fields, radiation_fields, microphysics_fields, & +&electric_fields, orography_fields, turbulence_fields, convection_fields, cloud_fields, surface_fields, soil_fields, snow_fields, & +&chemistry_fields, aerosol_fields, model_clock, cast_dt, dtemp_encorr, mesh, twod_mesh) + CALL invoke_update_from_slow_physics(self % state_after_slow(igh_t), self % dtheta, self % state_after_slow(igh_u), self % du) + IF (self % output_cld_incs) THEN + CALL cld_incs_init(cloud_fields, dcfl_adv, dcff_adv, dbcf_adv, sec_adv, suffix_adv) + END IF + END IF + CALL rhs_alg(self % rhs_n, varbeta * cast_dt, self % state_after_slow, self % state_n, moist_dyn, compute_eos = .FALSE., & +&compute_rhs_t_d = .TRUE., dlayer_rhs = .FALSE., model_clock = model_clock) + CALL copy_bundle(self % state_after_slow, self % advected_state, bundle_size) + mr_to_adv => self % mr_after_slow + CALL predictors_alg(self % advected_state, self % state_n(igh_u), self % rhs_n(igh_u), varbeta, model_clock) + CALL conditional_collection_copy(adv_tracer_all_outer_after_slow, generic_fields_to_copy = adv_tracer_all_outer, & +&field_list = adv_tracer_all_outer) + CALL conditional_collection_copy(adv_tracer_last_outer_after_slow, generic_fields_to_copy = adv_tracer_last_outer, & +&field_list = adv_tracer_last_outer) + CALL conditional_collection_copy(con_tracer_all_outer_after_slow, generic_fields_to_copy = con_tracer_all_outer, & +&field_list = con_tracer_all_outer) + CALL conditional_collection_copy(con_tracer_last_outer_after_slow, generic_fields_to_copy = con_tracer_last_outer, & +&field_list = con_tracer_last_outer) + CALL invoke_4(self % wind_prev, self % state(igh_u)) + outer_dynamics_loop:DO outer = 1, outer_iterations + IF (use_wavedynamics) THEN + CALL gungho_transport_control_alg(self % rhs_adv, self % advected_state, self % state(igh_u), self % state_n(igh_u), mr, & +&mr_to_adv, model_clock, outer, cheap_update, self % adv_inc_prev, self % wind_prev, self % state_after_slow(igh_d), & +&self % total_dry_flux, adv_tracer_all_outer, adv_tracer_all_outer_after_slow, adv_tracer_last_outer, & +&adv_tracer_last_outer_after_slow, con_tracer_all_outer, con_tracer_all_outer_after_slow, con_tracer_last_outer, & +&con_tracer_last_outer_after_slow) + IF (cheap_update .AND. (outer < outer_iterations)) THEN + CALL copy_bundle(self % rhs_adv, self % adv_inc_prev, bundle_size) + CALL mass_matrix_solver_alg(self % du, self % rhs_adv(igh_u)) + CALL invoke_5(self % advected_state(igh_d), self % rhs_adv(igh_d), self % advected_state(igh_t), self % rhs_adv(igh_t), & +&self % advected_state(igh_u), self % du) + CALL adv_tracer_all_outer_after_slow % clear + CALL con_tracer_all_outer_after_slow % clear + CALL conditional_collection_copy(adv_tracer_all_outer_after_slow, generic_fields_to_copy = adv_tracer_all_outer, & +&field_list = adv_tracer_all_outer) + CALL conditional_collection_copy(con_tracer_all_outer_after_slow, generic_fields_to_copy = con_tracer_all_outer, & +&field_list = con_tracer_all_outer) + IF (self % use_moisture) THEN + mr_to_adv => self % mr_after_adv + END IF + CALL invoke_6(self % wind_prev, self % state(igh_u)) + END IF + CALL invoke_7(self % dtheta, self % rhs_adv(igh_t), mm_wt) + CALL rhs_alg(self % rhs_np1, - varalpha * cast_dt, self % state, self % state, moist_dyn, compute_eos = .TRUE., & +&compute_rhs_t_d = .TRUE., dlayer_rhs = dlayer_on, model_clock = model_clock) + ELSE + IF (self % use_moisture) CALL copy_bundle(self % mr_after_slow, mr, nummr) + END IF + IF (self % use_moisture) CALL copy_bundle(mr, self % mr_after_adv, nummr) + IF (use_physics) THEN + IF (blayer_placement == blayer_placement_fast .OR. convection_placement == convection_placement_fast .OR. & +&stochastic_physics_placement == stochastic_physics_placement_fast) THEN + CALL calc_phys_predictors_alg(derived_fields, self % rhs_np1, self % rhs_adv, self % rhs_n, self % state, & +&self % state_after_slow, lbc_fields, model_clock) + END IF + IF (outer == outer_iterations .AND. self % output_cld_incs) THEN + CALL cld_incs_output(cloud_fields, dcfl_adv, dcff_adv, dbcf_adv, sec_adv, suffix_adv) + END IF + CALL fast_physics(self % du, self % dtheta, mr, self % state_n(igh_t), self % state_n(igh_d), self % state_n(igh_u), & +&self % state_n(igh_p), self % mr_n, derived_fields, radiation_fields, microphysics_fields, orography_fields, turbulence_fields, & +&convection_fields, cloud_fields, surface_fields, soil_fields, snow_fields, chemistry_fields, aerosol_fields, stph_fields, outer, & +&model_clock, cast_dt) + IF (smagorinsky .AND. smagorinsky_placement == smagorinsky_placement_outer) THEN + CALL smagorinsky_alg(self % dtheta, self % du, mr, self % state(igh_t), self % state(igh_u), derived_fields, & +&self % state(igh_d), cast_dt) + END IF + IF (use_wavedynamics) THEN + CALL set_bundle_scalar(0.0_r_def, self % rhs_phys, bundle_size) + CALL invoke_update_rhs_phys_from_fast_physics(self % rhs_phys(igh_t), self % dtheta, mm_wt, self % rhs_phys(igh_u), & +&self % du, mm_vel) + END IF + END IF + IF (use_wavedynamics) THEN + IF (guess_np1) THEN + IF (self % use_moisture) CALL moist_dyn_factors_alg(moist_dyn, mr) + CALL update_prognostic_scalars_alg(self % state, self % rhs_n, self % rhs_adv, self % rhs_phys, moist_dyn(gas_law)) + END IF + inner_dynamics_loop:DO inner = 1, inner_iterations + WRITE(log_scratch_space, '(A,2I3)') 'loop indices (o, i): ', outer, inner + CALL log_event(log_scratch_space, LOG_LEVEL_INFO) + IF (inner > 1) CALL rhs_alg(self % rhs_np1, - varalpha * cast_dt, self % state, self % state, moist_dyn, compute_eos = & +&.TRUE., compute_rhs_t_d = .FALSE., dlayer_rhs = dlayer_on, model_clock = model_clock) + IF (limited_area .AND. inner == 1 .AND. outer == 1) THEN + CALL lam_solver_lbc(self % state(igh_u), lbc_fields, prime_mesh_name) + CALL calc_rhs_lbc(self % rhs_lbc, lbc_fields, model_clock, prime_mesh_name, tau_r) + END IF + CALL bundle_axpy(- 1.0_r_def, self % rhs_np1, self % rhs_n, self % rhs_np1, bundle_size) + CALL add_bundle(self % rhs_np1, self % rhs_adv, self % rhs_np1, bundle_size) + CALL add_bundle(self % rhs_np1, self % rhs_phys, self % rhs_np1, bundle_size) + IF (limited_area) THEN + IF (inner == 1 .AND. outer == 1) THEN + CALL add_bundle(self % rhs_np1, self % rhs_lbc, self % rhs_np1, bundle_size) + END IF + CALL apply_mask_rhs(self % rhs_np1, prime_mesh_name) + END IF + IF (inner > 1) THEN + CALL invoke_9(self % rhs_np1(igh_d), self % rhs_np1(igh_t)) + END IF + write_moisture_diag = write_conservation_diag .AND. outer == outer_iterations .AND. inner == inner_iterations .AND. self & +&% use_moisture + CALL semi_implicit_solver_alg_step(self % state, self % rhs_np1, moist_dyn(gas_law), mr, write_moisture_diag, & +&first_iteration = (inner == 1)) + IF (.NOT. guess_np1 .AND. self % use_moisture) CALL moist_dyn_factors_alg(moist_dyn, mr) + IF (exner_from_eos) THEN + CALL derive_exner_from_eos(self % state, moist_dyn(gas_law)) + END IF + IF (limited_area) THEN + IF ((blend_frequency == blend_frequency_inner) .OR. (blend_frequency == blend_frequency_outer .AND. inner == & +&inner_iterations) .OR. (blend_frequency == blend_frequency_final .AND. inner == inner_iterations .AND. outer == & +&outer_iterations)) THEN + CALL lam_blend_lbc(self % state(igh_u), self % state(igh_p), self % state(igh_d), self % state(igh_t), mr, & +&lbc_fields, microphysics_fields, aerosol_fields, lbc_option, microphysics_casim, murk_lbc, moisture_formulation, prime_mesh_name) + END IF + END IF + END DO inner_dynamics_loop + ELSE + CALL invoke_10(self % state(igh_t), self % state_after_slow(igh_t), self % dtheta, self % state(igh_u), & +&self % state_after_slow(igh_u), self % du) + END IF + END DO outer_dynamics_loop + CALL adv_tracer_all_outer_after_slow % clear + CALL adv_tracer_last_outer_after_slow % clear + CALL con_tracer_all_outer_after_slow % clear + CALL con_tracer_last_outer_after_slow % clear + IF (transport_ageofair) THEN + CALL con_tracer_last_outer % get_field('ageofair', ageofair) + CALL ageofair_update(ageofair, model_clock) + END IF + CALL mixing_alg(mr, self % state(igh_t), self % state(igh_u), derived_fields, self % state(igh_d), cast_dt) + IF (self % use_moisture) THEN + CALL moist_dyn_factors_alg(moist_dyn, mr) + END IF + IF (use_physics) THEN + CALL map_physics_fields_alg(self % state(igh_u), self % state(igh_p), self % state(igh_d), self % state(igh_t), moist_dyn, & +&derived_fields) + END IF + IF (use_physics .AND. self % output_cld_incs) THEN + CALL cld_incs_output(cloud_fields, dcfl_tot, dcff_tot, dbcf_tot, sec_tot, suffix_tot) + END IF + IF (write_diag .AND. use_xios_io .AND. MOD(model_clock % get_step(), diagnostic_frequency) == 0) THEN + CALL output_diags_for_si(self % state, self % state_n, self % state_after_slow, mr, self % mr_n, self % mr_after_slow, & +&self % mr_after_adv, derived_fields, self % du, self % dtheta, self % dtheta_cld) + END IF + CALL invoke_11(u, self % state(igh_u), theta, self % state(igh_t), rho, self % state(igh_d), exner, self % state(igh_p)) + NULLIFY(mm_wt, mm_vel) + IF (LPROF) CALL stop_timing(id, 'semi_implicit_timestep') + END SUBROUTINE run_step + SUBROUTINE semi_implicit_alg_final(self) + IMPLICIT NONE + CLASS(semi_implicit_timestep_type), INTENT(INOUT) :: self + CALL semi_implicit_solver_alg_final + CALL final_si_operators + IF (ALLOCATED(self % state)) DEALLOCATE(self % state) + IF (ALLOCATED(self % state_n)) DEALLOCATE(self % state_n) + IF (ALLOCATED(self % state_after_slow)) DEALLOCATE(self % state_after_slow) + IF (ALLOCATED(self % advected_state)) DEALLOCATE(self % advected_state) + IF (ALLOCATED(self % rhs_n)) DEALLOCATE(self % rhs_n) + IF (ALLOCATED(self % rhs_np1)) DEALLOCATE(self % rhs_np1) + IF (ALLOCATED(self % rhs_adv)) DEALLOCATE(self % rhs_adv) + IF (ALLOCATED(self % rhs_phys)) DEALLOCATE(self % rhs_phys) + IF (ALLOCATED(self % rhs_lbc)) DEALLOCATE(self % rhs_lbc) + IF (ALLOCATED(self % adv_inc_prev)) DEALLOCATE(self % adv_inc_prev) + IF (ALLOCATED(self % mr_n)) DEALLOCATE(self % mr_n) + IF (ALLOCATED(self % mr_after_slow)) DEALLOCATE(self % mr_after_slow) + IF (ALLOCATED(self % mr_after_adv)) DEALLOCATE(self % mr_after_adv) + CALL self % dtheta % field_final + CALL self % du % field_final + CALL self % total_dry_flux % field_final + RETURN + END SUBROUTINE semi_implicit_alg_final + SUBROUTINE conditional_collection_copy(generic_fields_copied, generic_fields_to_copy, field_list) + USE semi_implicit_timestep_alg_mod_psy, ONLY: invoke_12 + USE field_collection_mod, ONLY: field_collection_type + USE field_collection_iterator_mod, ONLY: field_collection_iterator_type + IMPLICIT NONE + TYPE(field_collection_type), INTENT(OUT) :: generic_fields_copied + TYPE(field_collection_type), INTENT(IN) :: generic_fields_to_copy + TYPE(field_collection_type), INTENT(IN) :: field_list + TYPE(field_collection_iterator_type) :: iterator + CLASS(field_parent_type), POINTER :: abstract_field_ptr + TYPE(field_type), POINTER :: single_generic_field + TYPE(field_type) :: copied_generic_field + LOGICAL(KIND = l_def) :: l_copy + NULLIFY(abstract_field_ptr) + NULLIFY(single_generic_field) + CALL generic_fields_copied % initialise(name = 'fields_copied') + IF (generic_fields_to_copy % get_length() > 0) THEN + CALL iterator % initialise(generic_fields_to_copy) + DO + IF (.NOT. iterator % has_next()) EXIT + abstract_field_ptr => iterator % next() + SELECT TYPE(abstract_field_ptr) + TYPE IS (field_type) + single_generic_field => abstract_field_ptr + END SELECT + l_copy = field_list % field_exists(single_generic_field % get_name()) + IF (l_copy) THEN + CALL single_generic_field % copy_field_properties(copied_generic_field) + CALL invoke_12(copied_generic_field, single_generic_field) + CALL generic_fields_copied % add_field(copied_generic_field) + END IF + END DO + END IF + END SUBROUTINE conditional_collection_copy +END MODULE semi_implicit_timestep_alg_mod \ No newline at end of file diff --git a/science/gungho/source/psy/semi_implicit_timestep_alg_mod_psy.f90 b/science/gungho/source/psy/semi_implicit_timestep_alg_mod_psy.f90 new file mode 100644 index 000000000..c3313a13a --- /dev/null +++ b/science/gungho/source/psy/semi_implicit_timestep_alg_mod_psy.f90 @@ -0,0 +1,1090 @@ +module semi_implicit_timestep_alg_mod_psy + use constants_mod + use field_mod, only : field_proxy_type, field_proxy_type_10=>field_proxy_type, field_proxy_type_11=>field_proxy_type, & +&field_proxy_type_12=>field_proxy_type, field_proxy_type_1=>field_proxy_type, field_proxy_type_2=>field_proxy_type, & +&field_proxy_type_3=>field_proxy_type, field_proxy_type_4=>field_proxy_type, field_proxy_type_5=>field_proxy_type, & +&field_proxy_type_6=>field_proxy_type, field_proxy_type_7=>field_proxy_type, field_proxy_type_8=>field_proxy_type, & +&field_proxy_type_9=>field_proxy_type, field_type + use r_tran_field_mod, only : r_tran_field_proxy_type, r_tran_field_type + use operator_mod, only : operator_proxy_type, operator_type + implicit none + public + + contains + subroutine invoke_0(self_dtheta_cld) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_dtheta_cld + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_dtheta_cld_data => null() + type(field_proxy_type) :: self_dtheta_cld_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + + ! Initialise field and/or operator proxies + self_dtheta_cld_proxy = self_dtheta_cld%get_proxy() + self_dtheta_cld_data => self_dtheta_cld_proxy%data + + ! Create a mesh object + mesh => self_dtheta_cld_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_dtheta_cld_proxy%vspace%get_last_dof_halo(1) + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_c (set a real-valued field to a real scalar value) + self_dtheta_cld_data(df) = 0.0_r_def + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_dtheta_cld_proxy%set_dirty() + call self_dtheta_cld_proxy%set_clean(1) + + end subroutine invoke_0 + subroutine invoke_copy_init_fields_to_state(self_state, u, self_state_1, theta, self_state_2, rho, self_state_3, exner, & +&self_total_dry_flux) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_state + type(field_type), intent(in) :: u + type(field_type), intent(in) :: self_state_1 + type(field_type), intent(in) :: theta + type(field_type), intent(in) :: self_state_2 + type(field_type), intent(in) :: rho + type(field_type), intent(in) :: self_state_3 + type(field_type), intent(in) :: exner + type(r_tran_field_type), intent(in) :: self_total_dry_flux + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_state_data => null() + real(kind=r_def), pointer, dimension(:) :: u_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_1_data => null() + real(kind=r_def), pointer, dimension(:) :: theta_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_2_data => null() + real(kind=r_def), pointer, dimension(:) :: rho_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_3_data => null() + real(kind=r_def), pointer, dimension(:) :: exner_data => null() + real(kind=r_tran), pointer, dimension(:) :: self_total_dry_flux_data => null() + type(field_proxy_type_1) :: self_state_proxy + type(field_proxy_type_1) :: u_proxy + type(field_proxy_type_1) :: self_state_1_proxy + type(field_proxy_type_1) :: theta_proxy + type(field_proxy_type_1) :: self_state_2_proxy + type(field_proxy_type_1) :: rho_proxy + type(field_proxy_type_1) :: self_state_3_proxy + type(field_proxy_type_1) :: exner_proxy + type(r_tran_field_proxy_type) :: self_total_dry_flux_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + integer(kind=i_def) :: loop2_start + integer(kind=i_def) :: loop2_stop + integer(kind=i_def) :: loop3_start + integer(kind=i_def) :: loop3_stop + integer(kind=i_def) :: loop4_start + integer(kind=i_def) :: loop4_stop + + ! Initialise field and/or operator proxies + self_state_proxy = self_state%get_proxy() + self_state_data => self_state_proxy%data + u_proxy = u%get_proxy() + u_data => u_proxy%data + self_state_1_proxy = self_state_1%get_proxy() + self_state_1_data => self_state_1_proxy%data + theta_proxy = theta%get_proxy() + theta_data => theta_proxy%data + self_state_2_proxy = self_state_2%get_proxy() + self_state_2_data => self_state_2_proxy%data + rho_proxy = rho%get_proxy() + rho_data => rho_proxy%data + self_state_3_proxy = self_state_3%get_proxy() + self_state_3_data => self_state_3_proxy%data + exner_proxy = exner%get_proxy() + exner_data => exner_proxy%data + self_total_dry_flux_proxy = self_total_dry_flux%get_proxy() + self_total_dry_flux_data => self_total_dry_flux_proxy%data + + ! Create a mesh object + mesh => self_state_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_state_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = self_state_1_proxy%vspace%get_last_dof_annexed() + loop2_start = 1 + loop2_stop = self_state_2_proxy%vspace%get_last_dof_annexed() + loop3_start = 1 + loop3_stop = self_state_3_proxy%vspace%get_last_dof_annexed() + loop4_start = 1 + loop4_stop = self_total_dry_flux_proxy%vspace%get_last_dof_halo(1) + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_state_data(df) = u_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_state_1_data(df) = theta_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_1_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop2_start, loop2_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_state_2_data(df) = rho_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_2_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop3_start, loop3_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_state_3_data(df) = exner_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_3_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop4_start, loop4_stop, 1 + ! Built-in: setval_c (set a real-valued field to a real scalar value) + self_total_dry_flux_data(df) = 0.0_r_tran + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_total_dry_flux_proxy%set_dirty() + call self_total_dry_flux_proxy%set_clean(1) + + end subroutine invoke_copy_init_fields_to_state + subroutine invoke_2(theta_ref, self_state, rho_ref, self_state_1, exner_ref, self_state_2) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: theta_ref + type(field_type), intent(in) :: self_state + type(field_type), intent(in) :: rho_ref + type(field_type), intent(in) :: self_state_1 + type(field_type), intent(in) :: exner_ref + type(field_type), intent(in) :: self_state_2 + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: theta_ref_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_data => null() + real(kind=r_def), pointer, dimension(:) :: rho_ref_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_1_data => null() + real(kind=r_def), pointer, dimension(:) :: exner_ref_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_2_data => null() + type(field_proxy_type_2) :: theta_ref_proxy + type(field_proxy_type_2) :: self_state_proxy + type(field_proxy_type_2) :: rho_ref_proxy + type(field_proxy_type_2) :: self_state_1_proxy + type(field_proxy_type_2) :: exner_ref_proxy + type(field_proxy_type_2) :: self_state_2_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + integer(kind=i_def) :: loop2_start + integer(kind=i_def) :: loop2_stop + + ! Initialise field and/or operator proxies + theta_ref_proxy = theta_ref%get_proxy() + theta_ref_data => theta_ref_proxy%data + self_state_proxy = self_state%get_proxy() + self_state_data => self_state_proxy%data + rho_ref_proxy = rho_ref%get_proxy() + rho_ref_data => rho_ref_proxy%data + self_state_1_proxy = self_state_1%get_proxy() + self_state_1_data => self_state_1_proxy%data + exner_ref_proxy = exner_ref%get_proxy() + exner_ref_data => exner_ref_proxy%data + self_state_2_proxy = self_state_2%get_proxy() + self_state_2_data => self_state_2_proxy%data + + ! Create a mesh object + mesh => theta_ref_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = theta_ref_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = rho_ref_proxy%vspace%get_last_dof_annexed() + loop2_start = 1 + loop2_stop = exner_ref_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + theta_ref_data(df) = self_state_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call theta_ref_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + rho_ref_data(df) = self_state_1_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call rho_ref_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop2_start, loop2_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + exner_ref_data(df) = self_state_2_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call exner_ref_proxy%set_dirty() + + end subroutine invoke_2 + subroutine invoke_update_from_slow_physics(self_state_after_slow, self_dtheta, self_state_after_slow_1, self_du) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_state_after_slow + type(field_type), intent(in) :: self_dtheta + type(field_type), intent(in) :: self_state_after_slow_1 + type(field_type), intent(in) :: self_du + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_state_after_slow_data => null() + real(kind=r_def), pointer, dimension(:) :: self_dtheta_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_after_slow_1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_du_data => null() + type(field_proxy_type_3) :: self_state_after_slow_proxy + type(field_proxy_type_3) :: self_dtheta_proxy + type(field_proxy_type_3) :: self_state_after_slow_1_proxy + type(field_proxy_type_3) :: self_du_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + + ! Initialise field and/or operator proxies + self_state_after_slow_proxy = self_state_after_slow%get_proxy() + self_state_after_slow_data => self_state_after_slow_proxy%data + self_dtheta_proxy = self_dtheta%get_proxy() + self_dtheta_data => self_dtheta_proxy%data + self_state_after_slow_1_proxy = self_state_after_slow_1%get_proxy() + self_state_after_slow_1_data => self_state_after_slow_1_proxy%data + self_du_proxy = self_du%get_proxy() + self_du_data => self_du_proxy%data + + ! Create a mesh object + mesh => self_state_after_slow_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_state_after_slow_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = self_state_after_slow_1_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: inc_X_plus_Y (increment a real-valued field) + self_state_after_slow_data(df) = self_state_after_slow_data(df) + self_dtheta_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_after_slow_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: inc_X_plus_Y (increment a real-valued field) + self_state_after_slow_1_data(df) = self_state_after_slow_1_data(df) + self_du_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_after_slow_1_proxy%set_dirty() + + end subroutine invoke_update_from_slow_physics + subroutine invoke_4(self_wind_prev, self_state) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_wind_prev + type(field_type), intent(in) :: self_state + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_wind_prev_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_data => null() + type(field_proxy_type_4) :: self_wind_prev_proxy + type(field_proxy_type_4) :: self_state_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + + ! Initialise field and/or operator proxies + self_wind_prev_proxy = self_wind_prev%get_proxy() + self_wind_prev_data => self_wind_prev_proxy%data + self_state_proxy = self_state%get_proxy() + self_state_data => self_state_proxy%data + + ! Create a mesh object + mesh => self_wind_prev_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_wind_prev_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_wind_prev_data(df) = self_state_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_wind_prev_proxy%set_dirty() + + end subroutine invoke_4 + subroutine invoke_5(self_advected_state, self_rhs_adv, self_advected_state_1, self_rhs_adv_1, self_advected_state_2, self_du) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_advected_state + type(field_type), intent(in) :: self_rhs_adv + type(field_type), intent(in) :: self_advected_state_1 + type(field_type), intent(in) :: self_rhs_adv_1 + type(field_type), intent(in) :: self_advected_state_2 + type(field_type), intent(in) :: self_du + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_advected_state_data => null() + real(kind=r_def), pointer, dimension(:) :: self_rhs_adv_data => null() + real(kind=r_def), pointer, dimension(:) :: self_advected_state_1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_rhs_adv_1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_advected_state_2_data => null() + real(kind=r_def), pointer, dimension(:) :: self_du_data => null() + type(field_proxy_type_5) :: self_advected_state_proxy + type(field_proxy_type_5) :: self_rhs_adv_proxy + type(field_proxy_type_5) :: self_advected_state_1_proxy + type(field_proxy_type_5) :: self_rhs_adv_1_proxy + type(field_proxy_type_5) :: self_advected_state_2_proxy + type(field_proxy_type_5) :: self_du_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + integer(kind=i_def) :: loop2_start + integer(kind=i_def) :: loop2_stop + + ! Initialise field and/or operator proxies + self_advected_state_proxy = self_advected_state%get_proxy() + self_advected_state_data => self_advected_state_proxy%data + self_rhs_adv_proxy = self_rhs_adv%get_proxy() + self_rhs_adv_data => self_rhs_adv_proxy%data + self_advected_state_1_proxy = self_advected_state_1%get_proxy() + self_advected_state_1_data => self_advected_state_1_proxy%data + self_rhs_adv_1_proxy = self_rhs_adv_1%get_proxy() + self_rhs_adv_1_data => self_rhs_adv_1_proxy%data + self_advected_state_2_proxy = self_advected_state_2%get_proxy() + self_advected_state_2_data => self_advected_state_2_proxy%data + self_du_proxy = self_du%get_proxy() + self_du_data => self_du_proxy%data + + ! Create a mesh object + mesh => self_advected_state_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_advected_state_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = self_advected_state_1_proxy%vspace%get_last_dof_annexed() + loop2_start = 1 + loop2_stop = self_advected_state_2_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: inc_X_plus_Y (increment a real-valued field) + self_advected_state_data(df) = self_advected_state_data(df) + self_rhs_adv_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_advected_state_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: inc_X_plus_Y (increment a real-valued field) + self_advected_state_1_data(df) = self_advected_state_1_data(df) + self_rhs_adv_1_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_advected_state_1_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop2_start, loop2_stop, 1 + ! Built-in: inc_X_plus_Y (increment a real-valued field) + self_advected_state_2_data(df) = self_advected_state_2_data(df) + self_du_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_advected_state_2_proxy%set_dirty() + + end subroutine invoke_5 + subroutine invoke_6(self_wind_prev, self_state) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_wind_prev + type(field_type), intent(in) :: self_state + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_wind_prev_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_data => null() + type(field_proxy_type_6) :: self_wind_prev_proxy + type(field_proxy_type_6) :: self_state_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + + ! Initialise field and/or operator proxies + self_wind_prev_proxy = self_wind_prev%get_proxy() + self_wind_prev_data => self_wind_prev_proxy%data + self_state_proxy = self_state%get_proxy() + self_state_data => self_state_proxy%data + + ! Create a mesh object + mesh => self_wind_prev_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_wind_prev_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_wind_prev_data(df) = self_state_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_wind_prev_proxy%set_dirty() + + end subroutine invoke_6 + subroutine invoke_7(self_dtheta, self_rhs_adv, mm_wt) + use mesh_mod, only : mesh_type + use dg_inc_matrix_vector_kernel_mod, only : dg_inc_matrix_vector_code + type(field_type), intent(in) :: self_dtheta + type(field_type), intent(in) :: self_rhs_adv + type(operator_type), intent(in) :: mm_wt + integer(kind=i_def) :: df + integer(kind=i_def) :: cell + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_dtheta_data => null() + real(kind=r_def), pointer, dimension(:) :: self_rhs_adv_data => null() + real(kind=r_def), pointer, dimension(:,:,:) :: mm_wt_local_stencil => null() + integer(kind=i_def) :: nlayers_self_rhs_adv + integer(kind=i_def) :: ndf_aspc1_self_dtheta + integer(kind=i_def) :: undf_aspc1_self_dtheta + integer(kind=i_def) :: ndf_aspc1_self_rhs_adv + integer(kind=i_def) :: undf_aspc1_self_rhs_adv + integer(kind=i_def) :: ndf_adspc1_self_rhs_adv + integer(kind=i_def) :: undf_adspc1_self_rhs_adv + integer(kind=i_def), pointer :: map_adspc1_self_rhs_adv(:,:) => null() + integer(kind=i_def), pointer :: map_aspc1_self_dtheta(:,:) => null() + type(field_proxy_type_7) :: self_dtheta_proxy + type(field_proxy_type_7) :: self_rhs_adv_proxy + type(operator_proxy_type) :: mm_wt_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + integer(kind=i_def) :: loop2_start + integer(kind=i_def) :: loop2_stop + + ! Initialise field and/or operator proxies + self_dtheta_proxy = self_dtheta%get_proxy() + self_dtheta_data => self_dtheta_proxy%data + self_rhs_adv_proxy = self_rhs_adv%get_proxy() + self_rhs_adv_data => self_rhs_adv_proxy%data + mm_wt_proxy = mm_wt%get_proxy() + mm_wt_local_stencil => mm_wt_proxy%local_stencil + + ! Initialise number of layers + nlayers_self_rhs_adv = self_rhs_adv_proxy%vspace%get_nlayers() + + ! Create a mesh object + mesh => self_dtheta_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Look-up dofmaps for each function space + map_adspc1_self_rhs_adv => self_rhs_adv_proxy%vspace%get_whole_dofmap() + map_aspc1_self_dtheta => self_dtheta_proxy%vspace%get_whole_dofmap() + + ! Initialise number of DoFs for aspc1_self_dtheta + ndf_aspc1_self_dtheta = self_dtheta_proxy%vspace%get_ndf() + undf_aspc1_self_dtheta = self_dtheta_proxy%vspace%get_undf() + + ! Initialise number of DoFs for aspc1_self_rhs_adv + ndf_aspc1_self_rhs_adv = self_rhs_adv_proxy%vspace%get_ndf() + undf_aspc1_self_rhs_adv = self_rhs_adv_proxy%vspace%get_undf() + + ! Initialise number of DoFs for adspc1_self_rhs_adv + ndf_adspc1_self_rhs_adv = self_rhs_adv_proxy%vspace%get_ndf() + undf_adspc1_self_rhs_adv = self_rhs_adv_proxy%vspace%get_undf() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_dtheta_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = self_rhs_adv_proxy%vspace%get_last_dof_halo(1) + loop2_start = 1 + loop2_stop = mesh%get_last_edge_cell() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + self_dtheta_data(df) = self_rhs_adv_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_dtheta_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: setval_c (set a real-valued field to a real scalar value) + self_rhs_adv_data(df) = 0.0_r_def + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_adv_proxy%set_dirty() + call self_rhs_adv_proxy%set_clean(1) + !$omp parallel default(shared) private(cell) + !$omp do schedule(static) + do cell = loop2_start, loop2_stop, 1 + call dg_inc_matrix_vector_code(cell, nlayers_self_rhs_adv, self_rhs_adv_data, self_dtheta_data, mm_wt_proxy%ncell_3d, & +&mm_wt_local_stencil, ndf_adspc1_self_rhs_adv, undf_adspc1_self_rhs_adv, map_adspc1_self_rhs_adv(:,cell), ndf_aspc1_self_dtheta, & +&undf_aspc1_self_dtheta, map_aspc1_self_dtheta(:,cell)) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_adv_proxy%set_dirty() + + end subroutine invoke_7 + subroutine invoke_update_rhs_phys_from_fast_physics(self_rhs_phys, self_dtheta, mm_wt, self_rhs_phys_1, self_du, mm_vel) + use mesh_mod, only : mesh_type + use dg_inc_matrix_vector_kernel_mod, only : dg_inc_matrix_vector_code + use matrix_vector_kernel_mod, only : matrix_vector_code + use sci_enforce_bc_kernel_mod, only : enforce_bc_code + type(field_type), intent(in) :: self_rhs_phys + type(field_type), intent(in) :: self_dtheta + type(operator_type), intent(in) :: mm_wt + type(field_type), intent(in) :: self_rhs_phys_1 + type(field_type), intent(in) :: self_du + type(operator_type), intent(in) :: mm_vel + integer(kind=i_def) :: cell + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_rhs_phys_data => null() + real(kind=r_def), pointer, dimension(:) :: self_dtheta_data => null() + real(kind=r_def), pointer, dimension(:) :: self_rhs_phys_1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_du_data => null() + real(kind=r_def), pointer, dimension(:,:,:) :: mm_wt_local_stencil => null() + real(kind=r_def), pointer, dimension(:,:,:) :: mm_vel_local_stencil => null() + integer(kind=i_def) :: nlayers_self_rhs_phys + integer(kind=i_def) :: nlayers_self_rhs_phys_1 + integer(kind=i_def) :: colour + integer(kind=i_def), pointer :: cmap(:,:) + integer(kind=i_def) :: ncolour + integer(kind=i_def), allocatable, dimension(:,:) :: last_halo_cell_all_colours + integer(kind=i_def) :: ndf_adspc1_self_rhs_phys + integer(kind=i_def) :: undf_adspc1_self_rhs_phys + integer(kind=i_def) :: ndf_aspc1_self_dtheta + integer(kind=i_def) :: undf_aspc1_self_dtheta + integer(kind=i_def) :: ndf_aspc1_self_rhs_phys_1 + integer(kind=i_def) :: undf_aspc1_self_rhs_phys_1 + integer(kind=i_def) :: ndf_aspc2_self_du + integer(kind=i_def) :: undf_aspc2_self_du + integer(kind=i_def), pointer :: map_adspc1_self_rhs_phys(:,:) => null() + integer(kind=i_def), pointer :: map_aspc1_self_dtheta(:,:) => null() + integer(kind=i_def), pointer :: map_aspc1_self_rhs_phys_1(:,:) => null() + integer(kind=i_def), pointer :: map_aspc2_self_du(:,:) => null() + integer(kind=i_def), pointer :: boundary_dofs_self_rhs_phys_1(:,:) => null() + type(field_proxy_type_8) :: self_rhs_phys_proxy + type(field_proxy_type_8) :: self_dtheta_proxy + type(field_proxy_type_8) :: self_rhs_phys_1_proxy + type(field_proxy_type_8) :: self_du_proxy + type(operator_proxy_type) :: mm_wt_proxy + type(operator_proxy_type) :: mm_vel_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + integer(kind=i_def) :: loop2_start + integer(kind=i_def) :: loop3_start + integer(kind=i_def) :: loop3_stop + integer(kind=i_def) :: loop4_start + + ! Initialise field and/or operator proxies + self_rhs_phys_proxy = self_rhs_phys%get_proxy() + self_rhs_phys_data => self_rhs_phys_proxy%data + self_dtheta_proxy = self_dtheta%get_proxy() + self_dtheta_data => self_dtheta_proxy%data + mm_wt_proxy = mm_wt%get_proxy() + mm_wt_local_stencil => mm_wt_proxy%local_stencil + self_rhs_phys_1_proxy = self_rhs_phys_1%get_proxy() + self_rhs_phys_1_data => self_rhs_phys_1_proxy%data + self_du_proxy = self_du%get_proxy() + self_du_data => self_du_proxy%data + mm_vel_proxy = mm_vel%get_proxy() + mm_vel_local_stencil => mm_vel_proxy%local_stencil + + ! Initialise number of layers + nlayers_self_rhs_phys = self_rhs_phys_proxy%vspace%get_nlayers() + nlayers_self_rhs_phys_1 = self_rhs_phys_1_proxy%vspace%get_nlayers() + + ! Create a mesh object + mesh => self_rhs_phys_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Get the colourmap + ncolour = mesh%get_ncolours() + cmap => mesh%get_colour_map() + + ! Look-up dofmaps for each function space + map_adspc1_self_rhs_phys => self_rhs_phys_proxy%vspace%get_whole_dofmap() + map_aspc1_self_dtheta => self_dtheta_proxy%vspace%get_whole_dofmap() + map_aspc1_self_rhs_phys_1 => self_rhs_phys_1_proxy%vspace%get_whole_dofmap() + map_aspc2_self_du => self_du_proxy%vspace%get_whole_dofmap() + boundary_dofs_self_rhs_phys_1 => self_rhs_phys_1_proxy%vspace%get_boundary_dofs() + + ! Initialise number of DoFs for adspc1_self_rhs_phys + ndf_adspc1_self_rhs_phys = self_rhs_phys_proxy%vspace%get_ndf() + undf_adspc1_self_rhs_phys = self_rhs_phys_proxy%vspace%get_undf() + + ! Initialise number of DoFs for aspc1_self_dtheta + ndf_aspc1_self_dtheta = self_dtheta_proxy%vspace%get_ndf() + undf_aspc1_self_dtheta = self_dtheta_proxy%vspace%get_undf() + + ! Initialise number of DoFs for aspc1_self_rhs_phys_1 + ndf_aspc1_self_rhs_phys_1 = self_rhs_phys_1_proxy%vspace%get_ndf() + undf_aspc1_self_rhs_phys_1 = self_rhs_phys_1_proxy%vspace%get_undf() + + ! Initialise number of DoFs for aspc2_self_du + ndf_aspc2_self_du = self_du_proxy%vspace%get_ndf() + undf_aspc2_self_du = self_du_proxy%vspace%get_undf() + last_halo_cell_all_colours = mesh%get_last_halo_cell_all_colours() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = mesh%get_last_edge_cell() + loop1_start = 1 + loop1_stop = ncolour + loop2_start = 1 + loop3_start = 1 + loop3_stop = ncolour + loop4_start = 1 + + ! Call kernels and communication routines + !$omp parallel default(shared) private(cell) + !$omp do schedule(static) + do cell = loop0_start, loop0_stop, 1 + call dg_inc_matrix_vector_code(cell, nlayers_self_rhs_phys, self_rhs_phys_data, self_dtheta_data, mm_wt_proxy%ncell_3d, & +&mm_wt_local_stencil, ndf_adspc1_self_rhs_phys, undf_adspc1_self_rhs_phys, map_adspc1_self_rhs_phys(:,cell), & +&ndf_aspc1_self_dtheta, undf_aspc1_self_dtheta, map_aspc1_self_dtheta(:,cell)) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_phys_proxy%set_dirty() + if (self_du_proxy%is_dirty(depth=1)) then + call self_du_proxy%halo_exchange(depth=1) + end if + do colour = loop1_start, loop1_stop, 1 + !$omp parallel default(shared) private(cell) + !$omp do schedule(static) + do cell = loop2_start, last_halo_cell_all_colours(colour,1), 1 + call matrix_vector_code(cmap(colour,cell), nlayers_self_rhs_phys_1, self_rhs_phys_1_data, self_du_data, & +&mm_vel_proxy%ncell_3d, mm_vel_local_stencil, ndf_aspc1_self_rhs_phys_1, undf_aspc1_self_rhs_phys_1, & +&map_aspc1_self_rhs_phys_1(:,cmap(colour,cell)), ndf_aspc2_self_du, undf_aspc2_self_du, map_aspc2_self_du(:,cmap(colour,cell))) + enddo + !$omp end do + !$omp end parallel + enddo + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_phys_1_proxy%set_dirty() + do colour = loop3_start, loop3_stop, 1 + !$omp parallel default(shared) private(cell) + !$omp do schedule(static) + do cell = loop4_start, last_halo_cell_all_colours(colour,1), 1 + call enforce_bc_code(nlayers_self_rhs_phys_1, self_rhs_phys_1_data, ndf_aspc1_self_rhs_phys_1, undf_aspc1_self_rhs_phys_1, & +&map_aspc1_self_rhs_phys_1(:,cmap(colour,cell)), boundary_dofs_self_rhs_phys_1) + enddo + !$omp end do + !$omp end parallel + enddo + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_phys_1_proxy%set_dirty() + + end subroutine invoke_update_rhs_phys_from_fast_physics + subroutine invoke_9(self_rhs_np1, self_rhs_np1_1) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_rhs_np1 + type(field_type), intent(in) :: self_rhs_np1_1 + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_rhs_np1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_rhs_np1_1_data => null() + type(field_proxy_type_9) :: self_rhs_np1_proxy + type(field_proxy_type_9) :: self_rhs_np1_1_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + + ! Initialise field and/or operator proxies + self_rhs_np1_proxy = self_rhs_np1%get_proxy() + self_rhs_np1_data => self_rhs_np1_proxy%data + self_rhs_np1_1_proxy = self_rhs_np1_1%get_proxy() + self_rhs_np1_1_data => self_rhs_np1_1_proxy%data + + ! Create a mesh object + mesh => self_rhs_np1_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_rhs_np1_proxy%vspace%get_last_dof_halo(1) + loop1_start = 1 + loop1_stop = self_rhs_np1_1_proxy%vspace%get_last_dof_halo(1) + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_c (set a real-valued field to a real scalar value) + self_rhs_np1_data(df) = 0.0_r_def + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_np1_proxy%set_dirty() + call self_rhs_np1_proxy%set_clean(1) + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: setval_c (set a real-valued field to a real scalar value) + self_rhs_np1_1_data(df) = 0.0_r_def + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_rhs_np1_1_proxy%set_dirty() + call self_rhs_np1_1_proxy%set_clean(1) + + end subroutine invoke_9 + subroutine invoke_10(self_state, self_state_after_slow, self_dtheta, self_state_1, self_state_after_slow_1, self_du) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: self_state + type(field_type), intent(in) :: self_state_after_slow + type(field_type), intent(in) :: self_dtheta + type(field_type), intent(in) :: self_state_1 + type(field_type), intent(in) :: self_state_after_slow_1 + type(field_type), intent(in) :: self_du + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: self_state_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_after_slow_data => null() + real(kind=r_def), pointer, dimension(:) :: self_dtheta_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_after_slow_1_data => null() + real(kind=r_def), pointer, dimension(:) :: self_du_data => null() + type(field_proxy_type_10) :: self_state_proxy + type(field_proxy_type_10) :: self_state_after_slow_proxy + type(field_proxy_type_10) :: self_dtheta_proxy + type(field_proxy_type_10) :: self_state_1_proxy + type(field_proxy_type_10) :: self_state_after_slow_1_proxy + type(field_proxy_type_10) :: self_du_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + + ! Initialise field and/or operator proxies + self_state_proxy = self_state%get_proxy() + self_state_data => self_state_proxy%data + self_state_after_slow_proxy = self_state_after_slow%get_proxy() + self_state_after_slow_data => self_state_after_slow_proxy%data + self_dtheta_proxy = self_dtheta%get_proxy() + self_dtheta_data => self_dtheta_proxy%data + self_state_1_proxy = self_state_1%get_proxy() + self_state_1_data => self_state_1_proxy%data + self_state_after_slow_1_proxy = self_state_after_slow_1%get_proxy() + self_state_after_slow_1_data => self_state_after_slow_1_proxy%data + self_du_proxy = self_du%get_proxy() + self_du_data => self_du_proxy%data + + ! Create a mesh object + mesh => self_state_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = self_state_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = self_state_1_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: X_plus_Y (add real-valued fields) + self_state_data(df) = self_state_after_slow_data(df) + self_dtheta_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: X_plus_Y (add real-valued fields) + self_state_1_data(df) = self_state_after_slow_1_data(df) + self_du_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call self_state_1_proxy%set_dirty() + + end subroutine invoke_10 + subroutine invoke_11(u, self_state, theta, self_state_1, rho, self_state_2, exner, self_state_3) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: u + type(field_type), intent(in) :: self_state + type(field_type), intent(in) :: theta + type(field_type), intent(in) :: self_state_1 + type(field_type), intent(in) :: rho + type(field_type), intent(in) :: self_state_2 + type(field_type), intent(in) :: exner + type(field_type), intent(in) :: self_state_3 + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: u_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_data => null() + real(kind=r_def), pointer, dimension(:) :: theta_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_1_data => null() + real(kind=r_def), pointer, dimension(:) :: rho_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_2_data => null() + real(kind=r_def), pointer, dimension(:) :: exner_data => null() + real(kind=r_def), pointer, dimension(:) :: self_state_3_data => null() + type(field_proxy_type_11) :: u_proxy + type(field_proxy_type_11) :: self_state_proxy + type(field_proxy_type_11) :: theta_proxy + type(field_proxy_type_11) :: self_state_1_proxy + type(field_proxy_type_11) :: rho_proxy + type(field_proxy_type_11) :: self_state_2_proxy + type(field_proxy_type_11) :: exner_proxy + type(field_proxy_type_11) :: self_state_3_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + integer(kind=i_def) :: loop1_start + integer(kind=i_def) :: loop1_stop + integer(kind=i_def) :: loop2_start + integer(kind=i_def) :: loop2_stop + integer(kind=i_def) :: loop3_start + integer(kind=i_def) :: loop3_stop + + ! Initialise field and/or operator proxies + u_proxy = u%get_proxy() + u_data => u_proxy%data + self_state_proxy = self_state%get_proxy() + self_state_data => self_state_proxy%data + theta_proxy = theta%get_proxy() + theta_data => theta_proxy%data + self_state_1_proxy = self_state_1%get_proxy() + self_state_1_data => self_state_1_proxy%data + rho_proxy = rho%get_proxy() + rho_data => rho_proxy%data + self_state_2_proxy = self_state_2%get_proxy() + self_state_2_data => self_state_2_proxy%data + exner_proxy = exner%get_proxy() + exner_data => exner_proxy%data + self_state_3_proxy = self_state_3%get_proxy() + self_state_3_data => self_state_3_proxy%data + + ! Create a mesh object + mesh => u_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = u_proxy%vspace%get_last_dof_annexed() + loop1_start = 1 + loop1_stop = theta_proxy%vspace%get_last_dof_annexed() + loop2_start = 1 + loop2_stop = rho_proxy%vspace%get_last_dof_annexed() + loop3_start = 1 + loop3_stop = exner_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + u_data(df) = self_state_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call u_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop1_start, loop1_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + theta_data(df) = self_state_1_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call theta_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop2_start, loop2_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + rho_data(df) = self_state_2_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call rho_proxy%set_dirty() + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop3_start, loop3_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + exner_data(df) = self_state_3_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call exner_proxy%set_dirty() + + end subroutine invoke_11 + subroutine invoke_12(copied_generic_field, single_generic_field) + use mesh_mod, only : mesh_type + type(field_type), intent(in) :: copied_generic_field + type(field_type), intent(in) :: single_generic_field + integer(kind=i_def) :: df + type(mesh_type), pointer :: mesh => null() + integer(kind=i_def) :: max_halo_depth_mesh + real(kind=r_def), pointer, dimension(:) :: copied_generic_field_data => null() + real(kind=r_def), pointer, dimension(:) :: single_generic_field_data => null() + type(field_proxy_type_12) :: copied_generic_field_proxy + type(field_proxy_type_12) :: single_generic_field_proxy + integer(kind=i_def) :: loop0_start + integer(kind=i_def) :: loop0_stop + + ! Initialise field and/or operator proxies + copied_generic_field_proxy = copied_generic_field%get_proxy() + copied_generic_field_data => copied_generic_field_proxy%data + single_generic_field_proxy = single_generic_field%get_proxy() + single_generic_field_data => single_generic_field_proxy%data + + ! Create a mesh object + mesh => copied_generic_field_proxy%vspace%get_mesh() + max_halo_depth_mesh = mesh%get_halo_depth() + + ! Set-up all of the loop bounds + loop0_start = 1 + loop0_stop = copied_generic_field_proxy%vspace%get_last_dof_annexed() + + ! Call kernels and communication routines + !$omp parallel default(shared) private(df) + !$omp do schedule(static) + do df = loop0_start, loop0_stop, 1 + ! Built-in: setval_X (set a real-valued field equal to another such field) + copied_generic_field_data(df) = single_generic_field_data(df) + enddo + !$omp end do + !$omp end parallel + + ! Set halos dirty/clean for fields modified in the above loop(s) + call copied_generic_field_proxy%set_dirty() + + end subroutine invoke_12 + +end module semi_implicit_timestep_alg_mod_psy