Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Handles the type to compute averages during an MD
10 : !> \author Teodoro Laino [tlaino] - 03.2008 - University of Zurich
11 : ! **************************************************************************************************
12 : MODULE averages_types
13 : USE cell_types, ONLY: cell_type
14 : USE colvar_utils, ONLY: get_clv_force,&
15 : number_of_colvar
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_type,&
18 : cp_to_string
19 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
20 : cp_print_key_unit_nr
21 : USE force_env_types, ONLY: force_env_type
22 : USE input_section_types, ONLY: section_vals_get,&
23 : section_vals_get_subs_vals,&
24 : section_vals_remove_values,&
25 : section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: default_string_length,&
28 : dp
29 : USE md_ener_types, ONLY: md_ener_type
30 : USE virial_types, ONLY: virial_type
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : ! **************************************************************************************************
38 : TYPE average_quantities_type
39 : INTEGER :: ref_count = 0, itimes_start = -1
40 : LOGICAL :: do_averages = .FALSE.
41 : TYPE(section_vals_type), POINTER :: averages_section => NULL()
42 : ! Real Scalar Quantities
43 : REAL(KIND=dp) :: avetemp = 0.0_dp, avepot = 0.0_dp, avekin = 0.0_dp, &
44 : avevol = 0.0_dp, aveca = 0.0_dp, avecb = 0.0_dp, avecc = 0.0_dp
45 : REAL(KIND=dp) :: avetemp_baro = 0.0_dp, avehugoniot = 0.0_dp, avecpu = 0.0_dp
46 : REAL(KIND=dp) :: aveal = 0.0_dp, avebe = 0.0_dp, avega = 0.0_dp, avepress = 0.0_dp, &
47 : avekinc = 0.0_dp, avetempc = 0.0_dp, avepxx = 0.0_dp
48 : REAL(KIND=dp) :: avetemp_qm = 0.0_dp, avekin_qm = 0.0_dp, econs = 0.0_dp
49 : ! Virial
50 : TYPE(virial_type), POINTER :: virial => NULL()
51 : ! Colvar
52 : REAL(KIND=dp), POINTER, DIMENSION(:) :: avecolvar => NULL()
53 : REAL(KIND=dp), POINTER, DIMENSION(:) :: aveMmatrix => NULL()
54 : END TYPE average_quantities_type
55 :
56 : ! **************************************************************************************************
57 : INTERFACE get_averages
58 : MODULE PROCEDURE get_averages_rs, get_averages_rv, get_averages_rm
59 : END INTERFACE get_averages
60 :
61 : ! *** Public subroutines and data types ***
62 : PUBLIC :: average_quantities_type, create_averages, release_averages, &
63 : retain_averages, compute_averages
64 :
65 : ! *** Global parameters ***
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'averages_types'
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief Creates averages environment
72 : !> \param averages ...
73 : !> \param averages_section ...
74 : !> \param virial_avg ...
75 : !> \param force_env ...
76 : !> \author Teodoro Laino [tlaino] - 03.2008 - University of Zurich
77 : ! **************************************************************************************************
78 1789 : SUBROUTINE create_averages(averages, averages_section, virial_avg, force_env)
79 : TYPE(average_quantities_type), POINTER :: averages
80 : TYPE(section_vals_type), POINTER :: averages_section
81 : LOGICAL, INTENT(IN), OPTIONAL :: virial_avg
82 : TYPE(force_env_type), POINTER :: force_env
83 :
84 : INTEGER :: i, nint
85 : LOGICAL :: do_colvars
86 :
87 1789 : ALLOCATE (averages)
88 : ! Point to the averages section
89 1789 : averages%averages_section => averages_section
90 1789 : averages%ref_count = 1
91 1789 : CALL section_vals_val_get(averages_section, "_SECTION_PARAMETERS_", l_val=averages%do_averages)
92 1789 : IF (averages%do_averages) THEN
93 : ! Setup Virial if requested
94 1787 : IF (PRESENT(virial_avg)) THEN
95 20 : IF (virial_avg) THEN
96 4960 : ALLOCATE (averages%virial)
97 : END IF
98 : END IF
99 1787 : CALL section_vals_val_get(averages_section, "AVERAGE_COLVAR", l_val=do_colvars)
100 : ! Total number of COLVARs
101 1787 : nint = 0
102 1787 : IF (do_colvars) nint = number_of_colvar(force_env)
103 3576 : ALLOCATE (averages%avecolvar(nint))
104 3576 : ALLOCATE (averages%aveMmatrix(nint*nint))
105 1883 : DO i = 1, nint
106 1883 : averages%avecolvar(i) = 0.0_dp
107 : END DO
108 8182 : DO i = 1, nint*nint
109 6395 : averages%aveMmatrix(i) = 0.0_dp
110 : END DO
111 : END IF
112 1789 : END SUBROUTINE create_averages
113 :
114 : ! **************************************************************************************************
115 : !> \brief retains the given averages env
116 : !> \param averages ...
117 : !> \author Teodoro Laino [tlaino] - 03.2008 - University of Zurich
118 : ! **************************************************************************************************
119 20 : SUBROUTINE retain_averages(averages)
120 : TYPE(average_quantities_type), POINTER :: averages
121 :
122 20 : CPASSERT(ASSOCIATED(averages))
123 20 : CPASSERT(averages%ref_count > 0)
124 20 : averages%ref_count = averages%ref_count + 1
125 20 : END SUBROUTINE retain_averages
126 :
127 : ! **************************************************************************************************
128 : !> \brief releases the given averages env
129 : !> \param averages ...
130 : !> \author Teodoro Laino [tlaino] - 03.2008 - University of Zurich
131 : ! **************************************************************************************************
132 1809 : SUBROUTINE release_averages(averages)
133 : TYPE(average_quantities_type), POINTER :: averages
134 :
135 : TYPE(section_vals_type), POINTER :: work_section
136 :
137 1809 : IF (ASSOCIATED(averages)) THEN
138 1809 : CPASSERT(averages%ref_count > 0)
139 1809 : averages%ref_count = averages%ref_count - 1
140 1809 : IF (averages%ref_count == 0) THEN
141 1789 : IF (ASSOCIATED(averages%virial)) DEALLOCATE (averages%virial)
142 1789 : IF (ASSOCIATED(averages%avecolvar)) THEN
143 1787 : DEALLOCATE (averages%avecolvar)
144 : END IF
145 1789 : IF (ASSOCIATED(averages%aveMmatrix)) THEN
146 1787 : DEALLOCATE (averages%aveMmatrix)
147 : END IF
148 : ! Removes restart values from the corresponding restart section..
149 1789 : work_section => section_vals_get_subs_vals(averages%averages_section, "RESTART_AVERAGES")
150 1789 : CALL section_vals_remove_values(work_section)
151 1789 : NULLIFY (averages%averages_section)
152 : !
153 1789 : DEALLOCATE (averages)
154 : END IF
155 : END IF
156 :
157 1809 : END SUBROUTINE release_averages
158 :
159 : ! **************************************************************************************************
160 : !> \brief computes the averages
161 : !> \param averages ...
162 : !> \param force_env ...
163 : !> \param md_ener ...
164 : !> \param cell ...
165 : !> \param virial ...
166 : !> \param pv_scalar ...
167 : !> \param pv_xx ...
168 : !> \param used_time ...
169 : !> \param hugoniot ...
170 : !> \param abc ...
171 : !> \param cell_angle ...
172 : !> \param nat ...
173 : !> \param itimes ...
174 : !> \param time ...
175 : !> \param my_pos ...
176 : !> \param my_act ...
177 : !> \author Teodoro Laino [tlaino] - 03.2008 - University of Zurich
178 : ! **************************************************************************************************
179 81178 : SUBROUTINE compute_averages(averages, force_env, md_ener, cell, virial, &
180 : pv_scalar, pv_xx, used_time, hugoniot, abc, cell_angle, nat, itimes, &
181 : time, my_pos, my_act)
182 : TYPE(average_quantities_type), POINTER :: averages
183 : TYPE(force_env_type), POINTER :: force_env
184 : TYPE(md_ener_type), POINTER :: md_ener
185 : TYPE(cell_type), POINTER :: cell
186 : TYPE(virial_type), POINTER :: virial
187 : REAL(KIND=dp), INTENT(IN) :: pv_scalar, pv_xx
188 : REAL(KIND=dp), POINTER :: used_time
189 : REAL(KIND=dp), INTENT(IN) :: hugoniot
190 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: abc, cell_angle
191 : INTEGER, INTENT(IN) :: nat, itimes
192 : REAL(KIND=dp), INTENT(IN) :: time
193 : CHARACTER(LEN=default_string_length), INTENT(IN) :: my_pos, my_act
194 :
195 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_averages'
196 :
197 : CHARACTER(LEN=default_string_length) :: ctmp
198 : INTEGER :: delta_t, handle, i, nint, output_unit
199 : LOGICAL :: restart_averages
200 : REAL(KIND=dp) :: start_time
201 40589 : REAL(KIND=dp), DIMENSION(:), POINTER :: cvalues, Mmatrix, tmp
202 : TYPE(cp_logger_type), POINTER :: logger
203 : TYPE(section_vals_type), POINTER :: restart_section
204 :
205 40589 : CALL timeset(routineN, handle)
206 : CALL section_vals_val_get(averages%averages_section, "ACQUISITION_START_TIME", &
207 40589 : r_val=start_time)
208 40589 : IF (averages%do_averages) THEN
209 40569 : NULLIFY (cvalues, Mmatrix, logger)
210 40569 : logger => cp_get_default_logger()
211 : ! Determine the nr. of internal colvar (if any/requested)
212 40569 : nint = 0
213 40569 : IF (ASSOCIATED(averages%avecolvar)) nint = SIZE(averages%avecolvar)
214 :
215 : ! Evaluate averages if we collected enough statistics (user defined)
216 40569 : IF (time >= start_time) THEN
217 :
218 : ! Handling properly the restart
219 40555 : IF (averages%itimes_start == -1) THEN
220 1753 : restart_section => section_vals_get_subs_vals(averages%averages_section, "RESTART_AVERAGES")
221 1753 : CALL section_vals_get(restart_section, explicit=restart_averages)
222 1753 : IF (restart_averages) THEN
223 172 : CALL section_vals_val_get(restart_section, "ITIMES_START", i_val=averages%itimes_start)
224 172 : CALL section_vals_val_get(restart_section, "AVECPU", r_val=averages%avecpu)
225 172 : CALL section_vals_val_get(restart_section, "AVEHUGONIOT", r_val=averages%avehugoniot)
226 172 : CALL section_vals_val_get(restart_section, "AVETEMP_BARO", r_val=averages%avetemp_baro)
227 172 : CALL section_vals_val_get(restart_section, "AVEPOT", r_val=averages%avepot)
228 172 : CALL section_vals_val_get(restart_section, "AVEKIN", r_val=averages%avekin)
229 172 : CALL section_vals_val_get(restart_section, "AVETEMP", r_val=averages%avetemp)
230 172 : CALL section_vals_val_get(restart_section, "AVEKIN_QM", r_val=averages%avekin_qm)
231 172 : CALL section_vals_val_get(restart_section, "AVETEMP_QM", r_val=averages%avetemp_qm)
232 172 : CALL section_vals_val_get(restart_section, "AVEVOL", r_val=averages%avevol)
233 172 : CALL section_vals_val_get(restart_section, "AVECELL_A", r_val=averages%aveca)
234 172 : CALL section_vals_val_get(restart_section, "AVECELL_B", r_val=averages%avecb)
235 172 : CALL section_vals_val_get(restart_section, "AVECELL_C", r_val=averages%avecc)
236 172 : CALL section_vals_val_get(restart_section, "AVEALPHA", r_val=averages%aveal)
237 172 : CALL section_vals_val_get(restart_section, "AVEBETA", r_val=averages%avebe)
238 172 : CALL section_vals_val_get(restart_section, "AVEGAMMA", r_val=averages%avega)
239 172 : CALL section_vals_val_get(restart_section, "AVE_ECONS", r_val=averages%econs)
240 : ! Virial
241 172 : IF (virial%pv_availability) THEN
242 76 : CALL section_vals_val_get(restart_section, "AVE_PRESS", r_val=averages%avepress)
243 76 : CALL section_vals_val_get(restart_section, "AVE_PXX", r_val=averages%avepxx)
244 76 : IF (ASSOCIATED(averages%virial)) THEN
245 0 : CALL section_vals_val_get(restart_section, "AVE_PV_TOT", r_vals=tmp)
246 0 : averages%virial%pv_total = RESHAPE(tmp, [3, 3])
247 0 : CALL section_vals_val_get(restart_section, "AVE_PV_VIR", r_vals=tmp)
248 0 : averages%virial%pv_virial = RESHAPE(tmp, [3, 3])
249 0 : CALL section_vals_val_get(restart_section, "AVE_PV_KIN", r_vals=tmp)
250 0 : averages%virial%pv_kinetic = RESHAPE(tmp, [3, 3])
251 0 : CALL section_vals_val_get(restart_section, "AVE_PV_CNSTR", r_vals=tmp)
252 0 : averages%virial%pv_constraint = RESHAPE(tmp, [3, 3])
253 0 : CALL section_vals_val_get(restart_section, "AVE_PV_XC", r_vals=tmp)
254 0 : averages%virial%pv_xc = RESHAPE(tmp, [3, 3])
255 0 : CALL section_vals_val_get(restart_section, "AVE_PV_FOCK_4C", r_vals=tmp)
256 0 : averages%virial%pv_fock_4c = RESHAPE(tmp, [3, 3])
257 : END IF
258 : END IF
259 : ! Colvars
260 172 : IF (nint > 0) THEN
261 0 : CALL section_vals_val_get(restart_section, "AVE_COLVARS", r_vals=cvalues)
262 0 : CALL section_vals_val_get(restart_section, "AVE_MMATRIX", r_vals=Mmatrix)
263 0 : CPASSERT(nint == SIZE(cvalues))
264 0 : CPASSERT(nint*nint == SIZE(Mmatrix))
265 0 : averages%avecolvar = cvalues
266 0 : averages%aveMmatrix = Mmatrix
267 : END IF
268 : ELSE
269 1581 : averages%itimes_start = itimes
270 : END IF
271 : END IF
272 40555 : delta_t = itimes - averages%itimes_start + 1
273 :
274 : ! Perform averages
275 1587 : SELECT CASE (delta_t)
276 : CASE (1)
277 1587 : averages%avecpu = used_time
278 1587 : averages%avehugoniot = hugoniot
279 1587 : averages%avetemp_baro = md_ener%temp_baro
280 1587 : averages%avepot = md_ener%epot
281 1587 : averages%avekin = md_ener%ekin
282 1587 : averages%avetemp = md_ener%temp_part
283 1587 : averages%avekin_qm = md_ener%ekin_qm
284 1587 : averages%avetemp_qm = md_ener%temp_qm
285 1587 : averages%avevol = cell%deth
286 1587 : averages%aveca = abc(1)
287 1587 : averages%avecb = abc(2)
288 1587 : averages%avecc = abc(3)
289 1587 : averages%aveal = cell_angle(3)
290 1587 : averages%avebe = cell_angle(2)
291 1587 : averages%avega = cell_angle(1)
292 1587 : averages%econs = 0._dp
293 : ! Virial
294 1587 : IF (virial%pv_availability) THEN
295 224 : averages%avepress = pv_scalar
296 224 : averages%avepxx = pv_xx
297 224 : IF (ASSOCIATED(averages%virial)) THEN
298 520 : averages%virial%pv_total = virial%pv_total
299 520 : averages%virial%pv_virial = virial%pv_virial
300 520 : averages%virial%pv_kinetic = virial%pv_kinetic
301 520 : averages%virial%pv_constraint = virial%pv_constraint
302 520 : averages%virial%pv_xc = virial%pv_xc
303 520 : averages%virial%pv_fock_4c = virial%pv_fock_4c
304 : END IF
305 : END IF
306 : ! Colvars
307 1587 : IF (nint > 0) THEN
308 : CALL get_clv_force(force_env, nsize_xyz=nat*3, nsize_int=nint, &
309 2 : cvalues=averages%avecolvar, Mmatrix=averages%aveMmatrix)
310 : END IF
311 : CASE DEFAULT
312 38968 : CALL get_averages(averages%avecpu, used_time, delta_t)
313 38968 : CALL get_averages(averages%avehugoniot, hugoniot, delta_t)
314 38968 : CALL get_averages(averages%avetemp_baro, md_ener%temp_baro, delta_t)
315 38968 : CALL get_averages(averages%avepot, md_ener%epot, delta_t)
316 38968 : CALL get_averages(averages%avekin, md_ener%ekin, delta_t)
317 38968 : CALL get_averages(averages%avetemp, md_ener%temp_part, delta_t)
318 38968 : CALL get_averages(averages%avekin_qm, md_ener%ekin_qm, delta_t)
319 38968 : CALL get_averages(averages%avetemp_qm, md_ener%temp_qm, delta_t)
320 38968 : CALL get_averages(averages%avevol, cell%deth, delta_t)
321 38968 : CALL get_averages(averages%aveca, abc(1), delta_t)
322 38968 : CALL get_averages(averages%avecb, abc(2), delta_t)
323 38968 : CALL get_averages(averages%avecc, abc(3), delta_t)
324 38968 : CALL get_averages(averages%aveal, cell_angle(3), delta_t)
325 38968 : CALL get_averages(averages%avebe, cell_angle(2), delta_t)
326 38968 : CALL get_averages(averages%avega, cell_angle(1), delta_t)
327 38968 : CALL get_averages(averages%econs, md_ener%delta_cons, delta_t)
328 : ! Virial
329 38968 : IF (virial%pv_availability) THEN
330 4220 : CALL get_averages(averages%avepress, pv_scalar, delta_t)
331 4220 : CALL get_averages(averages%avepxx, pv_xx, delta_t)
332 4220 : IF (ASSOCIATED(averages%virial)) THEN
333 180 : CALL get_averages(averages%virial%pv_total, virial%pv_total, delta_t)
334 180 : CALL get_averages(averages%virial%pv_virial, virial%pv_virial, delta_t)
335 180 : CALL get_averages(averages%virial%pv_kinetic, virial%pv_kinetic, delta_t)
336 180 : CALL get_averages(averages%virial%pv_constraint, virial%pv_constraint, delta_t)
337 180 : CALL get_averages(averages%virial%pv_xc, virial%pv_xc, delta_t)
338 180 : CALL get_averages(averages%virial%pv_fock_4c, virial%pv_fock_4c, delta_t)
339 : END IF
340 : END IF
341 : ! Colvars
342 79523 : IF (nint > 0) THEN
343 72 : ALLOCATE (cvalues(nint))
344 72 : ALLOCATE (Mmatrix(nint*nint))
345 : CALL get_clv_force(force_env, nsize_xyz=nat*3, nsize_int=nint, cvalues=cvalues, &
346 24 : Mmatrix=Mmatrix)
347 24 : CALL get_averages(averages%avecolvar, cvalues, delta_t)
348 24 : CALL get_averages(averages%aveMmatrix, Mmatrix, delta_t)
349 24 : DEALLOCATE (cvalues)
350 24 : DEALLOCATE (Mmatrix)
351 : END IF
352 : END SELECT
353 : END IF
354 :
355 : ! Possibly print averages
356 : output_unit = cp_print_key_unit_nr(logger, averages%averages_section, "PRINT_AVERAGES", &
357 40569 : extension=".avg", file_position=my_pos, file_action=my_act)
358 40569 : IF (output_unit > 0) THEN
359 : WRITE (output_unit, FMT='(A15,1X,"=",1X,G15.9," NSTEP #",I15)') &
360 4 : "AVECPU", averages%avecpu, itimes, &
361 4 : "AVEHUGONIOT", averages%avehugoniot, itimes, &
362 4 : "AVETEMP_BARO", averages%avetemp_baro, itimes, &
363 4 : "AVEPOT", averages%avepot, itimes, &
364 4 : "AVEKIN", averages%avekin, itimes, &
365 4 : "AVETEMP", averages%avetemp, itimes, &
366 4 : "AVEKIN_QM", averages%avekin_qm, itimes, &
367 4 : "AVETEMP_QM", averages%avetemp_qm, itimes, &
368 4 : "AVEVOL", averages%avevol, itimes, &
369 4 : "AVECELL_A", averages%aveca, itimes, &
370 4 : "AVECELL_B", averages%avecb, itimes, &
371 4 : "AVECELL_C", averages%avecc, itimes, &
372 4 : "AVEALPHA", averages%aveal, itimes, &
373 4 : "AVEBETA", averages%avebe, itimes, &
374 4 : "AVEGAMMA", averages%avega, itimes, &
375 8 : "AVE_ECONS", averages%econs, itimes
376 : ! Print the virial
377 4 : IF (virial%pv_availability) THEN
378 : WRITE (output_unit, FMT='(A15,1X,"=",1X,G15.9," NSTEP #",I15)') &
379 0 : "AVE_PRESS", averages%avepress, itimes, &
380 0 : "AVE_PXX", averages%avepxx, itimes
381 0 : IF (ASSOCIATED(averages%virial)) THEN
382 : WRITE (output_unit, FMT='(A15,1X,"=",1X,G15.9," NSTEP #",I15)') &
383 0 : "AVE_PV_TOT", averages%virial%pv_total, itimes, &
384 0 : "AVE_PV_VIR", averages%virial%pv_virial, itimes, &
385 0 : "AVE_PV_KIN", averages%virial%pv_kinetic, itimes, &
386 0 : "AVE_PV_CNSTR", averages%virial%pv_constraint, itimes, &
387 0 : "AVE_PV_XC", averages%virial%pv_xc, itimes, &
388 0 : "AVE_PV_FOCK_4C", averages%virial%pv_fock_4c, itimes
389 : END IF
390 : END IF
391 196 : DO i = 1, nint
392 192 : ctmp = cp_to_string(i)
393 : WRITE (output_unit, FMT='(A15,1X,"=",1X,G15.9," NSTEP #",I15)') &
394 196 : TRIM("AVE_CV-"//ADJUSTL(ctmp)), averages%avecolvar(i), itimes
395 : END DO
396 4 : WRITE (output_unit, FMT='(/)')
397 : END IF
398 : CALL cp_print_key_finished_output(output_unit, logger, averages%averages_section, &
399 40569 : "PRINT_AVERAGES")
400 : END IF
401 40589 : CALL timestop(handle)
402 40589 : END SUBROUTINE compute_averages
403 :
404 : ! **************************************************************************************************
405 : !> \brief computes the averages - low level for REAL
406 : !> \param avg ...
407 : !> \param add ...
408 : !> \param delta_t ...
409 : !> \author Teodoro Laino [tlaino] - 03.2008 - University of Zurich
410 : ! **************************************************************************************************
411 631928 : SUBROUTINE get_averages_rs(avg, add, delta_t)
412 : REAL(KIND=dp), INTENT(INOUT) :: avg
413 : REAL(KIND=dp), INTENT(IN) :: add
414 : INTEGER, INTENT(IN) :: delta_t
415 :
416 631928 : avg = (avg*REAL(delta_t - 1, dp) + add)/REAL(delta_t, dp)
417 631928 : END SUBROUTINE get_averages_rs
418 :
419 : ! **************************************************************************************************
420 : !> \brief computes the averages - low level for REAL vector
421 : !> \param avg ...
422 : !> \param add ...
423 : !> \param delta_t ...
424 : !> \author Teodoro Laino [tlaino] - 10.2008 - University of Zurich
425 : ! **************************************************************************************************
426 48 : SUBROUTINE get_averages_rv(avg, add, delta_t)
427 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: avg
428 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: add
429 : INTEGER, INTENT(IN) :: delta_t
430 :
431 : INTEGER :: i
432 : LOGICAL :: check
433 :
434 48 : check = SIZE(avg) == SIZE(add)
435 48 : CPASSERT(check)
436 56496 : DO i = 1, SIZE(avg)
437 56496 : avg(i) = (avg(i)*REAL(delta_t - 1, dp) + add(i))/REAL(delta_t, dp)
438 : END DO
439 48 : END SUBROUTINE get_averages_rv
440 :
441 : ! **************************************************************************************************
442 : !> \brief computes the averages - low level for REAL matrix
443 : !> \param avg ...
444 : !> \param add ...
445 : !> \param delta_t ...
446 : !> \author Teodoro Laino [tlaino] - 10.2008 - University of Zurich
447 : ! **************************************************************************************************
448 1080 : SUBROUTINE get_averages_rm(avg, add, delta_t)
449 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: avg
450 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: add
451 : INTEGER, INTENT(IN) :: delta_t
452 :
453 : INTEGER :: i, j
454 : LOGICAL :: check
455 :
456 1080 : check = SIZE(avg, 1) == SIZE(add, 1)
457 1080 : CPASSERT(check)
458 1080 : check = SIZE(avg, 2) == SIZE(add, 2)
459 1080 : CPASSERT(check)
460 4320 : DO i = 1, SIZE(avg, 2)
461 14040 : DO j = 1, SIZE(avg, 1)
462 12960 : avg(j, i) = (avg(j, i)*REAL(delta_t - 1, dp) + add(j, i))/REAL(delta_t, dp)
463 : END DO
464 : END DO
465 1080 : END SUBROUTINE get_averages_rm
466 :
467 0 : END MODULE averages_types
|
