PLUME-MoM  1.0
Integralvolcanicplumemodel
 All Classes Files Functions Variables Pages
rise.f90
Go to the documentation of this file.
1 !********************************************************************************
3 !
9 !********************************************************************************
10 MODULE rise
11 
12  IMPLICIT NONE
13  SAVE
14 
15 CONTAINS
16 
17  !******************************************************************************
19  !
26  !******************************************************************************
27 
28  SUBROUTINE plumerise
29 
30  ! external variables
31  USE meteo_module, ONLY : rho_atm
32  USE mixture_module, ONLY : gas_mass_fraction , rho_mix, mass_flow_rate
33  USE particles_module, ONLY : n_part , mom0 , mom
34  USE particles_module, ONLY : distribution_variable
35  USE particles_module, ONLY : solid_partial_mass_fraction
36  USE plume_module, ONLY: s , w , z , vent_height , r , mag_u
37  USE solver_module, ONLY: ds, ds0, f, ftemp, rhs, rhstemp
38  USE solver_module, ONLY: f_stepold
39  USE variables, ONLY : verbose_level
40  USE variables, ONLY : dakota_flag , hysplit_flag
41  USE variables, ONLY : pi_g
42  USE variables, ONLY : height_weight , height_obj , mu_weight , &
43  mu_obj , sigma_weight , sigma_obj , skew_weight , skew_obj
44 
45  ! external procedures
46  USE inpout, ONLY: write_column , write_dakota , write_hysplit
47  USE meteo_module, ONLY: zmet, initialize_meteo
48  USE mixture_module, ONLY: initialize_mixture
49  USE particles_module, ONLY: initialize_particles
50  USE plume_module, ONLY: initialize_plume
51  USE solver_module, ONLY: rate, lump, marching, unlump
52 
53 
54  IMPLICIT NONE
55 
56  CHARACTER(LEN=20) :: description
57 
58  CHARACTER(len=8) :: x1 ! format descriptor
59 
60  INTEGER :: i_part
61 
62  REAL*8 :: mu(4)
63 
64  REAL*8 :: k1 , k2
65 
66  REAL*8 :: mu_phi , sigma_phi , skew_phi
67 
68  REAL*8 :: mass_fract
69 
70  REAL*8 :: solid_mass_flux , solid_mass_flux0
71 
72  REAL*8 :: solid_mass_flux_change
73 
74  REAL*8 :: obj_function
75 
76  REAL*8 :: w_old , w_oldold
77  REAL*8 :: w_minrel , w_maxrel
78  REAL*8 :: w_maxabs
79 
80  REAL*8 :: check_sb
81  REAL*8 :: eps_sb
82 
83  REAL*8 :: z_array(10000000)
84  REAL*8 :: w_array(10000000)
85 
86  REAL*8, ALLOCATABLE :: z_norm(:) , w_norm(:)
87  REAL*8, ALLOCATABLE :: first_der_right(:) , first_der_left(:)
88  REAL*8, ALLOCATABLE :: sec_der(:) , k(:)
89  REAL*8, ALLOCATABLE :: first_der_central(:)
90 
91  REAL*8 :: k_max
92 
93  REAL*8 :: column_regime
94 
95  INTEGER :: idx , max_idx
96 
97  REAL*8 :: delta_rho
98 
99  REAL*8 :: plume_height
100  REAL*8 :: height_nbl
101  REAL*8 :: deltarho_min
102  REAL*8 :: rho_nbl
103 
104  REAL*8 :: deltarho , deltarho_old
105 
106  !
107  ! ... Set initial conditions at the release height
108  !
109  CALL initialize_plume
110 
111  CALL initialize_meteo
112 
113  CALL initialize_particles
114 
115  !
116  ! ... Get meteo variables at release height
117  !
118  CALL zmet
119 
120  CALL initialize_mixture
121 
122  description = 'Initial MFR'
123 
124  CALL write_dakota(description,mass_flow_rate)
125 
126  w_old = w
127  w_oldold = w
128 
129  w_maxabs = w
130 
131  w_minrel = w
132  w_maxrel = w
133 
134  idx = 1
135 
136  z_array(idx) = z
137  w_array(idx) = w
138 
139  delta_rho = rho_mix - rho_atm
140 
141  DO i_part=1,n_part
142 
143  WRITE(x1,'(I2.2)') i_part ! converting integer to string using a 'internal file'
144 
145  description = 'Mean Diameter '//trim(x1)
146 
147  CALL write_dakota(description,mom0(i_part,1)/mom0(i_part,0))
148 
149  description = 'Sau. Mean Diam. '//trim(x1)
150 
151  CALL write_dakota(description,mom0(i_part,3)/mom0(i_part,2))
152 
153  !------ Convert from raw moments to central moments ----------------------
154  !------ http://mathworld.wolfram.com/CentralMoment.html ------------------
155 
156  mu(1) = mom0(i_part,1)/mom0(i_part,0)
157 
158  mu(2) = -(mom0(i_part,1)/mom0(i_part,0))**2 + (mom0(i_part,2) / &
159  mom0(i_part,0))
160 
161  mu(3) = 2.d0 * (mom0(i_part,1)/mom0(i_part,0))**3 &
162  - 3.d0 * (mom0(i_part,1)/mom0(i_part,0)) * (mom0(i_part,2) / &
163  mom0(i_part,0)) &
164  + (mom0(i_part,3)/mom0(i_part,0))
165 
166  mu(4) = -3.d0 * (mom0(i_part,1)/mom0(i_part,0))**4 &
167  + 6.d0 * (mom0(i_part,1)/mom0(i_part,0))**2 * (mom0(i_part,2) / &
168  mom0(i_part,0)) - 4.d0 * (mom0(i_part,1)/mom0(i_part,0)) * &
169  (mom0(i_part,3)/mom0(i_part,0)) + (mom0(i_part,4)/mom0(i_part,0))
170 
171  ! WRITE(*,*) 'mu',mu(1:4)
172  ! READ(*,*)
173 
174  ! description = 'Variance '//trim(x1)
175  !
176  ! CALL write_dakota(description,mu(2))
177  !
178  ! description = 'Skewness '//trim(x1)
179  !
180  ! CALL write_dakota( description , mu(3) / mu(2)**(3.D0/2.D0) )
181  !
182  ! description = 'Kurtosis '//trim(x1)
183  !
184  ! CALL write_dakota( description , mu(4)/ ( mu(2)**2 ) - 3.D0 )
185 
186  END DO
187 
188  DO i_part=1,n_part
189 
190  WRITE(x1,'(I2.2)') i_part ! converting integer to string using a 'internal file'
191 
192  IF ( distribution_variable .EQ. 'particles_number' ) THEN
193 
194  k1 = log( 1.d3 * mom0(i_part,1) / mom0(i_part,0) )
195  k2 = log( 1.d3 * mom0(i_part,3) / mom0(i_part,2) )
196 
197  mu_phi = - 0.25d0 * ( 5*k2 - k1 ) / log(2.d0)
198 
199  sigma_phi = sqrt( 0.5d0 * (k2-k1) ) / log(2.d0)
200 
201  ELSEIF ( distribution_variable .EQ. 'mass_fraction' ) THEN
202 
203  mu_phi = mom0(i_part,1)/mom0(i_part,0)
204 
205  sigma_phi = sqrt( -(mom0(i_part,1)/mom0(i_part,0))**2 + &
206  (mom0(i_part,2)/mom0(i_part,0)) )
207 
208  skew_phi = ( 2.d0 * (mom0(i_part,1)/mom0(i_part,0))**3 &
209  - 3.d0 * mu_phi * ( mom0(i_part,2) / mom0(i_part,0) ) &
210  + mom0(i_part,3) / mom0(i_part,0) ) / sigma_phi**3
211 
212  END IF
213 
214  description = 'Init Avg Diam '//trim(x1)
215 
216  CALL write_dakota(description,mu_phi)
217 
218  description = 'Init Var Diam '//trim(x1)
219 
220  CALL write_dakota(description,sigma_phi)
221 
222  description = 'Init Skw Diam '//trim(x1)
223 
224  CALL write_dakota(description,skew_phi)
225 
226  mass_fract = solid_partial_mass_fraction(i_part) * ( 1.d0 - &
227  gas_mass_fraction)
228 
229  description = 'Init Mass Fract '//trim(x1)
230 
231  CALL write_dakota(description,mass_fract)
232 
233  solid_mass_flux0 = solid_partial_mass_fraction(i_part) * ( 1.d0 - &
234  gas_mass_fraction) * rho_mix * pi_g * r**2 * mag_u
235 
236  description = 'Init Solid Flux '//trim(x1)
237 
238  CALL write_dakota(description,solid_mass_flux0)
239 
240  END DO
241 
242  !
243  ! ... Lump physical variables
244  !
245  CALL lump(f)
246  !
247  ! ----------------------------------------------------
248  ! ... assign initial stepping length and
249  ! ... start plume rise marching loop
250  ! ----------------------------------------------------
251  !
252  ds = ds0
253 
254  IF ( .NOT.dakota_flag ) CALL write_column
255 
256  IF ( hysplit_flag ) CALL write_hysplit(.false.)
257 
258  deltarho_min = 1000.d0
259 
260  deltarho_old = 0.d0
261 
262  main_loop: DO
263 
264  f_stepold = f
265 
266  !
267  ! ---------------------------------------------------------
268  ! ... compute meteo conditions and the rhs of the equations
269  ! ... for this location
270  ! ---------------------------------------------------------
271  !
272  IF ( verbose_level .GE. 2 ) THEN
273 
274  WRITE(*,*)
275  WRITE(*,*) '**************** BEFORE PREDICTOR STEP *****************'
276  WRITE(*,*)
277 
278  END IF
279 
280  CALL unlump(f)
281 
282  w_oldold = w_old
283  w_old = w
284 
285  CALL rate(rhs)
286  !
287  ! ... predictor step (compute temporary quantities)
288  CALL marching(f,ftemp,rhs)
289 
290  IF ( verbose_level .GE. 2 ) THEN
291 
292  WRITE(*,*)
293  WRITE(*,*) '**************** BEFORE CORRECTOR STEP *****************'
294  WRITE(*,*)
295 
296  END IF
297 
298  CALL unlump(ftemp)
299  CALL rate(rhstemp)
300  !
301  ! ... corrector step
302  rhs = 0.5d0 * ( rhstemp - rhs )
303  CALL marching(ftemp,f,rhs)
304 
305  IF ( verbose_level .GE. 2 ) THEN
306 
307  WRITE(*,*)
308  WRITE(*,*) '**************** AFTER CORRECTOR STEP *****************'
309  WRITE(*,*)
310 
311  END IF
312 
313 
314  CALL unlump(f)
315 
316  !
317  ! Reduce step-size condition
318  !
319  IF ( w .LE. 0.d0) THEN
320 
321  ds = 0.5d0 * ds
322  f = f_stepold
323 
324  ELSE
325 
326  idx = idx + 1
327 
328  z_array(idx) = z
329  w_array(idx) = w
330 
331  IF ( ( w_old .LT. w ) .AND. ( w_old .LT. w_oldold ) ) THEN
332 
333  w_minrel = w_old
334 
335  !WRITE(*,*) 'minrel',w_oldold,w_old,w
336 
337  END IF
338 
339  IF ( w .GT. w_maxabs ) w_maxabs = w
340 
341  IF ( ( w_old .GT. w ) .AND. ( w_old .GT. w_oldold ) ) THEN
342 
343  w_maxrel = w_old
344 
345  !WRITE(*,*) 'maxrel',w_oldold,w_old,w
346 
347  END IF
348 
349  delta_rho = min( delta_rho , rho_mix - rho_atm )
350 
351  ! used to define the neutral buoyancy level
352  deltarho = rho_mix - rho_atm
353 
354 
355  IF ( deltarho * deltarho_old .LT. 0.d0 ) THEN
356 
357  rho_nbl = rho_mix
358  height_nbl = z - vent_height
359 
360  END IF
361 
362  s = s + ds
363 
364  deltarho_old = deltarho
365 
366  IF ( verbose_level .GE. 2 ) THEN
367 
368  DO i_part=1,n_part
369 
370  WRITE(*,*) '**',mom(i_part,1)/mom(i_part,0)
371 
372  END DO
373 
374  READ(*,*)
375 
376  END IF
377 
378  IF ( .NOT.dakota_flag ) CALL write_column
379 
380  IF ( hysplit_flag ) CALL write_hysplit(.false.)
381 
382  ! ... Exit condition
383 
384  IF ( w .LE. 1.d-5 ) THEN
385 
386  EXIT main_loop
387 
388  END IF
389 
390  END IF
391 
392  END DO main_loop
393 
394  IF ( hysplit_flag ) CALL write_hysplit(.true.)
395 
396 
397  max_idx = idx
398 
399  ALLOCATE( z_norm(max_idx) , w_norm(idx) )
400  ALLOCATE( first_der_right(max_idx-1) , first_der_left(max_idx-1) )
401  ALLOCATE( sec_der(max_idx-2) , k(max_idx-2) )
402  ALLOCATE( first_der_central(max_idx-2) )
403 
404  z_norm = ( z_array(1:max_idx) - minval( z_array(1:max_idx) ) ) &
405  / ( maxval( z_array(1:max_idx) ) - minval( z_array(1:max_idx) ) )
406 
407  w_norm = w_array(1:max_idx) / maxval( w_array(1:max_idx) )
408 
409  first_der_right(1:max_idx-2) = ( z_norm(3:max_idx) - z_norm(2:max_idx-1) ) &
410  / ( w_norm(3:max_idx) - w_norm(2:max_idx-1) )
411 
412  first_der_left(1:max_idx-2) = ( z_norm(2:max_idx-1) - z_norm(1:max_idx-2) ) &
413  / ( w_norm(2:max_idx-1) - w_norm(1:max_idx-2) )
414 
415  sec_der(1:max_idx-2) = ( first_der_right(1:max_idx-2) - first_der_left(1:max_idx-2) ) &
416  / ( 0.5 * ( w_norm(3:max_idx) - w_norm(1:max_idx-2 ) ) )
417 
418  first_der_central = ( z_norm(3:max_idx) - z_norm(1:max_idx-2) ) &
419  / ( w_norm(3:max_idx) - w_norm(1:max_idx-2) )
420 
421  k(1:max_idx-2) = sec_der(1:max_idx-2) / ( ( 1.0 + first_der_central(1:max_idx-2)**2 )**(1.5d0) )
422 
423  k_max = maxval( k(1:max_idx-2) )
424 
425  check_sb = ( w_maxrel - w_minrel ) / w_maxabs
426 
427  eps_sb = 0.05
428 
429 !!$ IF ( check_sb .GT. eps_sb ) THEN
430 !!$
431 !!$ WRITE(*,*) 'w_minrel,w_maxrel,w_maxabs',w_minrel,w_maxrel,w_maxabs
432 !!$
433 !!$ WRITE(*,*) 'superbuoyant'
434 !!$
435 !!$ column_regime = 2
436 !!$
437 !!$ ELSE
438 !!$
439 !!$ WRITE(*,*) 'k_max',k_max
440 !!$
441 !!$ IF ( k_max < 1.d0 ) THEN
442 !!$
443 !!$ WRITE(*,*) 'collapsing'
444 !!$
445 !!$ column_regime = 3
446 !!$
447 !!$ ELSE
448 !!$
449 !!$ WRITE(*,*) 'buoyant'
450 !!$
451 !!$ column_regime = 1
452 !!$
453 !!$ END IF
454 
455 
456  ! modified criterium for collapsing
457 
458  IF ( delta_rho .GT. 0.d0 ) THEN
459 
460  WRITE(*,*) 'collapsing'
461 
462  column_regime = 3
463 
464  ELSE
465 
466  IF ( check_sb .GT. eps_sb ) THEN
467 
468  !WRITE(*,*) 'w_minrel,w_maxrel,w_maxabs',w_minrel,w_maxrel,w_maxabs
469 
470  WRITE(*,*) 'superbuoyant'
471 
472  column_regime = 2
473 
474  ELSE
475 
476  WRITE(*,*) 'buoyant'
477 
478  column_regime = 1
479 
480  END IF
481 
482  END IF
483 
484 
485  description = 'Column regime'
486 
487  CALL write_dakota(description,column_regime)
488 
489  description = 'NBL height (atv)'
490 
491  CALL write_dakota(description,height_nbl)
492 
493  description = 'Plume height (atv)'
494  plume_height = z - vent_height
495 
496  CALL write_dakota(description,plume_height)
497 
498  DO i_part=1,n_part
499 
500  WRITE(x1,'(I2.2)') i_part ! convert int to string using an 'internal file'
501 
502  IF ( distribution_variable .EQ. 'particles_number' ) THEN
503 
504  k1 = log( 1.d3 * mom(i_part,1) / mom(i_part,0) )
505  k2 = log( 1.d3 * mom(i_part,3) / mom(i_part,2) )
506 
507  mu_phi = - 0.25d0 * ( 5*k2 - k1 ) / log(2.d0)
508 
509  sigma_phi = sqrt( 0.5d0 * (k2-k1) ) / log(2.d0)
510 
511  ELSEIF ( distribution_variable .EQ. 'mass_fraction' ) THEN
512 
513  mu_phi = mom(i_part,1)/mom(i_part,0)
514 
515  sigma_phi = sqrt( -(mom(i_part,1)/mom(i_part,0))**2 + &
516  (mom(i_part,2)/mom(i_part,0)) )
517 
518  skew_phi = ( 2.d0 * (mom0(i_part,1)/mom0(i_part,0))**3 &
519  - 3.d0 * mu_phi * ( mom0(i_part,2) / mom0(i_part,0) ) &
520  + mom0(i_part,3) / mom0(i_part,0) ) / sigma_phi**3
521 
522  END IF
523 
524  description = 'Final Avg Diam '//trim(x1)
525 
526  CALL write_dakota(description,mu_phi)
527 
528  description = 'Final Var Diam '//trim(x1)
529 
530  CALL write_dakota(description,sigma_phi)
531 
532  description = 'Final Skw Diam '//trim(x1)
533 
534  CALL write_dakota(description,skew_phi)
535 
536  mass_fract = solid_partial_mass_fraction(i_part) * ( 1.d0 - &
537  gas_mass_fraction)
538 
539  description = 'Final Mass Fract '//trim(x1)
540 
541  CALL write_dakota(description,mass_fract)
542 
543  solid_mass_flux = solid_partial_mass_fraction(i_part) * ( 1.d0 - &
544  gas_mass_fraction) * rho_mix * pi_g * r**2 * mag_u
545 
546  description = 'Final Mass Flux '//trim(x1)
547 
548  CALL write_dakota(description,solid_mass_flux)
549 
550  solid_mass_flux_change = 1.d0 - solid_mass_flux / solid_mass_flux0
551 
552  description = 'Solid Flux Lost '//trim(x1)
553 
554  CALL write_dakota(description,solid_mass_flux_change)
555 
556 
557  END DO
558 
559  description = 'Objective_function'
560 
561  obj_function = mu_weight * ( ( mu_phi-mu_obj ) / &
562  max( abs(mu_phi),abs(mu_obj) ) )**2 &
563  + sigma_weight * ( ( sigma_phi-sigma_obj ) / &
564  max( abs(sigma_phi),abs(sigma_obj) ) )**2 &
565  + skew_weight * ( ( skew_phi-skew_obj ) / &
566  max( abs(skew_phi),abs(skew_obj) ) )**2 &
567  + height_weight * ( ( plume_height-height_obj ) / &
568  max( abs(plume_height),abs(height_obj) ) )**2
569 
570  CALL write_dakota(description,obj_function)
571 
572  WRITE(*,*) 'plume_height,obj_function' , plume_height , obj_function
573  WRITE(*,*) 'height_nbl',height_nbl
574 
575 
576  RETURN
577 
578  END SUBROUTINE plumerise
579 
580 END MODULE rise
581 !----------------------------------------------------------------------
solver_module::rate
subroutine rate(rhs_)
Compute the right-hand side of the equations.
Definition: solver_rise.f90:87
meteo_module
Meteo module.
Definition: meteo.f90:11
inpout::write_hysplit
subroutine write_hysplit(last)
Write outputs.
Definition: inpout.f90:1639
meteo_module::zmet
subroutine zmet
Meteo parameters.
Definition: meteo.f90:131
variables
Global variables.
Definition: variables.f90:10
particles_module::initialize_particles
subroutine initialize_particles
Particles variables inizialization.
Definition: particles.f90:154
inpout::write_column
subroutine write_column
Write outputs.
Definition: inpout.f90:1548
mixture_module::initialize_mixture
subroutine initialize_mixture
Mixture properties initialization.
Definition: mixture.f90:85
inpout::write_dakota
subroutine write_dakota(description, value)
Dakota outputs.
Definition: inpout.f90:1740
mixture_module
Gas/particles mixture module.
Definition: mixture.f90:11
rise::plumerise
subroutine plumerise
Main subroutine for the integration.
Definition: rise.f90:28
plume_module::initialize_plume
subroutine initialize_plume
Plume variables initialization.
Definition: plume.f90:51
solver_module
Solver module.
Definition: solver_rise.f90:11
plume_module
Plume module.
Definition: plume.f90:11
solver_module::unlump
subroutine unlump(f_)
Calculate physical variables from lumped variables.
Definition: solver_rise.f90:409
solver_module::lump
subroutine lump(f_)
Calculate the lumped variables.
Definition: solver_rise.f90:301
inpout
Input/Output module.
Definition: inpout.f90:11
meteo_module::initialize_meteo
subroutine initialize_meteo
Meteo parameters initialization.
Definition: meteo.f90:109
rise
Predictor-corrector module.
Definition: rise.f90:10
particles_module
Particles module.
Definition: particles.f90:11
solver_module::marching
subroutine marching(fold, fnew, rate)
Marching s one step.
Definition: solver_rise.f90:361
Generated on Wed Jun 1 2016 10:07:08 for PLUME-MoM by   doxygen 1.8.7