PC2 response to horizontal mixing

interfaces/physics_schemes_interface/rose-meta/um-cloud/HEAD/rose-meta.conf

@@ -482,6 +482,26 @@ sort-key=Panel-A02b
 value-titles=Original,Smooth
 values='original','smooth'
 
+[namelist:cloud=pc2_turb_horiz]
+compulsory=true
+description=Do PC2 homogeneous forcing from horizontal turbulent mixing
+help=If this is false, the vertical component of turbulent mixing
+    =(the 1D boundary-layer scheme) is treated as a homogeneous forcing
+    =process to compute corresponding increments to liquid cloud fraction
+    =and water content, but the horizontal component of turbulent mixing
+    =(active when using the 3D Smagorinsky or blended turbulence scheme)
+    =is neglected.  Horizontal turbulent mixing is applied to vapour q,
+    =liquid water content qcl and theta, but liquid-cloud fraction is unaltered
+    =and so may become inconsistent with the other variables.
+    =
+    =If this is true, the horizontal turbulent mixing increments to
+    =q+qcl and Tl are computed and increments to qcl and cloud-fraction
+    =are calculated using the homogeneous forcing approach, consistent
+    =with the vertical mixing.
+ns=namelist/Science/UM Cloud
+sort-key=Panel-A02c
+type=logical
+
 [namelist:cloud=pc2ini]
 compulsory=true
 description=Method for doing PC2 initiation
@@ -572,6 +592,7 @@ trigger=namelist:cloud=mphys_erosion: 'pc2' ;
        =namelist:cloud=max_sigmas: 'pc2', 'bimodal' ;
        =namelist:cloud=min_sigx_ft: 'pc2', 'bimodal' ;
        =namelist:cloud=turb_var_fac_bm: 'pc2', 'bimodal' ;
+       =namelist:cloud=pc2_turb_horiz: 'pc2';
        =namelist:external_forcing=pc2_force_response: 'pc2';
 value-titles=Smith, PC2, Bimodal
 values='smith', 'pc2', 'bimodal'

interfaces/physics_schemes_interface/rose-meta/um-cloud/versions.py

@@ -31,3 +31,13 @@ def upgrade(self, config, meta_config=None):
         # Add settings
         return config, self.reports
 """
+class vn31_t464(MacroUpgrade):
+    # Upgrade macro for #464 by Ian Boutle
+
+    BEFORE_TAG = "vn3.1"
+    AFTER_TAG = "vn3.1_t464"
+
+    def upgrade(self, config, meta_config=None):
+        # Add settings
+        self.add_setting(config, ["namelist:cloud","pc2_turb_horiz"], ".false.")
+        return config, self.reports

rose-stem/app/lfric_atm/opt/rose-app-comatb.conf

@@ -1,3 +1,6 @@
+[namelist:cloud]
+pc2_turb_horiz=.true.
+
 [namelist:convection]
 !!cape_timescale=1800.0
 cv_scheme='comorph'
@@ -12,3 +15,12 @@ resdep_precipramp=.true.
 
 [namelist:microphysics]
 heavy_rain_evap_fac=5.0e7
+
+[namelist:mixing]
+method='2d_smag'
+mix_factor=0.2
+smag_l_calc='UseDx'
+smagorinsky=.true.
+
+[namelist:physics]
+smagorinsky_placement='fast'

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_comp_tran_ref_3d_l120-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         489A8379
-Inner product checksum theta =         51BA505B
-Inner product checksum u =         5F118CC2
-Inner product checksum mr1 =         4187C6A4
-Inner product checksum mr2 =         3AF39B2F
-Inner product checksum mr3 =         34D56FED
+Inner product checksum rho =         489A848F
+Inner product checksum theta =         51BA502C
+Inner product checksum u =         5F118A30
+Inner product checksum mr1 =         41879646
+Inner product checksum mr2 =         3AEB01C1
+Inner product checksum mr3 =         34AF0E7F
 Inner product checksum mr4 =                0
 Inner product checksum mr5 =                0
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_comorph_tb-C12_ex1a_cce_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         46D8528E
-Inner product checksum theta =         518BDFC0
-Inner product checksum u =         6A84BE51
-Inner product checksum mr1 =         3FCE1718
-Inner product checksum mr2 =         37F0E036
-Inner product checksum mr3 =         37746CB5
-Inner product checksum mr4 =         375B7784
-Inner product checksum mr5 =         35F5451D
+Inner product checksum rho =         46D85402
+Inner product checksum theta =         518BDFF0
+Inner product checksum u =         6A8532A1
+Inner product checksum mr1 =         3FCE1630
+Inner product checksum mr2 =         37F1B46B
+Inner product checksum mr3 =         3784215A
+Inner product checksum mr4 =         37703B68
+Inner product checksum mr5 =         361E09DD
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3-seuk_MG_ex1a_cce_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         483997DE
+Inner product checksum rho =         483997E4
 Inner product checksum theta =         51354909
-Inner product checksum u =         612F5E11
-Inner product checksum mr1 =         4090F5E4
-Inner product checksum mr2 =         358EB0BA
-Inner product checksum mr3 =         2FA519F0
-Inner product checksum mr4 =         33F4FD57
-Inner product checksum mr5 =          BE85FEC
+Inner product checksum u =         612F5E14
+Inner product checksum mr1 =         4090F5FF
+Inner product checksum mr2 =         359032DF
+Inner product checksum mr3 =         2FA50E6C
+Inner product checksum mr4 =         33F4FD6C
+Inner product checksum mr5 =          BE83A4C
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_ens-seuk_MG_ex1a_cce_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         483997E3
-Inner product checksum theta =         51354900
-Inner product checksum u =         612F5EB6
-Inner product checksum mr1 =         40910DBB
-Inner product checksum mr2 =         3600A9CC
-Inner product checksum mr3 =         2FA77A9F
-Inner product checksum mr4 =         33F4FFEC
-Inner product checksum mr5 =          BE04B21
+Inner product checksum rho =         483997E8
+Inner product checksum theta =         513548FE
+Inner product checksum u =         612F5EC2
+Inner product checksum mr1 =         40910DB9
+Inner product checksum mr2 =         35E48BF9
+Inner product checksum mr3 =         2FA740BC
+Inner product checksum mr4 =         33F4FFD2
+Inner product checksum mr5 =          BE04B62
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_mixmol-seuk_MG_ex1a_cce_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         48398C6A
+Inner product checksum rho =         48398C6F
 Inner product checksum theta =         51354302
-Inner product checksum u =         612EDCAC
-Inner product checksum mr1 =         40910338
-Inner product checksum mr2 =         35876274
-Inner product checksum mr3 =         30010BB7
-Inner product checksum mr4 =         3404AC2D
-Inner product checksum mr5 =          C14A22D
+Inner product checksum u =         612EDCC7
+Inner product checksum mr1 =         409102DC
+Inner product checksum mr2 =         358E0F08
+Inner product checksum mr3 =         30010952
+Inner product checksum mr4 =         3404AC25
+Inner product checksum mr5 =          C14A228
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         4830B653
-Inner product checksum theta =         513FD830
-Inner product checksum u =         5F74B058
-Inner product checksum mr1 =         4155BC4A
-Inner product checksum mr2 =         3A3F75A2
-Inner product checksum mr3 =         3088CA46
+Inner product checksum rho =         4830B622
+Inner product checksum theta =         513FD846
+Inner product checksum u =         5F74C1E3
+Inner product checksum mr1 =         4155AF32
+Inner product checksum mr2 =         3A507FC4
+Inner product checksum mr3 =         30A4D6EB
 Inner product checksum mr4 =                0
 Inner product checksum mr5 =                0
 Inner product checksum mr6 =                0

science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90

@@ -23,7 +23,7 @@ contains
   !> @param[in]     rho            Density
   !> @param[in]     dt             The model timestep length
   subroutine mixing_alg(mr, theta, u, derived_fields, &
-                        rho, dt)
+                        cloud_fields, rho, exner, dt)
 
     use constants_mod,                 only: i_def,r_def
     use log_mod,                       only: log_event,         &
@@ -47,7 +47,8 @@ contains
     type( field_type ), intent(inout) :: mr ( nummr )
 
     type( field_collection_type ), intent(in) :: derived_fields
-    type( field_type ), intent(in) :: rho
+    type( field_collection_type ), intent(in) :: cloud_fields
+    type( field_type ), intent(in) :: rho, exner
 
     ! Timestepping information
     real( r_def ),      intent(in) :: dt
@@ -68,7 +69,8 @@ contains
       call du%initialise( u%get_function_space() )
       call invoke( setval_c(dtheta, 0.0_r_def), &
                    setval_c(du, 0.0_r_def) )
-      call smagorinsky_alg(dtheta, du, mr, theta, u, derived_fields, rho, dt)
+      call smagorinsky_alg(dtheta, du, mr, theta, u, derived_fields, &
+                           cloud_fields, rho, exner, dt)
       call invoke( inc_X_plus_Y( theta, dtheta ), &
                    inc_X_plus_Y( u,     du ) )
 

science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 → science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.X90

@@ -18,7 +18,14 @@ module smagorinsky_alg_mod
                                        boundary_layer_um
   use mixing_config_mod,         only: mix_factor, smag_l_calc, &
                                        smag_l_calc_use_geo
-
+  use physics_config_mod,        only: smagorinsky_placement, &
+                                       smagorinsky_placement_outer, &
+                                       smagorinsky_placement_end
+  use cloud_config_mod,          only: pc2_turb_horiz
+  use physics_mappings_alg_mod,  only: map_physics_scalars
+#ifdef UM_PHYSICS
+  use pc2_force_response_alg_mod, only: pc2_force_response_alg
+#endif
   use tracer_smagorinsky_diff_kernel_mod, only: tracer_smagorinsky_diff_kernel_type
   use momentum_smagorinsky_kernel_mod,    only: momentum_smagorinsky_kernel_type
 
@@ -56,7 +63,7 @@ contains
   !> @param[in]     rho            Density
   !> @param[in]     dt             The model timestep length
 subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, &
-                           derived_fields, rho, dt )
+                           derived_fields, cloud_fields, rho, exner, dt )
 
     implicit none
 
@@ -65,8 +72,9 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, &
     type( field_type ), intent(inout) :: mr ( nummr )
 
     type( field_collection_type ), intent(in) :: derived_fields
+    type( field_collection_type ), intent(in) :: cloud_fields
 
-    type( field_type ), intent(in) :: rho, theta, u
+    type( field_type ), intent(in) :: rho, theta, u, exner
     real( r_def ),      intent(in) :: dt
 
     ! Increments on state fields
@@ -82,6 +90,13 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, &
     type( field_type ), pointer :: dx_at_w2 => null()
     type( field_type ), pointer :: panel_id => null()
     type( field_type ) :: dmr_v, dmr_cl
+#ifdef UM_PHYSICS
+    type( field_type ) :: theta_b4smag, cfl_force, bcf_force
+    type( field_type ), pointer :: liquid_fraction
+    type( field_type ), pointer :: bulk_fraction
+    type( field_type ), pointer :: theta_star
+    type( field_type ), pointer :: exner_wth
+#endif
 
     type( mesh_type ), pointer :: mesh => null()
 
@@ -217,14 +232,61 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, &
 
     end if ! With or without BL scheme
 
+    ! Scale by dt to obtain full increments
+    call invoke(inc_a_times_X(dt, dtheta), &
+                inc_a_times_X(dt, du))
     if (use_moisture) then
-      call invoke( inc_X_plus_bY( mr(imr_v), dt, dmr_v ),    &
-                   inc_X_plus_bY( mr(imr_cl), dt, dmr_cl ) )
+      call invoke(inc_a_times_X(dt, dmr_v), &
+                  inc_a_times_X(dt, dmr_cl) )
     end if
 
+#ifdef UM_PHYSICS
+    ! Call pc2 homogeneous forcing response to Smagorinsky
+    if (pc2_turb_horiz) then
+
+      ! Calculate theta before Smag - this depends on placement of Smag
+      call theta%copy_field_properties(theta_b4smag)
+      if (smagorinsky_placement == smagorinsky_placement_outer) then
+        ! In outer loop, need to get theta_star and add on the incoming inc
+        call derived_fields%get_field('theta_star', theta_star)
+        call invoke(X_plus_Y(theta_b4smag, theta_star, dtheta_io))
+      else
+        ! Just add incoming inc to incoming theta
+        call invoke(X_plus_Y(theta_b4smag, theta, dtheta_io))
+      end if
+
+      ! Set cloud fraction forcing to 0
+      call theta%copy_field_properties(cfl_force)
+      call theta%copy_field_properties(bcf_force)
+      call invoke(setval_c(cfl_force, 0.0_r_def),           &
+                  setval_c(bcf_force, 0.0_r_def))
+
+      ! If we're at the end of timestep, we need to update exner_wth
+      if (smagorinsky_placement == smagorinsky_placement_end) then
+        call derived_fields%get_field('exner_in_wth', exner_wth)
+        call map_physics_scalars(exner_wth, exner)
+      end if
+
+      ! dtheta and dmr will be updated with the pc2 response in the alg
+      call pc2_force_response_alg(mr, theta_b4smag, derived_fields,    &
+                                  cloud_fields, dtheta, dmr_v, dmr_cl, &
+                                  cfl_force, bcf_force, dt)
+
+      ! update cloud fractions with response
+      call cloud_fields%get_field('liquid_fraction', liquid_fraction)
+      call cloud_fields%get_field('bulk_fraction', bulk_fraction)
+      call invoke(inc_X_plus_Y(liquid_fraction, cfl_force), &
+                  inc_X_plus_Y(bulk_fraction, bcf_force))
+    end if
+#endif
+
     ! Apply increments from Smagorinsky mixing
-    call invoke( inc_X_plus_bY( dtheta_io, dt, dtheta ),      &
-                 inc_X_plus_bY( du_io,     dt, du ) )
+    call invoke( inc_X_plus_Y( dtheta_io, dtheta ),      &
+                 inc_X_plus_Y( du_io,     du ) )
+    if (use_moisture) then
+      call invoke( inc_X_plus_Y( mr(imr_v), dmr_v ),    &
+                   inc_X_plus_Y( mr(imr_cl), dmr_cl ) )
+    end if
 
     if (write_conservation_diag .and. use_moisture) &
       call moisture_conservation_alg( rho, mr, 'After mixing' )

science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90

@@ -392,7 +392,7 @@ contains
       call derived_fields%get_field('theta_star', theta_star)
       call derived_fields%get_field('u_star', u_star)
       call smagorinsky_alg(dtheta, du, mr, theta_star, u_star, &
-                           derived_fields, rho, dt)
+                           derived_fields, cloud_fields, rho, exner, dt)
     end if
 
     call print_field_stats_alg(dtheta, LOG_LEVEL_INFO, printmin=-20.0_r_def, &

science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90

@@ -925,10 +925,10 @@ contains
                           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)
+          call smagorinsky_alg(self%dtheta, self%du, mr, self%state(igh_t),    &
+                               self%state(igh_u),                              &
+                               derived_fields, cloud_fields, self%state(igh_d),&
+                               self%state(igh_p), cast_dt)
         end if
 
         if (use_wavedynamics) then
@@ -1083,8 +1083,8 @@ contains
     ! Apply mixing
     !--------------------------------------------------------------------
     call mixing_alg(mr, self%state(igh_t),                                     &
-                    self%state(igh_u), derived_fields,                         &
-                    self%state(igh_d), cast_dt )
+                    self%state(igh_u), derived_fields, cloud_fields,           &
+                    self%state(igh_d), self%state(igh_p), cast_dt )
 
     ! ----------------------------------------------------------------------
     ! Call cloud scheme to generate cloud and latent heating after pressure

science/linear/source/algorithm/timestepping/tl_si_timestep_alg_mod.x90

@@ -382,6 +382,7 @@ contains
     type( mesh_type ), intent( in ), pointer :: twod_mesh
 
     real( kind=r_def ) :: cast_dt
+    type( field_collection_type ) :: cloud_fields
 
     ! Quadrature object
     type( quadrature_xyoz_type ), pointer :: qr => null()
@@ -913,8 +914,8 @@ contains
     ! Apply mixing
     !--------------------------------------------------------------------
     call mixing_alg(mr, state(igh_t),                             &
-                    state(igh_u), derived_fields,                 &
-                    state(igh_d), cast_dt )
+                    state(igh_u), derived_fields, cloud_fields,   &
+                    state(igh_d), state(igh_p), cast_dt )
 
     ! Update derived variables for time level n+1
     if (use_moisture) then