MONC
subgrid_profile_diagnostics.F90
Go to the documentation of this file.
8  use state_mod, only : model_state_type
17 
18 
19  implicit none
20 
21 #ifndef TEST_MODE
22  private
23 #endif
24 
25  integer :: total_points, iqv=0, iql=0, iqr=0, iqi=0, iqs=0, &
26  iqg=0
27 
28 ! diagnostic variables
29  real(kind=default_precision), dimension(:), allocatable :: &
32  ! subgrid tke fluxes
36  ! subgrid moisture fluxes
39 ! These arrays should in due course be allocated conditionally,
40 ! but let's just get it working first.
41  real(kind=default_precision), dimension(:), allocatable :: &
44  real(kind=default_precision), dimension(:,:), allocatable :: &
45  subke_2d
46 ! Constants here provisionally
48  real(kind=default_precision) :: qlcrit
49 
51 
53 
54  logical :: l_lem_dissipation_rate = .true.
55 
57 
58 contains
59 
63  subgrid_profile_diagnostics_get_descriptor%name="subgrid_profile_diagnostics"
65 
68 
71 ! Note: Multiple copies of diagnostics are required for different time processing, so
72 ! duplicate the variables as often as necessary
73  allocate(subgrid_profile_diagnostics_get_descriptor%published_fields(35))
74 
75  subgrid_profile_diagnostics_get_descriptor%published_fields(1)="uwsg_total_local"
76  subgrid_profile_diagnostics_get_descriptor%published_fields(2)="vwsg_total_local"
77  subgrid_profile_diagnostics_get_descriptor%published_fields(3)="uusg_total_local"
78  subgrid_profile_diagnostics_get_descriptor%published_fields(4)="vvsg_total_local"
79  subgrid_profile_diagnostics_get_descriptor%published_fields(5)="wwsg_total_local"
80  subgrid_profile_diagnostics_get_descriptor%published_fields(6)="tkesg_total_local"
81  subgrid_profile_diagnostics_get_descriptor%published_fields(7)="wtsg_total_local"
82  subgrid_profile_diagnostics_get_descriptor%published_fields(8)="th2sg_total_local"
83  subgrid_profile_diagnostics_get_descriptor%published_fields(9)="wqsg_total_local"
84  !SF TKE diags
85  subgrid_profile_diagnostics_get_descriptor%published_fields(10)="buoysg_total_local"
86  subgrid_profile_diagnostics_get_descriptor%published_fields(11)="sed_total_local"
87  subgrid_profile_diagnostics_get_descriptor%published_fields(12)="ssub_total_local"
88  subgrid_profile_diagnostics_get_descriptor%published_fields(13)="dissipation_total_local"
89  !===========
90  subgrid_profile_diagnostics_get_descriptor%published_fields(14)="subke"
91  subgrid_profile_diagnostics_get_descriptor%published_fields(15)="wqv_sg_total_local"
92  subgrid_profile_diagnostics_get_descriptor%published_fields(16)="wql_sg_total_local"
93  subgrid_profile_diagnostics_get_descriptor%published_fields(17)="wqr_sg_total_local"
94  subgrid_profile_diagnostics_get_descriptor%published_fields(18)="wqi_sg_total_local"
95  subgrid_profile_diagnostics_get_descriptor%published_fields(19)="wqs_sg_total_local"
96  subgrid_profile_diagnostics_get_descriptor%published_fields(20)="wqg_sg_total_local"
97 
98 ! =====================================================
99 ! 2nd, provisionally instantaneous, stream
100 
101  subgrid_profile_diagnostics_get_descriptor%published_fields(21)="i_uwsg_total_local"
102  subgrid_profile_diagnostics_get_descriptor%published_fields(22)="i_vwsg_total_local"
103  subgrid_profile_diagnostics_get_descriptor%published_fields(23)="i_uusg_total_local"
104  subgrid_profile_diagnostics_get_descriptor%published_fields(24)="i_vvsg_total_local"
105  subgrid_profile_diagnostics_get_descriptor%published_fields(25)="i_wwsg_total_local"
106  subgrid_profile_diagnostics_get_descriptor%published_fields(26)="i_tkesg_total_local"
107  subgrid_profile_diagnostics_get_descriptor%published_fields(27)="i_wtsg_total_local"
108  subgrid_profile_diagnostics_get_descriptor%published_fields(28)="i_th2sg_total_local"
109  subgrid_profile_diagnostics_get_descriptor%published_fields(29)="i_wqsg_total_local"
110 ! =====================================================
111 
112  subgrid_profile_diagnostics_get_descriptor%published_fields(30)="wkesg_total_local"
113  subgrid_profile_diagnostics_get_descriptor%published_fields(31)="theta_dis_total_local"
114  subgrid_profile_diagnostics_get_descriptor%published_fields(32)="viscosity_coef_total_local"
115  subgrid_profile_diagnostics_get_descriptor%published_fields(33)="diffusion_coef_total_local"
116  subgrid_profile_diagnostics_get_descriptor%published_fields(34)="richardson_number_total_local"
117  subgrid_profile_diagnostics_get_descriptor%published_fields(35)="richardson_squared_total_local"
118 
119 
121 
122  subroutine initialisation_callback(current_state)
123  type(model_state_type), target, intent(inout) :: current_state
124 
125  integer :: k
126 
127  if (.not. is_component_enabled(current_state%options_database, "smagorinsky")) then
128  call log_master_log(log_error, "Subgrid model diags requested but subgrid model not enabled - check config")
129  end if
130 
131  vertical_grid=current_state%global_grid%configuration%vertical
132 
133  total_points=(current_state%local_grid%size(y_index) * current_state%local_grid%size(x_index))
134 
135  allocate(uwsg_tot(current_state%local_grid%size(z_index)) &
136  , vwsg_tot(current_state%local_grid%size(z_index)) &
137  , uusg_tot(current_state%local_grid%size(z_index)) &
138  , vvsg_tot(current_state%local_grid%size(z_index)) &
139  , wwsg_tot(current_state%local_grid%size(z_index)) &
140  , tkesg_tot(current_state%local_grid%size(z_index)) &
141  , wkesg_tot(current_state%local_grid%size(z_index)) &
142  , vis_coef_tot(current_state%local_grid%size(z_index)) &
143  , diff_coef_tot(current_state%local_grid%size(z_index)) &
144  , richardson_number_tot(current_state%local_grid%size(z_index)) &
145  , richardson_squared_tot(current_state%local_grid%size(z_index)) &
146  , dissipation_tot(current_state%local_grid%size(z_index) ) &
147  , ssub_tot(current_state%local_grid%size(z_index)) &
148  , sed_tot(current_state%local_grid%size(z_index)) &
149  , buoysg_tot(current_state%local_grid%size(z_index)) )
150 
151 ! Check consistency of any changes with the deallocations later.
152  allocate(ssq(current_state%local_grid%size(z_index)))
153  allocate(elamr_sq(current_state%local_grid%size(z_index)))
154  allocate(richardson_number(current_state%local_grid%size(z_index)))
155  allocate(subgrid_tke(current_state%local_grid%size(z_index)))
156  ! subgrid tke 2d scalar field, included in this component to prevent copying
157  ! of code. Need to revisit this
158  allocate(subke_2d(current_state%local_grid%size(y_index),current_state%local_grid%size(x_index)))
159  if (current_state%th%active) then
160  allocate(epsth(current_state%local_grid%size(z_index)) &
161  , theta_dis_tot(current_state%local_grid%size(z_index)) &
162  , wtsg_tot(current_state%local_grid%size(z_index)) &
163  , th2sg_tot(current_state%local_grid%size(z_index)))
164  endif
165 
166  if (.not. current_state%passive_q .and. current_state%number_q_fields .gt. 0) then
167  allocate(wqsg_tot(current_state%local_grid%size(z_index)))
168  iqv=get_q_index(standard_q_names%VAPOUR, 'subgrid_profile_diags')
169  iql=get_q_index(standard_q_names%CLOUD_LIQUID_MASS, 'subgrid_profile_diags')
170  qlcrit=options_get_real(current_state%options_database, "qlcrit")
171  allocate(wqv_sg_tot(current_state%local_grid%size(z_index)), &
172  wql_sg_tot(current_state%local_grid%size(z_index)))
173  endif
174  if (current_state%rain_water_mixing_ratio_index > 0) then
175  iqr = current_state%rain_water_mixing_ratio_index
176  allocate(wqr_sg_tot(current_state%local_grid%size(z_index)))
177  endif
178  if (current_state%ice_water_mixing_ratio_index > 0) then
179  iqi = current_state%ice_water_mixing_ratio_index
180  allocate(wqi_sg_tot(current_state%local_grid%size(z_index)))
181  endif
182  if (current_state%snow_water_mixing_ratio_index > 0) then
183  iqs = current_state%snow_water_mixing_ratio_index
184  allocate(wqs_sg_tot(current_state%local_grid%size(z_index)))
185  endif
186  if (current_state%graupel_water_mixing_ratio_index > 0) then
187  iqg = current_state%graupel_water_mixing_ratio_index
188  allocate(wqg_sg_tot(current_state%local_grid%size(z_index)))
189  endif
190 
191  diagnostic_generation_frequency=options_get_integer(current_state%options_database, "sampling_frequency")
192  l_lem_dissipation_rate=options_get_logical(current_state%options_database, "l_lem_dissipation_rate")
193 
194  end subroutine initialisation_callback
195 
196  subroutine timestep_callback(current_state)
197  type(model_state_type), target, intent(inout) :: current_state
198 
199  integer :: k, i, j
200  integer :: jcol, icol, target_x_index, target_y_index
201  integer :: top_index
202  real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX)) :: &
203  s11, s22, s33, s12, s23, s13, s33_on_p
204  real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX)) :: &
205  tau11,tau22, tau33, tau12, tau23, tau13, tau33_on_p
206  real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX)) :: u_i_prime_tau_i
207  real(kind=default_precision) :: vistmp, vis12, vis12a, visonp2, visonp2a, vis13, vis23
208  real(kind=default_precision), dimension(current_state%local_grid%size(Z_INDEX)) :: umean, wu_umean, vmean, wv_vmean, &
209  w_pprime_at_w, rke1, w_qvprime_at_w, w_qclprime_at_w, w_thprime_at_w, wq, rho, rec_rho, buoy_cof, &
210  uw_tot, vw_tot,qvprime_at_w,qclprime_at_w,thprime_at_w, sed_eq13, sed_eq23, sed33,uwsg,vwsg,sed_swap
211  real(kind=default_precision) :: c_virtual
212  real(kind=default_precision) :: upr_at_w,vpr_at_w, vprime_w_local, uprime_w_local, u_1, u_1_bar
213  real(kind=default_precision) :: buoy_prod, sg_shear_prod, dissipation_rate
214 
215  logical :: use_Ri_for_buoyant_prod=.true.
216 
217  c_virtual = (ratio_mol_wts-1.0_default_precision)
218  jcol=current_state%column_local_y
219  icol=current_state%column_local_x
220  target_y_index=jcol-current_state%local_grid%halo_size(y_index)
221  target_x_index=icol-current_state%local_grid%halo_size(x_index)
222 
223  if (mod(current_state%timestep, diagnostic_generation_frequency) == 0) then
224  if (current_state%first_timestep_column) then
225  sed_tot(:) = 0.0_default_precision
226  ssub_tot(:) = 0.0_default_precision
227  buoysg_tot(:) = 0.0_default_precision
228  sed_eq23(:) = 0.0_default_precision
229  sed_eq13(:) = 0.0_default_precision
230  sed_swap(:) = 0.0_default_precision
231  uwsg_tot(:) = 0.0_default_precision
232  vwsg_tot(:) = 0.0_default_precision
233  uusg_tot(:) = 0.0_default_precision
234  vvsg_tot(:) = 0.0_default_precision
235  wwsg_tot(:) = 0.0_default_precision
236  tkesg_tot(:) = 0.0_default_precision
237  vis_coef_tot(:) = 0.0_default_precision
238  diff_coef_tot(:)= 0.0_default_precision
239  dissipation_tot(:) = 0.0_default_precision
240  wkesg_tot(:) = 0.0_default_precision
241  richardson_number_tot(:) = 0.0_default_precision
242  richardson_squared_tot(:) = 0.0_default_precision
243  subke_2d(:,:) = 0.0_default_precision
244  if (current_state%th%active) then
245  wtsg_tot(:) = 0.0_default_precision
246  th2sg_tot(:) = 0.0_default_precision
247  theta_dis_tot(:)= 0.0_default_precision
248  endif
249  if (.not. current_state%passive_q .and. &
250  current_state%number_q_fields .gt. 0) then
251  wqsg_tot(:) = 0.0_default_precision
252  wqv_sg_tot(:)=0.0_default_precision
253  wql_sg_tot(:)=0.0_default_precision
254  if (iqr > 0) wqr_sg_tot(:) = 0.0_default_precision
255  if (iqi > 0) wqi_sg_tot(:) = 0.0_default_precision
256  if (iqs > 0) wqs_sg_tot(:) = 0.0_default_precision
257  if (iqg > 0) wqg_sg_tot(:) = 0.0_default_precision
258  endif
259  end if
260 
261  if (.not. current_state%halo_column) then
262 
263  do k=2,current_state%local_grid%size(z_index)-1
264 
265 ! Consistent with ssq calculation from calculate_half_squared_strain_rate
266 ! It could usefully be moved into the Smagorinsky component as a subroutine.
267 
268 #ifdef U_ACTIVE
269  ! Average over 2 points W but 0.5 cancels 2
270  s11(k)=current_state%global_grid%configuration%horizontal%cx*(&
271  (current_state%u%data(k+1, jcol, icol)-&
272  current_state%u%data(k+1, jcol, icol-1))+&
273  (current_state%u%data(k, jcol, icol)-&
274  current_state%u%data(k, jcol, icol-1)))
275 #else
276  s11(k)=0.0_default_precision
277 #endif
278 #ifdef V_ACTIVE
279  ! Average over 2 points W but 0.5 cancels 2
280  s22(k)=current_state%global_grid%configuration%horizontal%cy*(&
281  (current_state%v%data(k+1, jcol, icol)-&
282  current_state%v%data(k+1, jcol-1, icol))+&
283  (current_state%v%data(k, jcol, icol)-&
284  current_state%v%data(k, jcol-1, icol)))
285 #else
286  s22(k)=0.0_default_precision
287 #endif
288 #ifdef W_ACTIVE
289  ! Average over 2 points W but 0.5 cancels 2
290  s33(k)=((current_state%w%data(k, jcol, icol)-&
291  current_state%w%data(k-1, jcol, icol))*&
292  current_state%global_grid%configuration%vertical%rdz(k)) +&
293  ((current_state%w%data(k+1, jcol, icol)-&
294  current_state%w%data(k, jcol, icol))*&
295  current_state%global_grid%configuration%vertical%rdz(k+1))
296 
297  ! No average needed
298  s33_on_p(k) = 2.0_default_precision * &
299  ((current_state%w%data(k, jcol, icol)-&
300  current_state%w%data(k-1, jcol, icol))*&
301  current_state%global_grid%configuration%vertical%rdz(k))
302 #else
303  s33(k)=0.0_default_precision
304  s33_on_p(k) = 0.0_default_precision
305 #endif
306 #if defined(U_ACTIVE) && defined(W_ACTIVE)
307  ! Average over 2 points U and W
308  s13(k)=(((current_state%u%data(k+1, jcol, icol)-&
309  current_state%u%data(k, jcol, icol))*&
310  current_state%global_grid%configuration%vertical%rdzn(k+1)+&
311  (current_state%w%data(k, jcol, icol+1)-&
312  current_state%w%data(k, jcol, icol))*&
313  current_state%global_grid%configuration%horizontal%cx)+&
314  ((current_state%u%data(k+1, jcol, icol-1)-&
315  current_state%u%data(k, jcol, icol-1))*&
316  current_state%global_grid%configuration%vertical%rdzn(k+1)+&
317  (current_state%w%data(k, jcol, icol)-&
318  current_state%w%data(k, jcol, icol-1))*&
319  current_state%global_grid%configuration%horizontal%cx))*0.5_default_precision
320 #elif defined(U_ACTIVE) && !defined(W_ACTIVE)
321  s13(k)=(((current_state%u%data(k+1, jcol, icol)-&
322  current_state%u%data(k, jcol, icol))*&
323  current_state%global_grid%configuration%vertical%rdzn(k+1))+&
324  ((current_state%u%data(k+1, jcol, icol-1)-&
325  current_state%u%data(k, jcol, icol-1))*&
326  current_state%global_grid%configuration%vertical%rdzn(k+1)))*0.5_default_precision
327 #elif !defined(U_ACTIVE) && defined(W_ACTIVE)
328  s13(k)=(((current_state%w%data(k, jcol, icol+1)-&
329  current_state%w%data(k, jcol, icol))*&
330  current_state%global_grid%configuration%horizontal%cx)+&
331  ((current_state%w%data(k, jcol, icol)-&
332  current_state%w%data(k, jcol, icol-1))*&
333  current_state%global_grid%configuration%horizontal%cx))*0.5_default_precision
334 #else
335  s13(k)=0.0_default_precision
336 #endif
337 #if defined(W_ACTIVE) && defined(V_ACTIVE)
338  ! Average over 2 points W and V
339  s23(k)=(((current_state%w%data(k, jcol, icol) - &
340  current_state%w%data(k, jcol-1, icol)) * &
341  current_state%global_grid%configuration%horizontal%cy + &
342  (current_state%v%data(k+1, jcol-1, icol) - &
343  current_state%v%data(k, jcol-1, icol)) * &
344  current_state%global_grid%configuration%vertical%rdzn(k+1)) + &
345  ((current_state%w%data(k, jcol+1, icol) - &
346  current_state%w%data(k, jcol, icol)) * &
347  current_state%global_grid%configuration%horizontal%cy + &
348  (current_state%v%data(k+1, jcol, icol)-&
349  current_state%v%data(k, jcol, icol))*&
350  current_state%global_grid%configuration%vertical%rdzn(k+1)))*0.5_default_precision
351 #elif defined(W_ACTIVE) && !defined(V_ACTIVE)
352  s23(k)=(((current_state%w%data(k, jcol, icol)-&
353  current_state%w%data(k, jcol-1, icol))*&
354  current_state%global_grid%configuration%horizontal%cy)+&
355  ((current_state%w%data(k, jcol+1, icol)-&
356  current_state%w%data(k, jcol, icol))*&
357  current_state%global_grid%configuration%horizontal%cy))*0.5_default_precision
358 #elif !defined(W_ACTIVE) && defined(V_ACTIVE)
359  s23(k)=(((current_state%v%data(k+1, jcol-1, icol)-&
360  current_state%v%data(k, jcol-1, icol))*&
361  current_state%global_grid%configuration%vertical%rdzn(k+1))+&
362  ((current_state%v%data(k+1, jcol, icol)-&
363  current_state%v%data(k, jcol, icol))*&
364  current_state%global_grid%configuration%vertical%rdzn(k+1)))*0.5_default_precision
365 #else
366  s23(k)=0.0_default_precision
367 #endif
368 #if defined(U_ACTIVE) && defined(V_ACTIVE)
369  ! Average over 8 points from U and V
370  s12(k)=(((((current_state%u%data(k, jcol, icol-1)-&
371  current_state%u%data(k, jcol-1, icol-1))*&
372  current_state%global_grid%configuration%horizontal%cy+&
373  (current_state%v%data(k, jcol-1, icol)-&
374  current_state%v%data(k, jcol-1, icol-1))*&
375  current_state%global_grid%configuration%horizontal%cx) +&
376  ((current_state%u%data(k, jcol+1, icol-1)-&
377  current_state%u%data(k, jcol, icol-1))*&
378  current_state%global_grid%configuration%horizontal%cy+&
379  (current_state%v%data(k, jcol, icol)-&
380  current_state%v%data(k, jcol, icol-1))*&
381  current_state%global_grid%configuration%horizontal%cx)) +(&
382  ((current_state%u%data(k, jcol, icol)-&
383  current_state%u%data(k, jcol-1, icol))*&
384  current_state%global_grid%configuration%horizontal%cy+&
385  (current_state%v%data(k, jcol-1, icol+1)-&
386  current_state%v%data(k, jcol-1, icol))*&
387  current_state%global_grid%configuration%horizontal%cx) +&
388  ((current_state%u%data(k, jcol+1, icol)-&
389  current_state%u%data(k, jcol, icol))*&
390  current_state%global_grid%configuration%horizontal%cy+&
391  (current_state%v%data(k, jcol, icol+1)-&
392  current_state%v%data(k, jcol, icol))*&
393  current_state%global_grid%configuration%horizontal%cx)))+((&
394  ((current_state%u%data(k+1, jcol, icol-1)-&
395  current_state%u%data(k+1, jcol-1, icol-1))*&
396  current_state%global_grid%configuration%horizontal%cy+&
397  (current_state%v%data(k+1, jcol-1, icol)-&
398  current_state%v%data(k+1, jcol-1, icol-1))*&
399  current_state%global_grid%configuration%horizontal%cx)+&
400  ((current_state%u%data(k+1, jcol+1, icol-1)-&
401  current_state%u%data(k+1, jcol, icol-1))*&
402  current_state%global_grid%configuration%horizontal%cy+&
403  (current_state%v%data(k+1, jcol, icol)-&
404  current_state%v%data(k+1, jcol, icol-1))*&
405  current_state%global_grid%configuration%horizontal%cx))+(&
406  ((current_state%u%data(k+1, jcol, icol)-&
407  current_state%u%data(k+1, jcol-1, icol))*&
408  current_state%global_grid%configuration%horizontal%cy+&
409  (current_state%v%data(k+1, jcol-1, icol+1)-&
410  current_state%v%data(k+1, jcol-1, icol))*&
411  current_state%global_grid%configuration%horizontal%cx)+&
412  ((current_state%u%data(k+1, jcol+1, icol)-&
413  current_state%u%data(k+1, jcol, icol))*&
414  current_state%global_grid%configuration%horizontal%cy+&
415  (current_state%v%data(k+1, jcol, icol+1)-&
416  current_state%v%data(k+1, jcol, icol))*&
417  current_state%global_grid%configuration%horizontal%cx))))*0.125_default_precision
418 
419 #elif defined(U_ACTIVE) && !defined(V_ACTIVE)
420 
421  s12(k)=(((((current_state%u%data(k, jcol, icol-1)-&
422  current_state%u%data(k, jcol-1, icol-1))*&
423  current_state%global_grid%configuration%horizontal%cy) +&
424  ((current_state%u%data(k, jcol+1, icol-1)-&
425  current_state%u%data(k, jcol, icol-1))*&
426  current_state%global_grid%configuration%horizontal%cy)) +(&
427  ((current_state%u%data(k, jcol, icol)-&
428  current_state%u%data(k, jcol-1, icol))*&
429  current_state%global_grid%configuration%horizontal%cy) +&
430  ((current_state%u%data(k, jcol+1, icol)-&
431  current_state%u%data(k, jcol, icol))*&
432  current_state%global_grid%configuration%horizontal%cy)))+((&
433  ((current_state%u%data(k+1, jcol, icol-1)-&
434  current_state%u%data(k+1, jcol-1, icol-1))*&
435  current_state%global_grid%configuration%horizontal%cy)+&
436  ((current_state%u%data(k+1, jcol+1, icol-1)-&
437  current_state%u%data(k+1, jcol, icol-1))*&
438  current_state%global_grid%configuration%horizontal%cy))+(&
439  ((current_state%u%data(k+1, jcol, icol)-&
440  current_state%u%data(k+1, jcol-1, icol))*&
441  current_state%global_grid%configuration%horizontal%cy)+&
442  ((current_state%u%data(k+1, jcol+1, icol)-&
443  current_state%u%data(k+1, jcol, icol))*&
444  current_state%global_grid%configuration%horizontal%cy))))*0.125_default_precision
445 
446 #elif !defined(U_ACTIVE) && defined(V_ACTIVE)
447 
448  s12(k)=(((((current_state%v%data(k, jcol-1, icol)-&
449  current_state%v%data(k, jcol-1, icol-1))*&
450  current_state%global_grid%configuration%horizontal%cx) +&
451  ((current_state%v%data(k, jcol, icol)-&
452  current_state%v%data(k, jcol, icol-1))*&
453  current_state%global_grid%configuration%horizontal%cx)) +&
454  (((current_state%v%data(k, jcol-1, icol+1)-&
455  current_state%v%data(k, jcol-1, icol))*&
456  current_state%global_grid%configuration%horizontal%cx) +&
457  ((current_state%v%data(k, jcol, icol+1)-&
458  current_state%v%data(k, jcol, icol))*&
459  current_state%global_grid%configuration%horizontal%cx)))+((&
460  ((current_state%v%data(k+1, jcol-1, icol)-&
461  current_state%v%data(k+1, jcol-1, icol-1))*&
462  current_state%global_grid%configuration%horizontal%cx)+&
463  ((current_state%v%data(k+1, jcol, icol)-&
464  current_state%v%data(k+1, jcol, icol-1))*&
465  current_state%global_grid%configuration%horizontal%cx))+(&
466  ((current_state%v%data(k+1, jcol-1, icol+1)-&
467  current_state%v%data(k+1, jcol-1, icol))*&
468  current_state%global_grid%configuration%horizontal%cx)+&
469  ((current_state%v%data(k+1, jcol, icol+1)-&
470  current_state%v%data(k+1, jcol, icol))*&
471  current_state%global_grid%configuration%horizontal%cx))))*0.125_default_precision
472 #else
473  s12(k)=0.0_default_precision
474 #endif
475  enddo
476 
477 #if defined(U_ACTIVE)
478  ! Average over 2 points U and W
479  s13(1)=(current_state%u%data(2, jcol, icol) + &
480  current_state%u%data(2, jcol, icol-1)) * &
481  0.5_default_precision / current_state%global_grid%configuration%vertical%zn(2)
482 #else
483  s13(1)=0.0_default_precision
484 #endif
485 #if defined(V_ACTIVE)
486  ! Average over 2 points W and V
487  s23(1)=(current_state%v%data(2, jcol, icol) + &
488  current_state%v%data(2, jcol-1, icol)) * &
489  0.5_default_precision / current_state%global_grid%configuration%vertical%zn(2)
490 #else
491  s23(1)=0.0_default_precision
492 #endif
493  s12(1)=0.0_default_precision
494 
495  do k=1,current_state%local_grid%size(z_index)-1
496  tau11(k) = current_state%global_grid%configuration%vertical%rho(k) * &
497  current_state%vis_coefficient%data(k,jcol,icol) * &
498  s11(k)
499  tau22(k) = current_state%global_grid%configuration%vertical%rho(k) * &
500  current_state%vis_coefficient%data(k,jcol,icol) * &
501  s22(k)
502  tau33(k) = current_state%global_grid%configuration%vertical%rho(k) * &
503  current_state%vis_coefficient%data(k,jcol,icol) * &
504  s33(k)
505  tau12(k) = current_state%global_grid%configuration%vertical%rho(k) * &
506  current_state%vis_coefficient%data(k,jcol,icol) * &
507  s12(k)
508  tau13(k) = current_state%global_grid%configuration%vertical%rho(k) * &
509  current_state%vis_coefficient%data(k,jcol,icol) * &
510  s13(k)
511  tau23(k) = current_state%global_grid%configuration%vertical%rho(k) * &
512  current_state%vis_coefficient%data(k,jcol,icol) * &
513  s23(k)
514  enddo
515 
516  do k=2,current_state%local_grid%size(z_index)-1
517  tau33_on_p(k) = current_state%global_grid%configuration%vertical%rhon(k) * 0.5 *&
518  (current_state%vis_coefficient%data(k-1,jcol,icol) + &
519  current_state%vis_coefficient%data(k, jcol,icol)) * &
520  s33_on_p(k)
521  enddo
522 
523 ! Subgrid diagnostics
524  do k=1, current_state%local_grid%size(z_index)-1
525 !
526 ! Field on kth rho-level, so gradient evaluated from k+1st and kth rhon-levels. Horizontally aligned with
527 ! w-points.
528  uwsg_tot(k) = uwsg_tot(k) - 0.5 * &
529  current_state%vis_coefficient%data(k,jcol,icol-1) * &
530  ( current_state%u%data(k+1,jcol,icol-1) - &
531  current_state%u%data(k,jcol,icol-1) ) * &
532  vertical_grid%rdzn(k+1)
533  uwsg_tot(k) = uwsg_tot(k) - 0.5 * &
534  current_state%vis_coefficient%data(k,jcol,icol) * &
535  ( current_state%u%data(k+1,jcol,icol) - &
536  current_state%u%data(k,jcol,icol) ) * &
537  vertical_grid%rdzn(k+1)
538 !
539 ! Field on kth rho-level, so gradient evaluated from k+1st and kth rhon-levels. Horizontally aligned with
540 ! w-points.
541  vwsg_tot(k) = vwsg_tot(k) - 0.5 * &
542  current_state%vis_coefficient%data(k,jcol-1,icol) * &
543  ( current_state%v%data(k+1,jcol-1,icol) - &
544  current_state%v%data(k,jcol-1,icol) ) * &
545  vertical_grid%rdzn(k+1)
546  vwsg_tot(k) = vwsg_tot(k) - 0.5 * &
547  current_state%vis_coefficient%data(k,jcol,icol) * &
548  ( current_state%v%data(k+1,jcol,icol) - &
549  current_state%v%data(k,jcol,icol) ) * &
550  vertical_grid%rdzn(k+1)
551  !
552  ! viscosity and diffusion coefficients
553  vis_coef_tot(k) = vis_coef_tot(k) + &
554  current_state%vis_coefficient%data(k,jcol,icol)
555  diff_coef_tot(k) = diff_coef_tot(k) + &
556  current_state%diff_coefficient%data(k,jcol,icol)
557  enddo
558  ! Field on kth rho-level, so gradient evaluated from k+1st and kth rhon-levels. Horizontally aligned with
559  ! w-points.
560  if (current_state%th%active) then
561  do k=1, current_state%local_grid%size(z_index)-1
562  wtsg_tot(k) = wtsg_tot(k) - &
563  current_state%diff_coefficient%data(k,jcol,icol) * &
564  ( current_state%th%data(k+1,jcol,icol) + &
565  vertical_grid%thref(k+1) - &
566  current_state%th%data(k,jcol,icol) - &
567  vertical_grid%thref(k) ) * &
568  vertical_grid%rdzn(k+1)
569  !
570  enddo
571  endif
572 
573  if (current_state%th%active .and. .not. current_state%passive_q .and. current_state%number_q_fields .gt. 0) then
574 
575  do k=1, current_state%local_grid%size(z_index)-1
576  wqv_sg_tot(k) = wqv_sg_tot(k) + &
577  current_state%diff_coefficient%data(k, jcol, icol)* &
578  vertical_grid%rdzn(k+1) * &
579  (current_state%q(iqv)%data(k,jcol,icol) - current_state%q(iqv)%data(k+1,jcol,icol))
580  wql_sg_tot(k) = wql_sg_tot(k) + &
581  current_state%diff_coefficient%data(k, jcol, icol)* &
582  vertical_grid%rdzn(k+1) * &
583  (current_state%q(iql)%data(k,jcol,icol) - current_state%q(iql)%data(k+1,jcol,icol))
584  enddo
585  if (iqr > 0) then
586  do k=1, current_state%local_grid%size(z_index)-1
587  wqr_sg_tot(k) = wqr_sg_tot(k) + &
588  current_state%diff_coefficient%data(k, jcol, icol)* &
589  vertical_grid%rdzn(k+1) * &
590  (current_state%q(iqr)%data(k,jcol,icol) - current_state%q(iqr)%data(k+1,jcol,icol))
591  enddo
592  endif
593  if (iqi > 0) then
594  do k=1, current_state%local_grid%size(z_index)-1
595  wqi_sg_tot(k) = wqi_sg_tot(k) + &
596  current_state%diff_coefficient%data(k, jcol, icol)* &
597  vertical_grid%rdzn(k+1) * &
598  (current_state%q(iqi)%data(k,jcol,icol) - current_state%q(iqi)%data(k+1,jcol,icol))
599  enddo
600  endif
601  if (iqi > 0) then
602  do k=1, current_state%local_grid%size(z_index)-1
603  wqs_sg_tot(k) = wqs_sg_tot(k) + &
604  current_state%diff_coefficient%data(k, jcol, icol)* &
605  vertical_grid%rdzn(k+1) * &
606  (current_state%q(iqs)%data(k,jcol,icol) - current_state%q(iqs)%data(k+1,jcol,icol))
607  enddo
608  endif
609  if (iqg > 0) then
610  do k=1, current_state%local_grid%size(z_index)-1
611  wqg_sg_tot(k) = wqg_sg_tot(k) + &
612  current_state%diff_coefficient%data(k, jcol, icol)* &
613  vertical_grid%rdzn(k+1) * &
614  (current_state%q(iqg)%data(k,jcol,icol) - current_state%q(iqg)%data(k+1,jcol,icol))
615  enddo
616  endif
617  endif
618 
619 ! Level 1 wind speed - needed in various places
620 
621  u_1=sqrt(current_state%u%data(2,jcol,icol) * &
622  current_state%u%data(2,jcol,icol) + &
623  current_state%v%data(2,jcol,icol) * &
624  current_state%v%data(2,jcol,icol))
625  u_1_bar=sqrt(current_state%global_grid%configuration%vertical%olubar(2)*&
626  current_state%global_grid%configuration%vertical%olubar(2)+&
627  current_state%global_grid%configuration%vertical%olvbar(2)*&
628  current_state%global_grid%configuration%vertical%olvbar(2))
629 
630  ! Do p levels - calculate subgrid TKE diagnostics based on tau, interpolate to w levels
631  do k=1, current_state%local_grid%size(z_index)-1
632  rho(k) = current_state%global_grid%configuration%vertical%rho(k)
633  end do
634 
635 ! Buoyancy coefficient on u (zn) levels
636 ! note: this is not official MONC elsewhere
637  do k=1, current_state%local_grid%size(z_index)
638  buoy_cof(k)=g/current_state%global_grid%configuration%vertical%thref(k)
639  end do
640 
641 
642 ! =======================================================
643 ! Calculation of subgrid diagnostics dependent on the dissipation.
644 !
645 ! Rationalize this with smagorinsky.
646  ri_crit = 0.25_default_precision
647 !
648  ssq=calculate_half_squared_strain_rate(current_state, current_state%u, current_state%v, current_state%w)
649  richardson_number=calculate_richardson_number(current_state, ssq, current_state%th, current_state%q)
650 !
651  do k=2, current_state%local_grid%size(z_index)-1
652  ! Mimic LEM: I think this is the square of the mixing length (Brown et al. 94, Eq 7)
653  elamr_sq(k) = 0.0
654 
655  if ( ( &
656  current_state%vis_coefficient%data(k, jcol, icol) &
657  > 0.0) .and. (richardson_number(k) < ri_crit) ) then
658  elamr_sq(k) = &
659  current_state%vis_coefficient%data(k, jcol, icol) / &
660  sqrt( 1.0 - richardson_number(k) * &
661  current_state%diff_coefficient%data(k, jcol, icol) / &
662  current_state%vis_coefficient%data(k, jcol, icol) )
663  end if
664  if (l_lem_dissipation_rate) then
665  dissipation_rate = ssq(k) * ( &
666  current_state%vis_coefficient%data(k, jcol, icol) - &
667  richardson_number(k) * &
668  current_state%diff_coefficient%data(k, jcol, icol) )
669  else ! CIRCLE-A code
670  if (use_ri_for_buoyant_prod) then
671 
672  buoy_prod = -ssq(k) * richardson_number(k) * &
673  current_state%diff_coefficient%data(k, jcol, icol)
674 
675  else
676  if (current_state%th%active .and. &
677  .not. current_state%passive_q .and. current_state%number_q_fields .gt. 0) then
678  buoy_prod = &
679  -current_state%diff_coefficient%data(k,jcol,icol) * &
680  !d/dz - so inside needs to be on p points
681  current_state%global_grid%configuration%vertical%rdzn(k+1) * &
682  !G/th_ref on p
683  (buoy_cof(k+1) * & ! 1
684  !(th')
685  (current_state%th%data(k+1,jcol,icol) + & ! 2
686  current_state%global_grid%configuration%vertical%thref(k+1) + & ! 2
687  !th_ref * (C_virtual* qv'(k)-qcl'(k))
688  current_state%global_grid%configuration%vertical%thref(k+1) * & ! 2
689  (c_virtual * current_state%q(1)%data(k+1,jcol,icol) - & ! 3
690  current_state%q(2)%data(k+1,jcol,icol))) - & ! 1
691  !diff for d/dz
692  buoy_cof(k) * & ! 1
693  (current_state%th%data(k,jcol,icol) + & ! 2
694  current_state%global_grid%configuration%vertical%thref(k) + & ! 2
695  current_state%global_grid%configuration%vertical%thref(k) * & ! 2
696  (c_virtual * current_state%q(1)%data(k,jcol,icol) - & ! 3
697  current_state%q(2)%data(k,jcol,icol))) ) ! 0
698  else
699  call log_master_log(log_error, &
700  "Subgrid diags - buoy_prod calc needs theta and q active, STOP")
701  endif
702  endif
703  ! Circle-A dissipation rate includes buoy_prod term, which is different
704  ! to the original LEM version
705  dissipation_rate = ssq(k) * &
706  current_state%vis_coefficient%data(k, jcol, icol) + &
707  buoy_prod
708 
709  buoysg_tot(k)=buoysg_tot(k) + buoy_prod
710 
711  endif
712 
713  dissipation_tot(k) = dissipation_tot(k) + dissipation_rate
714 
715  ! PAC comment: Not clear the code below is correct (where doea a2_n fit in).
716  !
717  ! This constant needs to go somewhere accessible and where it can be set. The code in the LEM
718  ! seems to go round and round in circles here.
719  a2_n = 0.23
720  subgrid_tke(k) = ( dissipation_rate * dissipation_rate * elamr_sq(k) ) ** (1.0/3.0) / a2_n
721 
722  ! Assume isotropy.
723  tkesg_tot(k) = tkesg_tot(k) + subgrid_tke(k)
724  uusg_tot(k) = uusg_tot(k) + (2.0/3.0) * subgrid_tke(k)
725  vvsg_tot(k) = vvsg_tot(k) + (2.0/3.0) * subgrid_tke(k)
726  wwsg_tot(k) = wwsg_tot(k) + (2.0/3.0) * subgrid_tke(k)
727  wkesg_tot(k) = wkesg_tot(k) + subgrid_tke(k) * &
728  current_state%w%data(k,jcol,icol)
730  richardson_squared_tot(k) = richardson_squared_tot(k) + richardson_number(k)**2.0_default_precision
731 
732  ! Calculate the integrated subgrid TKE for the scalar.
733  ! BAD place to put this but convenient
734  subke_2d(target_y_index, target_x_index) = subke_2d(target_y_index, target_x_index) + &
735  (subgrid_tke(k)*0.5_default_precision* &
736  vertical_grid%dzn(k+1)* &
737  vertical_grid%rho(k))
738  enddo
739  ! *********************** Subgrid buoyant production***************************
740  ! Note - calculating on z levels (i.e. w)
741  ! So need dzn, rdzn.
742  ! zn(k)-zn(k-1) is dzn(k), so nothing is stored in k=1
743  !------
744  if (.not. l_lem_dissipation_rate) then
745  if (current_state%th%active .and. &
746  .not. current_state%passive_q .and. current_state%number_q_fields .gt. 0) then
747  k = 1
748  buoy_prod = &
749  -current_state%diff_coefficient%data(k,jcol,icol) * &
750  !d/dz - so inside needs to be on p points
751  current_state%global_grid%configuration%vertical%rdzn(k+1) * &
752  !G/th_ref on p
753  (buoy_cof(k+1) * & ! 1
754  !(th')
755  (current_state%th%data(k+1,jcol,icol) + & ! 2
756  current_state%global_grid%configuration%vertical%thref(k+1) + & ! 2
757  !th_ref * (C_virtual* qv'(k)-qcl'(k))
758  current_state%global_grid%configuration%vertical%thref(k+1) * & ! 2
759  (c_virtual * current_state%q(1)%data(k+1,jcol,icol) & ! 3
760  - current_state%q(2)%data(k+1,jcol,icol))) - & ! 1
761  !diff for d/dz
762  buoy_cof(k) * & ! 1
763  (current_state%th%data(k,jcol,icol) + & ! 2
764  current_state%global_grid%configuration%vertical%thref(k) + & ! 2
765  current_state%global_grid%configuration%vertical%thref(k) * & ! 2
766  (c_virtual * current_state%q(1)%data(k,jcol,icol) & ! 3
767  - current_state%q(2)%data(k,jcol,icol))) ) ! 0
768 
769  buoysg_tot(k)=buoysg_tot(k) + buoy_prod
770 
771 
772  dissipation_tot(1) = dissipation_tot(1) + &
773  current_state%vis_coefficient%data(1, jcol, icol) &
774  * u_1 * u_1 / &
775  (current_state%global_grid%configuration%vertical%zn(2) * &
776  current_state%global_grid%configuration%vertical%zn(2)) + buoy_prod
777  else
778  call log_master_log(log_error, &
779  "Subgrid diags - buoy_prod calc needs theta and q active, STOP")
780  endif
781  endif
782 
783  ! Needs buoyancy part added later
784 
785  ! divide subke by altitude to make it column mean (as in the LEM)
786  subke_2d(target_y_index, target_x_index) = subke_2d(target_y_index, target_x_index)/ &
787  vertical_grid%z(current_state%local_grid%size(z_index))
788  !
789  ! Thermal equivalents:
790  ! Again, this needs to go somewhere better. pr_n is used in Smagorinsky anyway (as a local variable).
791 
792  if (current_state%th%active) then
793  epsth=calculate_thermal_dissipation_rate(current_state, current_state%th)
794  ath2_n = 0.3
795  pr_n = 0.7
796  do k=2, current_state%local_grid%size(z_index)-1
797  ! Theta dissipation rate
798  theta_dis_tot(k) = theta_dis_tot(k) + epsth(k)
799  if (subgrid_tke(k) > 0.0) &
800  th2sg_tot(k) = th2sg_tot(k) + sqrt( a2_n * elamr_sq(k) / subgrid_tke(k) ) * &
801  epsth(k) / ( ath2_n**2 * pr_n)
802  enddo
803  endif
804 
805 
806  ! *********************** Subgrid shear production***************************
807  ! Note - calculating on z levels (i.e. w)
808  !
809  do k=2,current_state%local_grid%size(z_index)-1
810 
811 
812  !Subgrid shear-------
813  umean(k)=(current_state%global_grid%configuration%vertical%olubar(k+1) - &
814  current_state%global_grid%configuration%vertical%olubar(k)) * &
815  current_state%global_grid%configuration%vertical%rdzn(k+1)
816  vmean(k)=(current_state%global_grid%configuration%vertical%olvbar(k+1) - &
817  current_state%global_grid%configuration%vertical%olvbar(k)) * &
818  current_state%global_grid%configuration%vertical%rdzn(k+1)
819 
820  sg_shear_prod = current_state%vis_coefficient%data(k,jcol,icol)* &
821  (s13(k)*umean(k)+s23(k)*vmean(k))
822 
823  ssub_tot(k)=ssub_tot(k) + sg_shear_prod
824 
825  !tau is at p points - so convert to w points
826  uwsg(k)=0.5*(tau13(k)+tau13(k+1))/rho(k)
827  vwsg(k)=0.5*(tau23(k)+tau23(k+1))/rho(k)
828 
829  end do
830 
831  k=1
832 
833  sg_shear_prod = current_state%vis_coefficient%data(1, jcol, icol) * &
834  ( 0.5 * &
835  (current_state%u%data(2,jcol,icol-1) + &
836  current_state%u%data(2,jcol,icol) ) * &
837  current_state%global_grid%configuration%vertical%olubar(2) + &
838  0.5 * &
839  (current_state%v%data(2,jcol-1,icol) + &
840  current_state%v%data(2,jcol,icol) ) * &
841  current_state%global_grid%configuration%vertical%olvbar(2) ) / &
842  (current_state%global_grid%configuration%vertical%zn(2) * &
843  current_state%global_grid%configuration%vertical%zn(2))
844 
845  ssub_tot(1)=ssub_tot(1) + sg_shear_prod
846 
847  ! *********************** Subgrid turbulence transport***************************
848  ! Note - calculating on z levels (i.e. w)
849  ! so need u_i_prime_tau_i on p levels
850 
851  do k=2, current_state%local_grid%size(z_index)
852 
853  u_i_prime_tau_i(k) = ( 0.5_default_precision * &
854  (current_state%u%data(k,jcol,icol-1) + &
855  current_state%u%data(k,jcol,icol) ) - &
856  current_state%global_grid%configuration%vertical%olubar(k)) * &
857  0.5_default_precision * ( tau13(k-1)+tau13(k))
858 
859  u_i_prime_tau_i(k) = u_i_prime_tau_i(k) + ( 0.5_default_precision * &
860  (current_state%v%data(k,jcol-1,icol) + &
861  current_state%v%data(k,jcol,icol) ) - &
862  current_state%global_grid%configuration%vertical%olvbar(k)) * &
863  0.5_default_precision * ( tau23(k-1)+tau23(k))
864 
865  u_i_prime_tau_i(k) = u_i_prime_tau_i(k) + ( 0.5_default_precision * &
866  (current_state%w%data(k-1,jcol,icol)+&
867  current_state%w%data(k,jcol,icol))) *&
868  tau33_on_p(k)
869 
870  end do
871 
872  u_i_prime_tau_i(1) = &
873  ( 0.5_default_precision * &
874  (current_state%u%data(2,jcol,icol-1) + &
875  current_state%u%data(2,jcol,icol) ) - &
876  current_state%global_grid%configuration%vertical%olubar(2) ) * &
877  ( 0.5_default_precision * &
878  (current_state%u%data(2,jcol,icol-1) + &
879  current_state%u%data(2,jcol,icol) ) )
880 
881  u_i_prime_tau_i(1) = u_i_prime_tau_i(1) + &
882  ( 0.5_default_precision * &
883  (current_state%v%data(2,jcol-1,icol) + &
884  current_state%v%data(2,jcol,icol) ) - &
885  current_state%global_grid%configuration%vertical%olvbar(2) ) * &
886  ( 0.5_default_precision * &
887  (current_state%v%data(2,jcol-1,icol) + &
888  current_state%v%data(2,jcol,icol) ) )
889 
890  ! Note - here we use the surface vis_coefficient to ensure exact cancellation in
891  ! numerical budget
892  u_i_prime_tau_i(1) = u_i_prime_tau_i(1) * &
893  current_state%vis_coefficient%data(1,jcol,icol) / &
894  current_state%global_grid%configuration%vertical%zn(2)
895 
896 
897  do k=2, current_state%local_grid%size(z_index)-1
898  sed_tot(k)=sed_tot(k) + (u_i_prime_tau_i(k+1)-u_i_prime_tau_i(k)) * &
899  current_state%global_grid%configuration%vertical%rdzn(k+1) / &
900  current_state%global_grid%configuration%vertical%rho(k)
901  end do
902 
903  sed_tot(1)=sed_tot(1) + u_i_prime_tau_i(1) / &
904  current_state%global_grid%configuration%vertical%zn(2)
905 
906  ! =======================================================
907  endif ! (.not. current_state%halo_column)
908  endif ! ( if diagnostic frequency step )
909  end subroutine timestep_callback
910 
915  subroutine field_information_retrieval_callback(current_state, name, field_information)
916  type(model_state_type), target, intent(inout) :: current_state
917  character(len=*), intent(in) :: name
918  type(component_field_information_type), intent(out) :: field_information
919 
920  field_information%field_type=component_array_field_type
921  field_information%data_type=component_double_data_type
922  if (name .eq. 'subke') then
923  field_information%number_dimensions=2
924  field_information%dimension_sizes(1)=current_state%local_grid%size(y_index)
925  field_information%dimension_sizes(2)=current_state%local_grid%size(x_index)
926  else
927  field_information%number_dimensions=1
928  field_information%dimension_sizes(1)=current_state%local_grid%size(z_index)
929  endif
930 
931  if (name .eq. "th2sg_total_local" .or. name .eq. "wtsg_total_local" &
932  .or. name .eq. "theta_dis_total_local" ) then
933  field_information%enabled=current_state%th%active
934  else if (name .eq. "wqv_sg_total_local" .or. name .eq. "wql_sg_total_local") then
935  field_information%enabled= (.not.current_state%passive_q) .and. &
936  (current_state%number_q_fields .gt. 0)
937  else if (name .eq. "wqr_sg_total_local") then
938  field_information%enabled= current_state%rain_water_mixing_ratio_index .gt. 0
939  else if (name .eq. "wqi_sg_total_local") then
940  field_information%enabled= current_state%ice_water_mixing_ratio_index .gt. 0
941  else if (name .eq. "wqs_sg_total_local") then
942  field_information%enabled= current_state%snow_water_mixing_ratio_index .gt. 0
943  else if (name .eq. "wqg_sg_total_local") then
944  field_information%enabled= current_state%graupel_water_mixing_ratio_index .gt. 0
945 
946 ! ========================================================================
947 ! 2nd stream
948  else if (name .eq. "i_th2sg_total_local") then
949  field_information%enabled=current_state%th%active
950  else if (name .eq. "i_wqsg_total_local") then
951  field_information%enabled= (.not.current_state%passive_q) .and. &
952  (current_state%liquid_water_mixing_ratio_index > 0)
953 ! ========================================================================
954  else
955  field_information%enabled=.true.
956  end if
958 
963  subroutine field_value_retrieval_callback(current_state, name, field_value)
964  type(model_state_type), target, intent(inout) :: current_state
965  character(len=*), intent(in) :: name
966  type(component_field_value_type), intent(out) :: field_value
967 
968  integer :: k
969 
970  if (name .eq. "uwsg_total_local") then
971  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
972  do k = 1, current_state%local_grid%size(z_index)
973  field_value%real_1d_array(k)=uwsg_tot(k)
974  enddo
975  else if (name .eq. "sed_total_local") then
976  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
977  do k = 1, current_state%local_grid%size(z_index)
978  field_value%real_1d_array(k)=sed_tot(k)
979  enddo
980  else if (name .eq. "ssub_total_local") then
981  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
982  do k = 1, current_state%local_grid%size(z_index)
983  field_value%real_1d_array(k)=ssub_tot(k)
984  enddo
985  else if (name .eq. "buoysg_total_local") then
986  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
987  do k = 1, current_state%local_grid%size(z_index)
988  field_value%real_1d_array(k)=buoysg_tot(k)
989  enddo
990  else if (name .eq. "vwsg_total_local") then
991  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
992  do k = 1, current_state%local_grid%size(z_index)
993  field_value%real_1d_array(k)=vwsg_tot(k)
994  enddo
995  else if (name .eq. "uusg_total_local") then
996  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
997  do k = 1, current_state%local_grid%size(z_index)
998  field_value%real_1d_array(k)=uusg_tot(k)
999  enddo
1000  else if (name .eq. "vvsg_total_local") then
1001  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1002  do k = 1, current_state%local_grid%size(z_index)
1003  field_value%real_1d_array(k)=vvsg_tot(k)
1004  enddo
1005  else if (name .eq. "wwsg_total_local") then
1006  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1007  do k = 1, current_state%local_grid%size(z_index)
1008  field_value%real_1d_array(k)=wwsg_tot(k)
1009  enddo
1010  else if (name .eq. "tkesg_total_local") then
1011  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1012  do k = 1, current_state%local_grid%size(z_index)
1013  field_value%real_1d_array(k)=tkesg_tot(k)
1014  enddo
1015  else if (name .eq. "wkesg_total_local") then
1016  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1017  do k = 1, current_state%local_grid%size(z_index)
1018  field_value%real_1d_array(k)=wkesg_tot(k)
1019  enddo
1020  else if (name .eq. "viscosity_coef_total_local") then
1021  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1022  do k = 1, current_state%local_grid%size(z_index)
1023  field_value%real_1d_array(k)=vis_coef_tot(k)
1024  enddo
1025  else if (name .eq. "diffusion_coef_total_local") then
1026  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1027  do k = 1, current_state%local_grid%size(z_index)
1028  field_value%real_1d_array(k)=diff_coef_tot(k)
1029  enddo
1030  else if (name .eq. "richardson_number_total_local") then
1031  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1032  do k = 1, current_state%local_grid%size(z_index)
1033  field_value%real_1d_array(k)=richardson_number_tot(k)
1034  enddo
1035  else if (name .eq. "richardson_squared_total_local") then
1036  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1037  do k = 1, current_state%local_grid%size(z_index)
1038  field_value%real_1d_array(k)=richardson_squared_tot(k)
1039  enddo
1040  else if (name .eq. "dissipation_total_local") then
1041  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1042  do k = 1, current_state%local_grid%size(z_index)
1043  field_value%real_1d_array(k)=dissipation_tot(k)
1044  enddo
1045  else if (name .eq. "wtsg_total_local") then
1046  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1047  do k = 1, current_state%local_grid%size(z_index)
1048  field_value%real_1d_array(k)=wtsg_tot(k)
1049  enddo
1050  else if (name .eq. "th2sg_total_local") then
1051  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1052  do k = 1, current_state%local_grid%size(z_index)
1053  field_value%real_1d_array(k)=th2sg_tot(k)
1054  enddo
1055  else if (name .eq. "theta_dis_total_local") then
1056  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1057  do k = 1, current_state%local_grid%size(z_index)
1058  field_value%real_1d_array(k)=theta_dis_tot(k)
1059  enddo
1060  else if (name .eq. "wqv_sg_total_local") then
1061  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1062  do k = 1, current_state%local_grid%size(z_index)
1063  field_value%real_1d_array(k)=wqv_sg_tot(k)
1064  enddo
1065  else if (name .eq. "wql_sg_total_local") then
1066  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1067  do k = 1, current_state%local_grid%size(z_index)
1068  field_value%real_1d_array(k)=wql_sg_tot(k)
1069  enddo
1070  else if (name .eq. "wqr_sg_total_local") then
1071  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1072  do k = 1, current_state%local_grid%size(z_index)
1073  field_value%real_1d_array(k)=wqr_sg_tot(k)
1074  enddo
1075  else if (name .eq. "wqi_sg_total_local") then
1076  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1077  do k = 1, current_state%local_grid%size(z_index)
1078  field_value%real_1d_array(k)=wqi_sg_tot(k)
1079  enddo
1080  else if (name .eq. "wqs_sg_total_local") then
1081  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1082  do k = 1, current_state%local_grid%size(z_index)
1083  field_value%real_1d_array(k)=wqs_sg_tot(k)
1084  enddo
1085  else if (name .eq. "wqg_sg_total_local") then
1086  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1087  do k = 1, current_state%local_grid%size(z_index)
1088  field_value%real_1d_array(k)=wqg_sg_tot(k)
1089  enddo
1090 ! =====================================================
1091 ! 2nd stream
1092  else if (name .eq. "i_uwsg_total_local") then
1093  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1094  do k = 1, current_state%local_grid%size(z_index)
1095  field_value%real_1d_array(k)=uwsg_tot(k)
1096  enddo
1097  else if (name .eq. "i_vwsg_total_local") then
1098  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1099  do k = 1, current_state%local_grid%size(z_index)
1100  field_value%real_1d_array(k)=vwsg_tot(k)
1101  enddo
1102  else if (name .eq. "i_uusg_total_local") then
1103  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1104  do k = 1, current_state%local_grid%size(z_index)
1105  field_value%real_1d_array(k)=uusg_tot(k)
1106  enddo
1107  else if (name .eq. "i_vvsg_total_local") then
1108  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1109  do k = 1, current_state%local_grid%size(z_index)
1110  field_value%real_1d_array(k)=vvsg_tot(k)
1111  enddo
1112  else if (name .eq. "i_wwsg_total_local") then
1113  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1114  do k = 1, current_state%local_grid%size(z_index)
1115  field_value%real_1d_array(k)=wwsg_tot(k)
1116  enddo
1117  else if (name .eq. "i_tkesg_total_local") then
1118  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1119  do k = 1, current_state%local_grid%size(z_index)
1120  field_value%real_1d_array(k)=tkesg_tot(k)
1121  enddo
1122  else if (name .eq. "i_wtsg_total_local") then
1123  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1124  do k = 1, current_state%local_grid%size(z_index)
1125  field_value%real_1d_array(k)=wtsg_tot(k)
1126  enddo
1127  else if (name .eq. "i_th2sg_total_local") then
1128  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1129  do k = 1, current_state%local_grid%size(z_index)
1130  field_value%real_1d_array(k)=th2sg_tot(k)
1131  enddo
1132  else if (name .eq. "i_wqsg_total_local") then
1133  allocate(field_value%real_1d_array(current_state%local_grid%size(z_index)))
1134  do k = 1, current_state%local_grid%size(z_index)
1135  field_value%real_1d_array(k)=wqsg_tot(k)
1136  enddo
1137  ! =====================================================
1138  else if (name .eq. 'subke') then
1139  allocate(field_value%real_2d_array(current_state%local_grid%size(y_index), &
1140  current_state%local_grid%size(x_index)))
1141  field_value%real_2d_array(:,:)=subke_2d(:,:)
1142  end if
1143  end subroutine field_value_retrieval_callback
1144 
registry_mod::is_component_enabled
logical function, public is_component_enabled(options_database, component_name)
Determines whether or not a specific component is registered and enabled.
Definition: registry.F90:334
logging_mod::log_error
integer, parameter, public log_error
Only log ERROR messages.
Definition: logging.F90:11
subgrid_profile_diagnostics_mod::wtsg_tot
real(kind=default_precision), dimension(:), allocatable wtsg_tot
Definition: subgrid_profile_diagnostics.F90:29
smagorinsky_mod::calculate_richardson_number
real(kind=default_precision) function, dimension(current_state%local_grid%size(z_index)), public calculate_richardson_number(current_state, ssq, th, q)
Calculates the richardson number depending upon the setup of the model and the method selected.
Definition: smagorinsky.F90:205
smagorinsky_mod::calculate_half_squared_strain_rate
real(kind=default_precision) function, dimension(current_state%local_grid%size(z_index)), public calculate_half_squared_strain_rate(current_state, u, v, w)
Calculates the half squared strain rate on l_w-points which is used in determining the Richardson num...
Definition: smagorinsky.F90:411
subgrid_profile_diagnostics_mod::wql_sg_tot
real(kind=default_precision), dimension(:), allocatable wql_sg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::uusg_tot
real(kind=default_precision), dimension(:), allocatable uusg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::diff_coef_tot
real(kind=default_precision), dimension(:), allocatable diff_coef_tot
Definition: subgrid_profile_diagnostics.F90:29
logging_mod::log_warn
integer, parameter, public log_warn
Log WARNING and ERROR messages.
Definition: logging.F90:12
subgrid_profile_diagnostics_mod::ssub_tot
real(kind=default_precision), dimension(:), allocatable ssub_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::subgrid_profile_diagnostics_get_descriptor
type(component_descriptor_type) function, public subgrid_profile_diagnostics_get_descriptor()
Provides the component descriptor for the core to register.
Definition: subgrid_profile_diagnostics.F90:63
subgrid_profile_diagnostics_mod::wqs_sg_tot
real(kind=default_precision), dimension(:), allocatable wqs_sg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::elamr_sq
real(kind=default_precision), dimension(:), allocatable elamr_sq
Definition: subgrid_profile_diagnostics.F90:41
subgrid_profile_diagnostics_mod::richardson_number
real(kind=default_precision), dimension(:), allocatable richardson_number
Definition: subgrid_profile_diagnostics.F90:41
grids_mod::x_index
integer, parameter, public x_index
Definition: grids.F90:14
optionsdatabase_mod::options_get_integer
integer function, public options_get_integer(options_database, key, index)
Retrieves an integer value from the database that matches the provided key.
Definition: optionsdatabase.F90:217
grids_mod::y_index
integer, parameter, public y_index
Definition: grids.F90:14
subgrid_profile_diagnostics_mod::vertical_grid
type(vertical_grid_configuration_type) vertical_grid
Definition: subgrid_profile_diagnostics.F90:52
subgrid_profile_diagnostics_mod::wkesg_tot
real(kind=default_precision), dimension(:), allocatable wkesg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::l_lem_dissipation_rate
logical l_lem_dissipation_rate
Definition: subgrid_profile_diagnostics.F90:54
subgrid_profile_diagnostics_mod::theta_dis_tot
real(kind=default_precision), dimension(:), allocatable theta_dis_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::richardson_squared_tot
real(kind=default_precision), dimension(:), allocatable richardson_squared_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::th2sg_tot
real(kind=default_precision), dimension(:), allocatable th2sg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::field_information_retrieval_callback
subroutine field_information_retrieval_callback(current_state, name, field_information)
Field information retrieval callback, this returns information for a specific components published fi...
Definition: subgrid_profile_diagnostics.F90:916
subgrid_profile_diagnostics_mod::field_value_retrieval_callback
subroutine field_value_retrieval_callback(current_state, name, field_value)
Field value retrieval callback, this returns the value of a specific published field.
Definition: subgrid_profile_diagnostics.F90:964
subgrid_profile_diagnostics_mod::vis_coef_tot
real(kind=default_precision), dimension(:), allocatable vis_coef_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::ri_crit
real(kind=default_precision) ri_crit
Definition: subgrid_profile_diagnostics.F90:47
subgrid_profile_diagnostics_mod::timestep_callback
subroutine timestep_callback(current_state)
Definition: subgrid_profile_diagnostics.F90:197
monc_component_mod
Interfaces and types that MONC components must specify.
Definition: monc_component.F90:6
subgrid_profile_diagnostics_mod::subgrid_tke
real(kind=default_precision), dimension(:), allocatable subgrid_tke
Definition: subgrid_profile_diagnostics.F90:41
subgrid_profile_diagnostics_mod::total_points
integer total_points
Definition: subgrid_profile_diagnostics.F90:25
subgrid_profile_diagnostics_mod::wqg_sg_tot
real(kind=default_precision), dimension(:), allocatable wqg_sg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::wwsg_tot
real(kind=default_precision), dimension(:), allocatable wwsg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::vwsg_tot
real(kind=default_precision), dimension(:), allocatable vwsg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::iqi
integer iqi
Definition: subgrid_profile_diagnostics.F90:25
subgrid_profile_diagnostics_mod::a2_n
real(kind=default_precision) a2_n
Definition: subgrid_profile_diagnostics.F90:47
science_constants_mod
Scientific constant values used throughout simulations. Each has a default value and this can be over...
Definition: scienceconstants.F90:3
subgrid_profile_diagnostics_mod::qlcrit
real(kind=default_precision) qlcrit
Definition: subgrid_profile_diagnostics.F90:48
grids_mod::z_index
integer, parameter, public z_index
Grid index parameters.
Definition: grids.F90:14
subgrid_profile_diagnostics_mod::sed_tot
real(kind=default_precision), dimension(:), allocatable sed_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::pr_n
real(kind=default_precision) pr_n
Definition: subgrid_profile_diagnostics.F90:47
subgrid_profile_diagnostics_mod::diagnostic_generation_frequency
integer diagnostic_generation_frequency
Definition: subgrid_profile_diagnostics.F90:50
science_constants_mod::ratio_mol_wts
real(kind=default_precision), public ratio_mol_wts
Definition: scienceconstants.F90:13
q_indices_mod::standard_q_names
type(standard_q_names_type), public standard_q_names
Definition: q_indices.F90:59
monc_component_mod::component_field_value_type
Wrapper type for the value returned for a published field from a component.
Definition: monc_component.F90:22
state_mod::model_state_type
The ModelState which represents the current state of a run.
Definition: state.F90:39
subgrid_profile_diagnostics_mod::iqs
integer iqs
Definition: subgrid_profile_diagnostics.F90:25
subgrid_profile_diagnostics_mod::wqr_sg_tot
real(kind=default_precision), dimension(:), allocatable wqr_sg_tot
Definition: subgrid_profile_diagnostics.F90:29
monc_component_mod::component_field_information_type
Definition: monc_component.F90:31
subgrid_profile_diagnostics_mod::wqi_sg_tot
real(kind=default_precision), dimension(:), allocatable wqi_sg_tot
Definition: subgrid_profile_diagnostics.F90:29
optionsdatabase_mod::options_get_logical
logical function, public options_get_logical(options_database, key, index)
Retrieves a logical value from the database that matches the provided key.
Definition: optionsdatabase.F90:154
monc_component_mod::component_scalar_field_type
integer, parameter, public component_scalar_field_type
Definition: monc_component.F90:15
subgrid_profile_diagnostics_mod::iql
integer iql
Definition: subgrid_profile_diagnostics.F90:25
subgrid_profile_diagnostics_mod::uwsg_tot
real(kind=default_precision), dimension(:), allocatable uwsg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::epsth
real(kind=default_precision), dimension(:), allocatable epsth
Definition: subgrid_profile_diagnostics.F90:41
subgrid_profile_diagnostics_mod::richardson_number_tot
real(kind=default_precision), dimension(:), allocatable richardson_number_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::subke_2d
real(kind=default_precision), dimension(:,:), allocatable subke_2d
Definition: subgrid_profile_diagnostics.F90:44
q_indices_mod::get_q_index
integer function, public get_q_index(name, assigning_component)
Add in a new entry into the register if the name does not already exist or return the index of the pr...
Definition: q_indices.F90:112
smagorinsky_mod
Calculates the Smagorinsky eddy viscosity and diffusivity at l_w-points.
Definition: smagorinsky.F90:2
logging_mod
Logging utility.
Definition: logging.F90:2
subgrid_profile_diagnostics_mod::initialisation_callback
subroutine initialisation_callback(current_state)
Definition: subgrid_profile_diagnostics.F90:123
grids_mod::vertical_grid_configuration_type
The configuration of the grid vertically.
Definition: grids.F90:28
datadefn_mod
Contains common definitions for the data and datatypes used by MONC.
Definition: datadefn.F90:2
logging_mod::log_master_log
subroutine, public log_master_log(level, message)
Will log just from the master process.
Definition: logging.F90:47
subgrid_profile_diagnostics_mod::ssq
real(kind=default_precision), dimension(:), allocatable ssq
Definition: subgrid_profile_diagnostics.F90:41
subgrid_profile_diagnostics_mod::ath2_n
real(kind=default_precision) ath2_n
Definition: subgrid_profile_diagnostics.F90:47
subgrid_profile_diagnostics_mod::iqv
integer iqv
Definition: subgrid_profile_diagnostics.F90:25
q_indices_mod
This manages the Q variables and specifically the mapping between names and the index that they are s...
Definition: q_indices.F90:2
subgrid_profile_diagnostics_mod::tkesg_tot
real(kind=default_precision), dimension(:), allocatable tkesg_tot
Definition: subgrid_profile_diagnostics.F90:29
registry_mod
MONC component registry.
Definition: registry.F90:5
smagorinsky_mod::calculate_thermal_dissipation_rate
real(kind=default_precision) function, dimension(current_state%local_grid%size(z_index)), public calculate_thermal_dissipation_rate(current_state, th)
Definition: smagorinsky.F90:625
subgrid_profile_diagnostics_mod::iqg
integer iqg
Definition: subgrid_profile_diagnostics.F90:25
subgrid_profile_diagnostics_mod::wqv_sg_tot
real(kind=default_precision), dimension(:), allocatable wqv_sg_tot
Definition: subgrid_profile_diagnostics.F90:29
monc_component_mod::component_integer_data_type
integer, parameter, public component_integer_data_type
Definition: monc_component.F90:16
subgrid_profile_diagnostics_mod::vvsg_tot
real(kind=default_precision), dimension(:), allocatable vvsg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::dissipation_tot
real(kind=default_precision), dimension(:), allocatable dissipation_tot
Definition: subgrid_profile_diagnostics.F90:29
grids_mod
Functionality to support the different types of grid and abstraction between global grids and local o...
Definition: grids.F90:5
subgrid_profile_diagnostics_mod::iqr
integer iqr
Definition: subgrid_profile_diagnostics.F90:25
subgrid_profile_diagnostics_mod
Definition: subgrid_profile_diagnostics.F90:1
optionsdatabase_mod
Manages the options database. Contains administration functions and deduce runtime options from the c...
Definition: optionsdatabase.F90:7
monc_component_mod::component_descriptor_type
Description of a component.
Definition: monc_component.F90:42
monc_component_mod::component_double_data_type
integer, parameter, public component_double_data_type
Definition: monc_component.F90:16
optionsdatabase_mod::options_get_real
real(kind=default_precision) function, public options_get_real(options_database, key, index)
Retrieves a real value from the database that matches the provided key.
Definition: optionsdatabase.F90:91
datadefn_mod::default_precision
integer, parameter, public default_precision
MPI communication type which we use for the prognostic and calculation data.
Definition: datadefn.F90:17
state_mod
The model state which represents the current state of a run.
Definition: state.F90:2
monc_component_mod::component_array_field_type
integer, parameter, public component_array_field_type
Definition: monc_component.F90:15
subgrid_profile_diagnostics_mod::buoysg_tot
real(kind=default_precision), dimension(:), allocatable buoysg_tot
Definition: subgrid_profile_diagnostics.F90:29
subgrid_profile_diagnostics_mod::wqsg_tot
real(kind=default_precision), dimension(:), allocatable wqsg_tot
Definition: subgrid_profile_diagnostics.F90:29
science_constants_mod::g
real(kind=default_precision), public g
Definition: scienceconstants.F90:13