Merge branch 'main' into add-fortran-linter

Author: james-bruten-mo

applications/lfric_atm/metadata/axis_def_main.xml

@@ -33,10 +33,11 @@
   <axis id="sw_bands_radiation_levels" name="sw_bands_radiation_levels" />
   <axis id="lw_bands_radiation_levels" name="lw_bands_radiation_levels" />
   <axis id="photolysis_pathways" name="photolysis_pathways" />
-  <axis id="photol_species" name="photol_species" />  
+  <axis id="photol_species" name="photol_species" />
   <axis id="one_time_axis" name="one_time_axis" n_glo="1" />
 
   <axis id="random_seed_size" name="random_seed_size" />
+  <axis id="stph_spectral_dimensions" name="stph_spectral_dimensions" />
   <axis id="scalar_axis" name="scalar" n_glo="1" />
   <axis id="aod_wavel" name="aod_wavel" />
 

applications/lfric_atm/metadata/grid_def_main.xml

@@ -40,7 +40,7 @@
     <domain domain_ref="face"/>
     <axis axis_ref="vert_axis_full_levels"/>
     <axis axis_ref="photol_species"/>
-  </grid>  
+  </grid>
   <grid id="3D_monthly_grid">
     <domain domain_ref="face"/>
     <axis axis_ref="vert_axis_full_levels"/>
@@ -118,6 +118,9 @@
   <grid id="random_seed_grid">
     <axis axis_ref="random_seed_size"/>
   </grid>
+  <grid id="stochastic_physics_grid">
+    <axis axis_ref="stph_spectral_dimensions"/>
+  </grid>
   <!-- Zoom filters -->
   <grid id="full_level_0_face_zoom">
     <domain domain_ref="face"/>

applications/lfric_atm/metadata/lfric_dictionary.xml

@@ -351,10 +351,16 @@
   <field id="aer_lw_absorption" name="aer_lw_absorption" long_name="aer_lw_absorption" unit="1" grid_ref="full_level_face_grid_lw_bands_aero_modes" />
   <field id="aer_lw_scattering" name="aer_lw_scattering" long_name="aer_lw_scattering" unit="1" grid_ref="full_level_face_grid_lw_bands_aero_modes" />
   <field id="aer_lw_asymmetry" name="aer_lw_asymmetry" long_name="aer_lw_asymmetry" unit="1" grid_ref="full_level_face_grid_lw_bands_aero_modes" />
-  <!-- stochastic fields --> 
+  <!-- stochastic fields -->
   <field id="blpert_rand_fld" name="blpert_rand_fld" long_name="blpert_random_number" unit="1" domain_ref="face" />
   <field id="blpert_flag" name="blpert_flag" long_name="blpert_perturbed_flag" unit="1" domain_ref="face" />
   <field id="random_seed" name="random_seed" long_name="random_seed_for_stochastic_physics" unit="1" grid_ref="random_seed_grid" />
+  <field id="spt_spectral_coeffc" name="spt_spectral_coeffc" long_name="spt_real_spectral_coefficients_for_stochastic_physics" unit="1" grid_ref="stochastic_physics_grid" />
+  <field id="spt_spectral_coeffs" name="spt_spectral_coeffs" long_name="spt_imaginary_spectral_coefficients_for_stochastic_physics" unit="1" grid_ref="stochastic_physics_grid" />
+  <field id="spt_power_law" name="spt_power_law" long_name="spt_power_law_for_stochastic_physics" unit="1" grid_ref="stochastic_physics_grid" />
+  <field id="skeb_spectral_coeffc" name="skeb_spectral_coeffc" long_name="skeb_real_spectral_coefficients_for_stochastic_physics" unit="1" grid_ref="stochastic_physics_grid" />
+  <field id="skeb_spectral_coeffs" name="skeb_spectral_coeffs" long_name="skeb_imaginary_spectral_coefficients_for_stochastic_physics" unit="1" grid_ref="stochastic_physics_grid" />
+  <field id="skeb_power_law" name="skeb_power_law" long_name="skeb_power_law_for_stochastic_physics" unit="1" grid_ref="stochastic_physics_grid" />
   <!-- snow fields -->
   <field id="tile_snow_mass" name="tile_snow_mass" long_name="tile_snow_mass" unit="kg m-2" domain_ref="face" axis_ref="surface_tiles" />
   <field id="tile_snow_rgrain" name="tile_snow_rgrain" long_name="tile_snow_grain_size" unit="microns" domain_ref="face" axis_ref="surface_tiles" />

interfaces/physics_schemes_interface/source/algorithm/skeb_main_alg_mod.x90

@@ -21,6 +21,9 @@ module skeb_main_alg_mod
     ! use collections
     use function_space_collection_mod, only: function_space_collection
     use mesh_collection_mod,           only: mesh_collection
+    ! Physical constants
+    use planet_constants_mod, only: planet_radius
+
     ! xios output and timers
     use io_config_mod, only: write_diag, use_xios_io
     use timing_mod,    only: start_timing, stop_timing, tik, LPROF
@@ -38,8 +41,6 @@ module skeb_main_alg_mod
 
     private
 
-    ! Logical controlling whether spectral coeffs need calculating
-    logical(kind=l_def), save :: initialize_skeb_spectral_coeffs = .true.
     logical(kind=l_def ) :: du_rot_skeb_flag
     logical(kind=l_def ) :: dv_rot_skeb_flag
     logical(kind=l_def ) :: du_div_skeb_flag
@@ -56,7 +57,7 @@ module skeb_main_alg_mod
     real(kind=r_def), parameter :: TWOP_P1 = 2.0_r_def*P + 1.0_r_def
     real(kind=r_def), parameter :: FOUR_ON_VAR = 48.0_r_def
 
-    public skeb_main_alg
+    public skeb_main_alg, skeb_init_power_law
 
  contains
   !>@brief Run the Stochastic Kinetic Energy Backscatter (SKEB)
@@ -87,13 +88,18 @@ module skeb_main_alg_mod
   !>         8) Add wind increments to the flow and diagnostics
   !>
   !>         See UMD81 for full scheme details
-  !>@param[in]         du_stph            Stochastic Physics increments for winds
-  !>@param[in]         rho                Density on W3
-  !>@param[in]         u                  prognostic winds in W2
-  !>@param[in]         convection_fields  Fields from convection scheme
-  !>@param[in]         clock              Model time information
-
-  subroutine skeb_main_alg(du_stph, rho, u, convection_fields, derived_fields, clock)
+  !>@param[in]         du_stph               Stochastic Physics increments for winds
+  !>@param[in]         rho                   Density on W3
+  !>@param[in]         u                     prognostic winds in W2
+  !>@param[in]         convection_fields     Fields from convection scheme
+  !>@param[in,out]     skeb_spectral_coeffc  Real Spectral coefficients
+  !>@param[in,out]     skeb_spectral_coeffs  Imaginary Spectral coefficients
+  !>@param[in,out]     skeb_power_law        Spectral power law
+  !>@param[in]         clock                 Model time information
+
+  subroutine skeb_main_alg(du_stph, rho, u, convection_fields, derived_fields, &
+                       skeb_spectral_coeffc, skeb_spectral_coeffs,             &
+                       skeb_power_law, clock)
 
   ! SKEB namelist options
   use stochastic_physics_config_mod,  only:                                         &
@@ -174,9 +180,6 @@ module skeb_main_alg_mod
   ! extract u and v from SKEB winds for diagnostics
   use physics_mappings_alg_mod,      only: map_physics_winds
 
-  ! Physical constants
-  use planet_constants_mod, only: planet_radius
-
   use sci_enforce_upper_bound_kernel_mod, only: enforce_upper_bound_kernel_type
   use enforce_crit_kernel_mod, only: enforce_crit_kernel_type
 
@@ -193,6 +196,11 @@ module skeb_main_alg_mod
   type( field_collection_type ), intent(in) :: convection_fields
   type( field_collection_type ), intent(in) :: derived_fields
 
+  ! Spectral coefficients
+  real(r_def), intent(inout) :: skeb_spectral_coeffc(:)
+  real(r_def), intent(inout) :: skeb_spectral_coeffs(:)
+  real(r_def), intent(inout) :: skeb_power_law(:)
+
   ! classes
   class(clock_type), intent(in)   :: clock
 
@@ -274,27 +282,22 @@ module skeb_main_alg_mod
   !w2 intermediate field for scaling by cos(lat)
   type( field_type ) :: field_w2
 
-  ! Spectral coefficients
-  real(kind=r_def), allocatable, save :: skeb_spectral_coeffc(:)
-  real(kind=r_def), allocatable, save :: skeb_spectral_coeffs(:)
-  real(kind=r_def), allocatable, save :: skeb_power_law(:)
-
   ! time decorrelation parameter
-  real(kind=r_def), save :: skeb_alpha
+  real(kind=r_def) :: skeb_alpha
 
   ! Convective resolution factor to modulate conv mask
   real(kind=r_def) :: convective_resolution_factor
 
   ! scaling total dissipation by br/total_backscatter
   real(kind=r_def) :: factor_psif, max_backscatter, crit_backscatter
   ! parameters for the power law and decorrelation time alpha
-  real(kind=r_def) :: gamma, dt, spl_coeff
+  real(kind=r_def) :: dt
 
   ! mesh_id
   integer(kind=i_def) :: mesh_id
 
   ! Iterators in for loops
-  integer(i_def) :: m, n, n_row, stencil_extent
+  integer(i_def) :: n, stencil_extent
 
   ! SKEB calculations are performed at cell centres (W3 points). To index the
   ! "eta_theta_levels" array correctly at these points, set offset to 1
@@ -325,66 +328,13 @@ module skeb_main_alg_mod
                setval_c(fp_skeb, 0.0_r_def) )
 
   dt = real(clock%get_seconds_per_step(), r_def)
+  skeb_alpha=1-exp(-dt/skeb_decorrelation_time)
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!  1)  Create Forcing Pattern  !!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  ! Initialize spectral coefficients for the forcing pattern
-  if (initialize_skeb_spectral_coeffs) then
-
-    !allocate spectral coefficients and power law
-    allocate(skeb_spectral_coeffc(stph_spectral_dim))
-    allocate(skeb_spectral_coeffs(stph_spectral_dim))
-    allocate(skeb_power_law(stph_spectral_dim))
-
-    ! set them to zero (invokes don't work for non-fields types)
-    skeb_spectral_coeffc = 0.0_r_def
-    skeb_spectral_coeffs = 0.0_r_def
-    skeb_power_law = 0.0_r_def
-
-    !!!!!! 1.a compute power law
-
-    ! compute alpha for temporal decorrelation
-    skeb_alpha=1-exp(-dt/skeb_decorrelation_time)
-
-    ! This power law is what Glenn Shutts observed in the CRM simulations
-    ! (see Berner et al., 2009: J. Atmos. Sci, pp 603-626)
-    gamma = 0.0_r_def
-    do n = stph_n_min, stph_n_max
-      gamma = gamma + (n+1)*(2*n+1)*n**(TWOP_P1)
-    end do
-    gamma = gamma/skeb_alpha
-
-    ! Below we calculate the power law which is given as
-    ! g(n) = spl_coeff * n^p
-    ! where
-    ! spl_coeff = SQRT([4 * Pi * a^2 * skeb_total_backscatter]/[dt * var * gamma])
-    ! var is the variance of random numbers [-0.5; 0.5] = 1/12
-    ! tested by 1000 cases of random arrays of size = 1e9
-    ! note: (4/var) is pre-calculated as FOUR_ON_VAR = 48
-
-    ! Set n_row as the summatory of n
-    n_row = 0
-    ! add up those scales below the minimum wavenumber to the
-    ! row iterator
-    do n = 1, stph_n_min-1
-      n_row = n_row + n
-    end do
-
-    ! Build power law
-    spl_coeff = planet_radius * sqrt(FOUR_ON_VAR*pi*skeb_total_backscatter/  &
-                  (dt*gamma))
-    do n = stph_n_min, stph_n_max
-      n_row = n_row + n
-      do m = 0, n
-        skeb_power_law(n_row + m) = spl_coeff * n**P
-      end do
-    end do
-    initialize_skeb_spectral_coeffs = .false.
-  end if
-
-  !!!!!! 1.b call stph_fp_main to create forcing pattern for SKEB
+  !!!!!! Create forcing pattern for SKEB
   call stph_fp_main_alg(skeb_level_bottom-1_i_def, skeb_level_top-1_i_def,     &
                         lev_offset, stph_n_min, stph_n_max, stph_spectral_dim, &
                         skeb_alpha, skeb_power_law,                            &
@@ -701,4 +651,64 @@ module skeb_main_alg_mod
   end if ! end if write_diags and use_xios
 
   end subroutine skeb_main_alg
+
+  ! Initialise the power law array
+  ! @param[out] skeb_power_law Power law array
+  ! @param[in]  clock          Model clock
+  subroutine skeb_init_power_law(skeb_power_law, clock)
+    use stochastic_physics_config_mod, only: skeb_decorrelation_time, &
+                                             skeb_total_backscatter,  &
+                                             stph_n_min, stph_n_max
+
+    real(r_def), intent(out) :: skeb_power_law(:)
+    class(clock_type), intent(in)   :: clock
+
+    ! time decorrelation parameter
+    real(kind=r_def) :: skeb_alpha
+    real(kind=r_def) :: gamma, dt, spl_coeff
+
+    ! Iterators in for loops
+    integer(i_def) :: m, n, n_row
+
+    skeb_power_law = 0.0_r_def
+
+    ! compute alpha for temporal decorrelation
+    dt = real(clock%get_seconds_per_step(), r_def)
+    skeb_alpha=1-exp(-dt/skeb_decorrelation_time)
+
+    ! This power law is what Glenn Shutts observed in the CRM simulations
+    ! (see Berner et al., 2009: J. Atmos. Sci, pp 603-626)
+    gamma = 0.0_r_def
+    do n = stph_n_min, stph_n_max
+      gamma = gamma + (n+1)*(2*n+1)*n**(TWOP_P1)
+    end do
+    gamma = gamma/skeb_alpha
+
+    ! Below we calculate the power law which is given as
+    ! g(n) = spl_coeff * n^p
+    ! where
+    ! spl_coeff = SQRT([4 * Pi * a^2 * skeb_total_backscatter]/[dt * var * gamma])
+    ! var is the variance of random numbers [-0.5; 0.5] = 1/12
+    ! tested by 1000 cases of random arrays of size = 1e9
+    ! note: (4/var) is pre-calculated as FOUR_ON_VAR = 48
+
+    ! Set n_row as the summatory of n
+    n_row = 0
+    ! add up those scales below the minimum wavenumber to the
+    ! row iterator
+    do n = 1, stph_n_min-1
+      n_row = n_row + n
+    end do
+
+    ! Build power law
+    spl_coeff = planet_radius * sqrt(FOUR_ON_VAR*pi*skeb_total_backscatter/  &
+                  (dt*gamma))
+    do n = stph_n_min, stph_n_max
+      n_row = n_row + n
+      do m = 0, n
+        skeb_power_law(n_row + m) = spl_coeff * n**P
+      end do
+    end do
+
+  end subroutine skeb_init_power_law
 end module skeb_main_alg_mod