PLUME-MoM  1.0
Integralvolcanicplumemodel
 All Classes Files Functions Variables Pages
inpout.f90
Go to the documentation of this file.
1 !********************************************************************
3 !
9 !********************************************************************
10 
11 MODULE inpout
12 
13  USE variables
14 
15  USE plume_module, ONLY: vent_height, alpha_inp , beta_inp , particles_loss ,&
16  r0 , w0 , z , mfr_exp0
17 
18  USE particles_module, ONLY: n_part , n_mom , mom0 , rhop_mom
19 
20  USE particles_module, ONLY : solid_partial_mass_fraction , diam1 , rho1 , &
21  diam2 , rho2 , cp_part , settling_model , distribution , &
22  distribution_variable , solid_mass_fraction , shape_factor
23 
24  USE meteo_module, ONLY: gt , gs , p0 , t0 , h1 , h2 , u_atm0 , duatm_dz0 , &
25  visc_atm0 , rair , cpair , read_atm_profile , u_r , z_r , exp_wind , &
26  wind_mult_coeff
27 
28  USE solver_module, ONLY: ds0
29 
30  USE mixture_module, ONLY: tp0 , gas_mass_fraction0 , rwvapour , cpwvapour , &
31  initial_neutral_density
32 
33  IMPLICIT NONE
34 
36  INTEGER :: n_unit
37 
39  CHARACTER(LEN=30) :: inp_file
40 
42  CHARACTER(LEN=30) :: bak_file
43 
45  CHARACTER(LEN=30) :: run_name
46 
48  CHARACTER(LEN=30) :: col_file
49 
51  CHARACTER(LEN=30) :: hy_file
52 
54  CHARACTER(LEN=30) :: mom_file
55 
57  CHARACTER(LEN=30) :: dak_file
58 
60  CHARACTER(LEN=30) :: mat_file
61 
63  CHARACTER(LEN=50) :: atm_file
64 
66  INTEGER :: atm_unit
67 
68 
70  INTEGER :: bak_unit
71 
73  INTEGER :: mat_unit
74 
76  INTEGER :: inp_unit
77 
79  INTEGER :: col_unit
80 
82  INTEGER :: hy_unit
83 
85  INTEGER :: mom_unit
86 
88  INTEGER :: dak_unit
89 
90  REAL*8, ALLOCATABLE :: mu_lognormal(:) , sigma_lognormal(:)
91 
92  REAL*8 :: month
93  REAL*8 :: lat
94 
95  REAL*8 :: hy_deltaz , hy_z , hy_z_old , hy_x , hy_y , hy_x_old , hy_y_old
96 
97  REAL*8, ALLOCATABLE :: solid_mfr(:) , solid_mfr_old(:)
98 
99  namelist / control_parameters / run_name , verbose_level , dakota_flag , &
100  inversion_flag , hysplit_flag
101 
102  namelist / inversion_parameters / height_weight , height_obj , mu_weight , &
103  mu_obj , sigma_weight , sigma_obj , skew_weight , skew_obj
104 
105  namelist / plume_parameters / alpha_inp , beta_inp , particles_loss
106 
107  namelist / atm_parameters / visc_atm0 , rair , cpair , wind_mult_coeff , &
108  read_atm_profile , settling_model , shape_factor
109 
110  namelist / std_atm_parameters / gt , gs , p0 , t0 , h1 , h2 , u_atm0 , &
111  duatm_dz0 , u_r , z_r , exp_wind
112 
113  namelist / table_atm_parameters / month , lat , u_r , z_r , exp_wind
114 
115  namelist / initial_values / r0 , w0 , mfr_exp0 , tp0 , &
116  initial_neutral_density , gas_mass_fraction0 , vent_height , ds0 , &
117  n_part , distribution , distribution_variable , n_mom
118 
119  namelist / mixture_parameters / diam1 , rho1 , diam2 , rho2 , cp_part , &
120  rwvapour , cpwvapour
121 
122  namelist / lognormal_parameters / solid_partial_mass_fraction , &
123  mu_lognormal , sigma_lognormal
124 
125 
126  SAVE
127 
128 CONTAINS
129 
130 
131  !*****************************************************************************
133  !
140  !******************************************************************************
141 
142  SUBROUTINE initialize
143 
144  ! External procedures
145 
146  USE particles_module, ONLY: allocate_particles
147 
148  IMPLICIT NONE
149 
150  LOGICAL :: lexist
151 
152  gi = 9.81d0 ! Gravity acceleration
153  pi_g = 4.d0 * atan(1.d0)
154 
155  wind_mult_coeff = 1.d0
156 
157 
158  n_unit = 10
159 
160  inp_file = 'plume_model.inp'
161 
162  INQUIRE (file=inp_file,exist=lexist)
163 
164  IF (lexist .EQV. .false.) THEN
165 
166  !
167  !*** Initialization of variables readed in the input file (any version of the
168  !*** input file)
169  !
170 
171  !---------- parameters of the CONTROL_PARAMETERS namelist -------------------
172  run_name = 'default_run'
173  verbose_level = 0
174  dakota_flag = .false.
175  hysplit_flag = .false.
176  inversion_flag = .false.
177 
178  !---------- parameters of the INERSION_PARAMETERS namelist ------------------
179  height_weight = 0.d0
180  height_obj = 0.d0
181  mu_weight = 0.d0
182  mu_obj = 0.d0
183  sigma_weight = 0.d0
184  sigma_obj = 0.d0
185  skew_weight = 0.d0
186  skew_obj = 0.d0
187 
188  !---------- parameters of the PLUME_PARAMETERS namelist ---------------------
189  alpha_inp = 9.0d-2
190  beta_inp = 0.6d0
191  particles_loss = .true.
192 
193  !---------- parameters of the ATM_PARAMETERS namelist -----------------------
194  visc_atm0 = 1.8d-5
195  rair= 287.026
196  cpair= 998.000000
197  wind_mult_coeff = 1.d0
198  read_atm_profile = "standard"
199  settling_model = "textor"
200  shape_factor = 0.43
201 
202  !---------- parameters of the STD_ATM_PARAMETERS namelist -------------------
203  gt = -6.5d-3 ! Temp gradient Troposphere
204  gs = 1.0d-3 ! Temp gradient Stratosphere
205  p0 = 101325.d0 ! Pressure at sea level
206  t0 = 288.15d0 ! Temperature at sea level
207  h1 = 11.d3
208  h2 = 20.d3
209  u_atm0 = 0.d0
210  duatm_dz0 = 0.d0
211 
212  !---------- parameters of the INITIAL_VALUES namelist -----------------------
213 
214  r0 = 0.d0
215  w0 = 0.d0
216  mfr_exp0 = -1.0
217  tp0 = 1273.d0
218  initial_neutral_density = .false.
219  gas_mass_fraction0 = 3.0d-2
220  vent_height = 1500.d0
221  ds0 = 5.d0
222  n_part = 1
223  distribution = 'lognormal'
224  distribution_variable = 'particles_number'
225  n_mom = 6
226 
227  CALL allocate_particles
228 
229  !---------- parameters of the MIXTURE_PARAMETERS namelist -------------------
230 
231  diam1 = 8.d-6
232  rho1 = 2000.d0
233  diam2 = 2.d-3
234  rho2 = 2600.d0
235  cp_part = 1100.d0
236  rwvapour = 462.d0
237  cpwvapour = 1810.d0
238 
239  !---------- parameters of the LOGNORMAL_PARAMETERS namelist -----------------
240 
241  ALLOCATE( mu_lognormal(n_part) )
242  ALLOCATE( sigma_lognormal(n_part) )
243 
244  solid_partial_mass_fraction = 1.d0
245  mu_lognormal= 2.d0
246  sigma_lognormal= 1.6d0
247 
248  inp_unit = n_unit
249 
250  OPEN(inp_unit,file=inp_file,status='NEW')
251 
252  WRITE(inp_unit, control_parameters )
253  WRITE(inp_unit, plume_parameters )
254  WRITE(inp_unit, atm_parameters )
255  WRITE(inp_unit, std_atm_parameters )
256  WRITE(inp_unit, initial_values )
257  WRITE(inp_unit, mixture_parameters )
258  WRITE(inp_unit, lognormal_parameters )
259 
260  CLOSE(inp_unit)
261 
262  WRITE(*,*) 'Input file plume_model.inp not found'
263  WRITE(*,*) 'A new one with default values has been created'
264  stop
265 
266  END IF
267 
268  END SUBROUTINE initialize
269 
270  !******************************************************************************
272  !
278  !******************************************************************************
279 
280  SUBROUTINE read_inp
281 
282  ! External variables
283 
284  USE meteo_module, ONLY: rho_atm , ta , pa , atm_profile , n_atm_profile
285 
286  USE moments_module, ONLY : n_nodes
287 
288  USE solver_module, ONLY: ds0
289 
290  USE mixture_module, ONLY: gas_volume_fraction0 , rgasmix
291 
292  ! External procedures
293 
294  USE meteo_module, ONLY : zmet
295 
296  USE meteo_module, ONLY : h_levels
297 
298  USE meteo_module, ONLY : rho_atm_month_lat , pres_atm_month_lat , temp_atm_month_lat
299 
300  USE moments_module, ONLY : beta_function , wheeler_algorithm , coefficient
301 
302  USE particles_module, ONLY: particles_density , allocate_particles
303 
304  IMPLICIT NONE
305 
306  LOGICAL :: tend1
307  CHARACTER(LEN=80) :: card
308 
309  INTEGER :: i , k , j
310 
311  REAL*8, DIMENSION(max_n_part) :: solid_volume_fraction0
312  REAL*8, ALLOCATABLE :: d_max(:)
313  REAL*8, ALLOCATABLE :: p_beta(:) , q_beta(:)
314 
315  REAL*8, ALLOCATABLE :: mu_bar(:) , sigma_bar(:)
316  REAL*8, ALLOCATABLE :: diam_constant(:)
317  REAL*8, ALLOCATABLE :: diam_constant_phi(:)
318 
319  REAL*8 :: solid_tot_volume_fraction0
320 
321  REAL*8 :: c0
322 
323  REAL*8, DIMENSION(max_n_part) :: rho_solid_avg
324 
325  REAL*8 :: rho_solid_tot_avg
326 
327  REAL*8 :: rho_gas
328  REAL*8 :: rho_mix
329 
330  REAL*8 :: alfa_s
331 
332  REAL*8, ALLOCATABLE :: xi(:) , wi(:)
333  REAL*8, ALLOCATABLE :: part_dens_array(:)
334 
335  REAL*8, ALLOCATABLE :: atm_profile0(:,:)
336 
337  INTEGER :: i_part
338 
339  INTEGER*8 :: fact2
340 
341  INTEGER, ALLOCATABLE :: coeff(:,:)
342 
343  REAL*8, ALLOCATABLE :: rho_atm_month(:,:)
344 
345  REAL*8 :: rho_atm_jan(100,13)
346  REAL*8 :: rho_atm_apr(100,13)
347  REAL*8 :: rho_atm_jul(100,13)
348  REAL*8 :: rho_atm_oct(100,13)
349 
350  REAL*8, ALLOCATABLE :: pres_atm_month(:,:)
351 
352  REAL*8 :: pres_atm_jan(100,13)
353  REAL*8 :: pres_atm_apr(100,13)
354  REAL*8 :: pres_atm_jul(100,13)
355  REAL*8 :: pres_atm_oct(100,13)
356 
357  REAL*8, ALLOCATABLE :: temp_atm_month(:,:)
358 
359  REAL*8 :: temp_atm_jan(100,13)
360  REAL*8 :: temp_atm_apr(100,13)
361  REAL*8 :: temp_atm_jul(100,13)
362  REAL*8 :: temp_atm_oct(100,13)
363 
364  INTEGER :: atm_level
365 
366  INTEGER :: n_atm_levels
367 
368  REAL*8 :: coeff_lat
369 
370  INTEGER :: io
371 
372  namelist / beta_parameters / solid_partial_mass_fraction , p_beta , q_beta ,&
373  d_max
374 
375  namelist / constant_parameters / solid_partial_mass_fraction , &
376  diam_constant_phi
377 
378  namelist / hysplit_parameters / hy_deltaz
379 
380  WRITE(*,*) 'PLUME_MODEL: *** Starting the run ***'
381 
382  n_unit = n_unit + 1
383 
384  inp_unit = n_unit
385 
386  inp_file = 'plume_model.inp'
387 
388  OPEN(inp_unit,file=inp_file,status='old')
389 
390  READ(inp_unit, control_parameters)
391 
392  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read control_parameters: done'
393 
394  IF ( inversion_flag ) READ(inp_unit, inversion_parameters)
395 
396  READ(inp_unit, plume_parameters)
397 
398  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read plume_parameters: done'
399 
400  READ(inp_unit, atm_parameters)
401 
402  IF ( read_atm_profile .EQ. 'table' ) THEN
403 
404  read( inp_unit, table_atm_parameters )
405 
406  n_unit = n_unit + 1
407 
408  atm_unit = n_unit
409 
410  atm_file = '../AtmProfile_info/Density_April.txt'
411 
412  OPEN(atm_unit,file=atm_file,status='old')
413 
414  READ(atm_unit,*)
415 
416  atm_level = 0
417 
418  atm_read_levels_apr: DO
419 
420  atm_level = atm_level + 1
421 
422  ! READ(atm_unit,*, END = 220 ) rho_atm_apr(atm_level,1:8)
423 
424  READ(atm_unit,*,iostat=io ) rho_atm_apr(atm_level,1:8)
425 
426  IF ( io > 0 ) EXIT
427 
428  END DO atm_read_levels_apr
429 
430 220 CLOSE(atm_unit)
431 
432  n_unit = n_unit + 1
433 
434  atm_unit = n_unit
435 
436  atm_file = '../AtmProfile_info/Density_Jan.txt'
437 
438  OPEN(atm_unit,file=atm_file,status='old')
439 
440  READ(atm_unit,*)
441 
442  atm_level = 0
443 
444  atm_read_levels_jan: DO
445 
446  atm_level = atm_level + 1
447 
448  READ(atm_unit,*,iostat=io) rho_atm_jan(atm_level,1:8)
449 
450  IF ( io > 0 ) EXIT
451 
452  END DO atm_read_levels_jan
453 
454 221 CLOSE(atm_unit)
455 
456  n_unit = n_unit + 1
457 
458  atm_unit = n_unit
459 
460  atm_file = '../AtmProfile_info/Density_July.txt'
461 
462  OPEN(atm_unit,file=atm_file,status='old')
463 
464  READ(atm_unit,*)
465 
466  atm_level = 0
467 
468  atm_read_levels_jul: DO
469 
470  atm_level = atm_level + 1
471 
472  READ(atm_unit,*,iostat=io) rho_atm_jul(atm_level,1:8)
473 
474  IF ( io > 0 ) EXIT
475 
476  END DO atm_read_levels_jul
477 
478 222 CLOSE(atm_unit)
479 
480  n_unit = n_unit + 1
481 
482  atm_unit = n_unit
483 
484  atm_file = '../AtmProfile_info/Density_Oct.txt'
485 
486  OPEN(atm_unit,file=atm_file,status='old')
487 
488  READ(atm_unit,*)
489 
490  atm_level = 0
491 
492  atm_read_levels_oct: DO
493 
494  atm_level = atm_level + 1
495 
496  READ(atm_unit,*,iostat=io) rho_atm_oct(atm_level,1:8)
497 
498  IF ( io > 0 ) EXIT
499 
500  n_atm_levels = atm_level
501 
502  END DO atm_read_levels_oct
503 
504 223 CLOSE(atm_unit)
505 
506  ! ----- READ PRESSURES -------
507 
508  n_unit = n_unit + 1
509 
510  atm_unit = n_unit
511 
512  atm_file = '../AtmProfile_info/Pressure_April.txt'
513 
514  OPEN(atm_unit,file=atm_file,status='old')
515 
516  READ(atm_unit,*)
517 
518  atm_level = 0
519 
520  pres_read_levels_apr: DO
521 
522  atm_level = atm_level + 1
523 
524  READ(atm_unit,*,iostat=io) pres_atm_apr(atm_level,1:8)
525 
526  IF ( io > 0 ) EXIT
527 
528  END DO pres_read_levels_apr
529 
530 230 CLOSE(atm_unit)
531 
532  n_unit = n_unit + 1
533 
534  atm_unit = n_unit
535 
536  atm_file = '../AtmProfile_info/Pressure_Jan.txt'
537 
538  OPEN(atm_unit,file=atm_file,status='old')
539 
540  READ(atm_unit,*)
541 
542  atm_level = 0
543 
544  pres_read_levels_jan: DO
545 
546  atm_level = atm_level + 1
547 
548  READ(atm_unit,*,iostat=io) pres_atm_jan(atm_level,1:8)
549 
550  IF ( io > 0 ) EXIT
551 
552  END DO pres_read_levels_jan
553 
554 231 CLOSE(atm_unit)
555 
556  n_unit = n_unit + 1
557 
558  atm_unit = n_unit
559 
560  atm_file = '../AtmProfile_info/Pressure_July.txt'
561 
562  OPEN(atm_unit,file=atm_file,status='old')
563 
564  READ(atm_unit,*)
565 
566  atm_level = 0
567 
568  pres_read_levels_jul: DO
569 
570  atm_level = atm_level + 1
571 
572  READ(atm_unit,*,iostat=io) pres_atm_jul(atm_level,1:8)
573 
574  IF ( io > 0 ) EXIT
575 
576  END DO pres_read_levels_jul
577 
578 232 CLOSE(atm_unit)
579 
580 
581  n_unit = n_unit + 1
582 
583  atm_unit = n_unit
584 
585  atm_file = '../AtmProfile_info/Pressure_Oct.txt'
586 
587  OPEN(atm_unit,file=atm_file,status='old')
588 
589  READ(atm_unit,*)
590 
591  atm_level = 0
592 
593  pres_read_levels_oct: DO
594 
595  atm_level = atm_level + 1
596 
597  READ(atm_unit,*,iostat=io) pres_atm_oct(atm_level,1:8)
598 
599  IF ( io > 0 ) EXIT
600 
601  n_atm_levels = atm_level
602 
603  END DO pres_read_levels_oct
604 
605 233 CLOSE(atm_unit)
606 
607 
608 
609  ! ----- READ TEMPERATURES -------
610 
611  n_unit = n_unit + 1
612 
613  atm_unit = n_unit
614 
615  atm_file = '../AtmProfile_info/Temp_April.txt'
616 
617  OPEN(atm_unit,file=atm_file,status='old')
618 
619  READ(atm_unit,*)
620 
621  atm_level = 0
622 
623  temp_read_levels_apr: DO
624 
625  atm_level = atm_level + 1
626 
627  READ(atm_unit,*,iostat=io) temp_atm_apr(atm_level,1:8)
628 
629  IF ( io > 0 ) EXIT
630 
631  END DO temp_read_levels_apr
632 
633 240 CLOSE(atm_unit)
634 
635  n_unit = n_unit + 1
636 
637  atm_unit = n_unit
638 
639  atm_file = '../AtmProfile_info/Temp_Jan.txt'
640 
641  OPEN(atm_unit,file=atm_file,status='old')
642 
643  READ(atm_unit,*)
644 
645  atm_level = 0
646 
647  temp_read_levels_jan: DO
648 
649  atm_level = atm_level + 1
650 
651  READ(atm_unit,*,iostat=io) temp_atm_jan(atm_level,1:8)
652 
653  IF ( io > 0 ) EXIT
654 
655  END DO temp_read_levels_jan
656 
657 241 CLOSE(atm_unit)
658 
659  n_unit = n_unit + 1
660 
661  atm_unit = n_unit
662 
663  atm_file = '../AtmProfile_info/Temp_July.txt'
664 
665  OPEN(atm_unit,file=atm_file,status='old')
666 
667  READ(atm_unit,*)
668 
669  atm_level = 0
670 
671  temp_read_levels_jul: DO
672 
673  atm_level = atm_level + 1
674 
675  READ(atm_unit,*,iostat=io) temp_atm_jul(atm_level,1:8)
676 
677  IF ( io > 0 ) EXIT
678 
679  END DO temp_read_levels_jul
680 
681 242 CLOSE(atm_unit)
682 
683 
684  n_unit = n_unit + 1
685 
686  atm_unit = n_unit
687 
688  atm_file = '../AtmProfile_info/Temp_Oct.txt'
689 
690  OPEN(atm_unit,file=atm_file,status='old')
691 
692  READ(atm_unit,*)
693 
694  atm_level = 0
695 
696  temp_read_levels_oct: DO
697 
698  atm_level = atm_level + 1
699 
700  READ(atm_unit,*,iostat=io) temp_atm_oct(atm_level,1:8)
701 
702  IF ( io > 0 ) EXIT
703 
704  n_atm_levels = atm_level
705 
706  END DO temp_read_levels_oct
707 
708 243 CLOSE(atm_unit)
709 
710  ALLOCATE( h_levels(n_atm_levels) )
711 
712  ALLOCATE( rho_atm_month_lat(n_atm_levels) , rho_atm_month(n_atm_levels,8) )
713  ALLOCATE( pres_atm_month_lat(n_atm_levels) , pres_atm_month(n_atm_levels,8) )
714  ALLOCATE( temp_atm_month_lat(n_atm_levels) , temp_atm_month(n_atm_levels,8) )
715 
716  IF ((month .GE. 0.d0) .and. (month .LE. 1.d0)) THEN
717  WRITE(*,*) 'winter'
718  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_jan(1:n_atm_levels,1:8)
719  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_jan(1:n_atm_levels,1:8)
720  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_jan(1:n_atm_levels,1:8)
721 
722  ELSEIF ((month .GT. 1.d0) .and. (month .LE. 2.d0)) THEN
723  WRITE(*,*) 'spring'
724  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_apr(1:n_atm_levels,1:8)
725  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_apr(1:n_atm_levels,1:8)
726  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_apr(1:n_atm_levels,1:8)
727 
728  ELSEIF ((month .GT. 2.d0) .and. (month .LE. 3.d0)) THEN
729  WRITE(*,*) 'summer'
730  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_jul(1:n_atm_levels,1:8)
731  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_jul(1:n_atm_levels,1:8)
732  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_jul(1:n_atm_levels,1:8)
733 
734  ELSEIF ((month .GT. 3.d0) .and. (month .LE. 4.d0)) THEN
735  WRITE(*,*) 'autumn'
736  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_apr(1:n_atm_levels,1:8)
737  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_apr(1:n_atm_levels,1:8)
738  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_apr(1:n_atm_levels,1:8)
739 
740  END IF
741 
742  IF ( ( lat .GE. 0.d0 ) .AND. ( lat .LE. 15.d0 ) ) THEN
743 
744  coeff_lat = 1.d0 - ( lat - 0.d0 ) / ( 15.d0 - 0.d0 )
745 
746  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * rho_atm_month(1:n_atm_levels,2) &
747  + ( 1.d0 - coeff_lat ) * rho_atm_month(1:n_atm_levels,3)
748 
749  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * pres_atm_month(1:n_atm_levels,2) &
750  + ( 1.d0 - coeff_lat ) * pres_atm_month(1:n_atm_levels,3)
751 
752  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * temp_atm_month(1:n_atm_levels,2) &
753  + ( 1.d0 - coeff_lat ) * temp_atm_month(1:n_atm_levels,3)
754 
755  ELSEIF ( ( lat .GT. 15.d0 ) .AND. ( lat .LE. 30.d0 ) ) THEN
756 
757  coeff_lat = 1.d0 - ( lat - 15.d0 ) / ( 30.d0 - 15.d0 )
758 
759  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * rho_atm_month(1:n_atm_levels,3) &
760  + ( 1.d0 - coeff_lat ) * rho_atm_month(1:n_atm_levels,4)
761 
762  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * pres_atm_month(1:n_atm_levels,3) &
763  + ( 1.d0 - coeff_lat ) * pres_atm_month(1:n_atm_levels,5)
764 
765  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * temp_atm_month(1:n_atm_levels,3) &
766  + ( 1.d0 - coeff_lat ) * temp_atm_month(1:n_atm_levels,5)
767 
768  ELSEIF ( ( lat .GT. 30.d0 ) .AND. ( lat .LE. 45.d0 ) ) THEN
769 
770  coeff_lat = 1.d0 - ( lat - 30.d0 ) / ( 45.d0 - 30.d0 )
771 
772  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * rho_atm_month(1:n_atm_levels,4) &
773  + ( 1.d0 - coeff_lat ) * rho_atm_month(1:n_atm_levels,5)
774 
775  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * pres_atm_month(1:n_atm_levels,4) &
776  + ( 1.d0 - coeff_lat ) * pres_atm_month(1:n_atm_levels,5)
777 
778  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * temp_atm_month(1:n_atm_levels,4) &
779  + ( 1.d0 - coeff_lat ) * temp_atm_month(1:n_atm_levels,5)
780 
781  ELSEIF ( ( lat .GT. 45.d0 ) .AND. ( lat .LE. 60.d0 ) ) THEN
782 
783  coeff_lat = 1.d0 - ( lat - 45.d0 ) / ( 60.d0 - 45.d0 )
784 
785  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * rho_atm_month(1:n_atm_levels,5) &
786  + ( 1.d0 - coeff_lat ) * rho_atm_month(1:n_atm_levels,6)
787 
788  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * pres_atm_month(1:n_atm_levels,5) &
789  + ( 1.d0 - coeff_lat ) * pres_atm_month(1:n_atm_levels,6)
790 
791  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * temp_atm_month(1:n_atm_levels,5) &
792  + ( 1.d0 - coeff_lat ) * temp_atm_month(1:n_atm_levels,6)
793 
794  ELSEIF ( ( lat .GT. 60.d0 ) .AND. ( lat .LE. 75.d0 ) ) THEN
795 
796  coeff_lat = 1.d0 - ( lat - 60.d0 ) / ( 75.d0 - 60.d0 )
797 
798  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * rho_atm_month(1:n_atm_levels,6) &
799  + ( 1.d0 - coeff_lat ) * rho_atm_month(1:n_atm_levels,7)
800 
801  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * pres_atm_month(1:n_atm_levels,6) &
802  + ( 1.d0 - coeff_lat ) * pres_atm_month(1:n_atm_levels,7)
803 
804  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * temp_atm_month(1:n_atm_levels,6) &
805  + ( 1.d0 - coeff_lat ) * temp_atm_month(1:n_atm_levels,7)
806 
807  ELSEIF ( ( lat .GT. 75.d0 ) .AND. ( lat .LE. 90.d0 ) ) THEN
808 
809  coeff_lat = 1.d0 - ( lat - 75.d0 ) / ( 90.d0 - 75.d0 )
810 
811  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * rho_atm_month(1:n_atm_levels,7) &
812  + ( 1.d0 - coeff_lat ) * rho_atm_month(1:n_atm_levels,8)
813 
814  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * pres_atm_month(1:n_atm_levels,7) &
815  + ( 1.d0 - coeff_lat ) * pres_atm_month(1:n_atm_levels,8)
816 
817  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * temp_atm_month(1:n_atm_levels,7) &
818  + ( 1.d0 - coeff_lat ) * temp_atm_month(1:n_atm_levels,8)
819 
820  END IF
821 
822  pres_atm_month_lat(1:n_atm_levels) = 100.d0 * pres_atm_month_lat(1:n_atm_levels)
823 
824  h_levels(1:n_atm_levels) = 1000.d0 * temp_atm_month(1:n_atm_levels,1)
825 
826  !WRITE(*,*) 'rho_atm_month_lat(1:n_atm_levels)', rho_atm_month_lat(1:n_atm_levels)
827 
828  ELSEIF ( read_atm_profile .EQ. 'card' ) THEN
829 
830  tend1 = .false.
831 
832  WRITE(*,*) 'search atm_profile'
833 
834  atm_profile_search: DO
835 
836  READ(inp_unit,*, END = 200 ) card
837 
838  IF( trim(card) == 'ATM_PROFILE' ) THEN
839 
840  EXIT atm_profile_search
841 
842  END IF
843 
844  END DO atm_profile_search
845 
846  READ(inp_unit,*) n_atm_profile
847 
848  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'n_atm_profile',n_atm_profile
849 
850  ALLOCATE( atm_profile(7,n_atm_profile) )
851  ALLOCATE( atm_profile0(7,n_atm_profile) )
852 
853  DO i = 1, n_atm_profile
854 
855  READ(inp_unit,*) atm_profile0(1:7,i)
856 
857 
858  atm_profile(1:7,i) = atm_profile0(1:7,i)
859  ! convert from km to meters
860  atm_profile(1,i) = atm_profile(1,i) * 1000.d0
861 
862  ! convert from hPa to Pa
863  atm_profile(3,i) = atm_profile(3,i) * 100.d0
864 
865  atm_profile(6,i) = atm_profile(6,i) * wind_mult_coeff
866  atm_profile(7,i) = atm_profile(7,i) * wind_mult_coeff
867 
868  IF ( verbose_level .GE. 1 ) WRITE(*,*) i,atm_profile(1,i)
869 
870  END DO
871 
872  goto 210
873 200 tend1 = .true.
874 210 CONTINUE
875 
876  rewind(inp_unit)
877 
878  ELSEIF ( read_atm_profile .EQ. 'standard' ) THEN
879 
880  read( inp_unit,std_atm_parameters )
881 
882  END IF
883 
884  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read atm_parameters: done'
885 
886  READ(inp_unit, initial_values)
887 
888  IF ( ( mfr_exp0 .LT. 0.d0 ) .AND. ( r0 .EQ. 0.d0 ) .AND. ( w0 .GT. 0.d0 ) ) THEN
889 
890  WRITE(*,*) 'WARNING: initial radius calculated from MER and velocity'
891 
892  END IF
893 
894  IF ( ( mfr_exp0 .LT. 0.d0 ) .AND. ( r0 .EQ. 0.d0 ) .AND. ( w0 .EQ. 0.d0 ) ) THEN
895 
896  WRITE(*,*) 'WARNING: initial radius and velocity calculated from MER and gas mass fraction'
897  stop
898 
899  END IF
900 
901  IF ( ( mfr_exp0 .GT. 0.d0 ) .AND. ( w0 .GT. 0.d0 ) .AND. ( r0 .GT. 0.d0 ) ) THEN
902 
903  WRITE(*,*) 'ERROR: too many unknown input parameters: input mfr_exp0 or w0 and r0'
904  stop
905 
906  END IF
907 
908  IF ( distribution .EQ. 'constant' ) THEN
909 
910  n_mom = 4
911  n_nodes = 1
912 
913  ELSE
914 
915  n_nodes = 0.5d0 * n_mom
916 
917  END IF
918 
919  IF ( hysplit_flag ) THEN
920 
921  IF ( distribution .NE. 'constant' ) THEN
922 
923  WRITE(*,*) 'hysplit run requires constant distribution'
924  stop
925 
926  ELSE
927 
928  READ(inp_unit, hysplit_parameters)
929 
930  ALLOCATE( solid_mfr(n_part) , solid_mfr_old(n_part) )
931 
932  hy_z = vent_height + hy_deltaz
933  hy_z_old = vent_height
934  hy_x_old = 0.d0
935  hy_y_old = 0.d0
936 
937  END IF
938 
939  END IF
940 
941  z = vent_height
942 
943  CALL allocate_particles
944 
945  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read initial_parameters: done'
946 
947  READ(inp_unit, mixture_parameters)
948 
949  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read mixture_parameters: done'
950 
951 
952  IF ( distribution .EQ. 'beta' ) THEN
953 
954  ALLOCATE( p_beta(n_part) )
955  ALLOCATE( q_beta(n_part) )
956  ALLOCATE( d_max(n_part) )
957 
958  READ(inp_unit, beta_parameters)
959 
960  ELSEIF ( distribution .EQ. 'lognormal' ) THEN
961 
962  ALLOCATE( mu_lognormal(n_part) )
963  ALLOCATE( sigma_lognormal(n_part) )
964 
965  ALLOCATE( mu_bar(n_part) )
966  ALLOCATE( sigma_bar(n_part) )
967 
968  READ(inp_unit, lognormal_parameters)
969 
970  mu_bar = -log( 2.d0 ) * mu_lognormal
971  sigma_bar = log( 2.d0 ) * sigma_lognormal
972 
973  WRITE(*,*) 'mu_bar',mu_bar
974  WRITE(*,*) 'sigma_bar',sigma_bar
975 
976 
977  ELSEIF ( distribution .EQ. 'constant' ) THEN
978 
979  ALLOCATE( diam_constant(n_part) )
980  ALLOCATE( diam_constant_phi(n_part) )
981 
982  READ(inp_unit, constant_parameters)
983 
984  diam_constant = 1.d-3 * 2.d0**(-diam_constant_phi)
985 
986  END IF
987 
988  IF ( sum( solid_partial_mass_fraction(1:n_part) ) .NE. 1.d0 ) THEN
989 
990  WRITE(*,*) 'WARNING: Sum of solid mass fractions > 1'
991 
992  END IF
993 
994  solid_partial_mass_fraction(n_part) = 1.d0 - &
995  sum( solid_partial_mass_fraction(1:n_part-1) )
996 
997  solid_mass_fraction(1:n_part) = ( 1.d0 - gas_mass_fraction0 ) * &
998  solid_partial_mass_fraction(1:n_part)
999 
1000  ALLOCATE( xi(n_nodes) )
1001  ALLOCATE( wi(n_nodes) )
1002  ALLOCATE( part_dens_array(n_nodes) )
1003 
1004  ! evaluate the moments from the parameters of the beta distribution
1005  ! these moments have to be corrected for the mass fractions give in
1006  ! the input file
1007 
1008  IF ( distribution_variable .EQ. 'mass_fraction' ) THEN
1009 
1010  ALLOCATE( coeff(0:n_mom,0:n_mom) )
1011  CALL coefficient(n_mom,coeff)
1012 
1013  END IF
1014 
1015  DO i_part = 1,n_part
1016 
1017  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'i_part',i_part
1018 
1019  DO i = 0, n_mom-1
1020 
1021  IF ( distribution_variable .EQ. 'particles_number' ) THEN
1022 
1023  IF ( distribution .EQ. 'beta' ) THEN
1024 
1025  mom0(i_part,i) = d_max(i_part)**i *beta_function(p_beta(i_part) &
1026  +i,q_beta(i_part)) / beta_function(p_beta(i_part) , &
1027  q_beta(i_part))
1028 
1029  ELSEIF ( distribution .EQ. 'lognormal' ) THEN
1030 
1031  mom0(i_part,i) = 6.d0 / pi_g * 10.d0**(-3*(i-3)) * &
1032  exp( ( i-3.d0 ) * mu_bar(i_part) + ( i - 3.d0 )**2 / 2.d0 *&
1033  sigma_bar(i_part)**2 )
1034 
1035  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'before correction', &
1036  mom0(i_part,i)
1037 
1038  ELSEIF ( distribution .EQ. 'constant' ) THEN
1039 
1040  mom0(i_part,i) = diam_constant(i_part)**i
1041 
1042  END IF
1043 
1044  ELSEIF ( distribution_variable .EQ. 'mass_fraction' ) THEN
1045 
1046  IF ( distribution .EQ. 'lognormal' ) THEN
1047 
1048  IF ( mu_lognormal(i_part) .EQ. 0.d0) mu_lognormal(i_part) = 1.d-5
1049 
1050  mom0(i_part,i) = 0.d0
1051 
1052  DO k=0,floor(0.5d0 * i)
1053 
1054  fact2 = product((/(j, j = 2*k-1,1,-2)/))
1055 
1056  mom0(i_part,i) = mom0(i_part,i) + coeff(i,2*k) * fact2 &
1057  * sigma_lognormal(i_part)**(2*k) * &
1058  mu_lognormal(i_part)**(i-2*k)
1059 
1060  IF ( verbose_level .GT. 3 ) THEN
1061 
1062  WRITE(*,*) 'i,k,coeff_bin,fact2',i,k,coeff(i,2*k),fact2
1063  WRITE(*,*) 'mom0',mom0(i_part,i)
1064  READ(*,*)
1065 
1066  END IF
1067 
1068  END DO
1069 
1070  ELSE
1071 
1072  END IF
1073 
1074  END IF
1075 
1076  END DO
1077 
1078  CALL zmet
1079 
1080  IF ( initial_neutral_density ) THEN
1081 
1082  rgasmix = rair
1083 
1084  ELSE
1085 
1086  rgasmix = rwvapour
1087 
1088  END IF
1089 
1090  rho_gas = pa / ( rgasmix * tp0 )
1091 
1092  IF ( distribution .EQ. 'constant' ) THEN
1093 
1094  CALL wheeler_algorithm( mom0(i_part,0:1) , xi , wi )
1095 
1096  ELSE
1097 
1098  CALL wheeler_algorithm( mom0(i_part,0:n_mom-1) , xi , wi )
1099 
1100  END IF
1101 
1102 
1103  DO i=1,n_nodes
1104 
1105  part_dens_array(i) = particles_density( i_part , xi(i) )
1106 
1107  END DO
1108 
1109  IF ( verbose_level .GE. 1 ) THEN
1110 
1111  WRITE(*,*) 'part_dens_array'
1112  WRITE(*,*) 'xi',xi
1113  WRITE(*,*) 'wi',wi
1114  WRITE(*,*) 'rho',part_dens_array
1115 
1116  END IF
1117 
1118  ! the density of the particles phases are evaluated here. It is
1119  ! independent from the mass fraction of the particles phases, so
1120  ! it is possible to evaluate them with the "uncorrected" moments
1121 
1122  IF ( distribution_variable .EQ. 'particles_number' ) THEN
1123 
1124  rho_solid_avg(i_part) = sum( part_dens_array * wi * xi**3 ) / &
1125  mom0(i_part,3)
1126 
1127  ELSEIF ( distribution_variable .EQ. 'mass_fraction' ) THEN
1128 
1129  rho_solid_avg(i_part) = 1.d0 / ( sum( wi / part_dens_array ) / &
1130  mom0(i_part,0) )
1131 
1132  ELSE
1133 
1134  WRITE(*,*) 'input_file: distribution_variable',distribution_variable
1135  stop
1136 
1137  END IF
1138 
1139  IF ( verbose_level .GE. 1 ) THEN
1140 
1141  WRITE(*,*) 'rho avg',rho_solid_avg(i_part)
1142 
1143  END IF
1144 
1145  END DO
1146 
1147  ! the average solid density is evaluated through the mass fractions and
1148  ! the densities of the particles phases
1149 
1150  rho_solid_tot_avg = 1.d0 / sum( solid_partial_mass_fraction(1:n_part) / &
1151  rho_solid_avg(1:n_part) )
1152 
1153  IF ( verbose_level .GE. 1 ) THEN
1154 
1155  WRITE(*,*)
1156  WRITE(*,*) '******* CHECK ON MASS AND VOLUME FRACTIONS *******'
1157  WRITE(*,*) 'rho solid avg', rho_solid_tot_avg
1158 
1159  END IF
1160 
1161  IF ( initial_neutral_density ) THEN
1162 
1163  rho_mix = rho_atm
1164 
1165  solid_tot_volume_fraction0 = ( rho_mix - rho_gas ) / &
1166  ( rho_solid_tot_avg - rho_gas )
1167 
1168  gas_volume_fraction0 = 1.d0 - solid_tot_volume_fraction0
1169 
1170  gas_mass_fraction0 = gas_volume_fraction0 * rho_gas / rho_mix
1171 
1172  ELSE
1173 
1174  gas_volume_fraction0 = rho_solid_tot_avg / ( rho_gas * ( 1.d0 / &
1175  gas_mass_fraction0 - 1.d0 ) + rho_solid_tot_avg )
1176 
1177  solid_tot_volume_fraction0 = 1.d0 - gas_volume_fraction0
1178 
1179  rho_mix = gas_volume_fraction0 * rho_gas + solid_tot_volume_fraction0 &
1180  * rho_solid_tot_avg
1181 
1182  END IF
1183 
1184  IF ( verbose_level .GE. 1 ) THEN
1185 
1186  WRITE(*,*) 'gas_volume_fraction0',gas_volume_fraction0
1187  WRITE(*,*) 'solid_tot_volume_fraction0',solid_tot_volume_fraction0
1188  WRITE(*,*) 'rho_gas',rho_gas
1189  WRITE(*,*) 'rho_mix',rho_mix
1190  WRITE(*,*)
1191 
1192  END IF
1193 
1194  DO i_part = 1,n_part
1195 
1196  ! the volume fraction of the particle phases ( with respect to the
1197  ! solid phase only) is evaluated
1198 
1199  alfa_s = solid_partial_mass_fraction(i_part) * rho_solid_tot_avg / &
1200  rho_solid_avg(i_part)
1201 
1202  ! this is the volume fraction of the particles phases in the mixture
1203 
1204  solid_volume_fraction0(i_part) = solid_tot_volume_fraction0 * alfa_s
1205 
1206  ! the coefficient C0 (=mom0) for the particles size distribution is
1207  ! evaluated in order to have the corrected volume or mass fractions
1208 
1209  IF ( distribution_variable .EQ. 'particles_number' ) THEN
1210 
1211  c0 = 6.d0 / pi_g * solid_volume_fraction0(i_part) / mom0(i_part,3)
1212 
1213  ELSEIF ( distribution_variable .EQ. 'mass_fraction' ) THEN
1214 
1215  c0 = ( 1.d0 - gas_mass_fraction0 ) / mom0(i_part,0) * &
1216  solid_partial_mass_fraction(i_part)
1217 
1218  END IF
1219 
1220  ! the moments are corrected with the factor C0
1221 
1222  DO i = 0, n_mom-1
1223 
1224  mom0(i_part,i) = c0 * mom0(i_part,i)
1225 
1226  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'mom',i,mom0(i_part,i)
1227 
1228  END DO
1229 
1230  IF ( verbose_level .GE. 1 ) THEN
1231 
1232  WRITE(*,*) 'i_part =',i_part
1233  WRITE(*,*) 'alfa_s',i_part,alfa_s
1234  WRITE(*,*) 'solid_volume_fraction0',solid_volume_fraction0(i_part)
1235  WRITE(*,*) 'solid_partial_mass_fract', &
1236  solid_partial_mass_fraction(i_part)
1237  WRITE(*,*) 'solid_mass_fract', solid_mass_fraction(i_part)
1238  WRITE(*,*)
1239 
1240  END IF
1241 
1242  END DO
1243 
1244  gas_volume_fraction0 = 1.d0 - solid_tot_volume_fraction0
1245 
1246  gas_mass_fraction0 = gas_volume_fraction0 * rho_gas / rho_mix
1247 
1248  WRITE(*,*) 'solid volume fraction',solid_tot_volume_fraction0
1249  WRITE(*,*) 'solid total mass_fraction', solid_tot_volume_fraction0 * &
1250  rho_solid_tot_avg / rho_mix
1251 
1252  IF ( distribution_variable .EQ. 'particles_number' ) THEN
1253 
1254  WRITE(*,*) 'solid_mass_fractions', mom0(1:n_part,3) * &
1255  rho_solid_avg(1:n_part) / ( sum( mom0(1:n_part,3) * &
1256  rho_solid_avg(1:n_part)) )
1257 
1258  ELSEIF ( distribution_variable .EQ. 'mass_fraction' ) THEN
1259 
1260  WRITE(*,*) 'solid_mass_fractions', mom0(1:n_part,0)
1261 
1262  END IF
1263 
1264  WRITE(*,*) 'gas volume fraction', gas_volume_fraction0
1265  WRITE(*,*) 'gas mass fraction', gas_mass_fraction0
1266 
1267  ! the parameters of the particles phases distributions are saved in a file
1268  ! readable by Matlab
1269 
1270  IF ( .NOT.dakota_flag ) THEN
1271 
1272  n_unit = n_unit + 1
1273 
1274  mat_unit = n_unit
1275 
1276  mat_file = trim(run_name)//'.m'
1277 
1278  OPEN(mat_unit,file=mat_file,status='unknown')
1279 
1280  WRITE(mat_unit,*) 'n_part = ',n_part,';'
1281  WRITE(mat_unit,*) 'gas_volume_fraction = ',gas_volume_fraction0,';'
1282 
1283  IF ( distribution .EQ. 'beta' ) THEN
1284 
1285  WRITE(mat_unit,*) 'p = [',p_beta(1:n_part),'];'
1286  WRITE(mat_unit,*) 'q = [',q_beta(1:n_part),'];'
1287  WRITE(mat_unit,*) 'd_max = [',d_max(1:n_part),'];'
1288 
1289 
1290  ELSEIF ( distribution .EQ. 'lognormal' ) THEN
1291 
1292  WRITE(mat_unit,*) 'mu = [',mu_lognormal(1:n_part),'];'
1293  WRITE(mat_unit,*) 'sigma = [',sigma_lognormal(1:n_part),'];'
1294 
1295  ELSEIF ( distribution .EQ. 'constant' ) THEN
1296 
1297  WRITE(mat_unit,*) 'diam = [',diam_constant(1:n_part),'];'
1298 
1299  END IF
1300 
1301  WRITE(mat_unit,*) 'solid_mass_fractions = [', &
1302  solid_partial_mass_fraction(1:n_part),'];'
1303 
1304  WRITE(mat_unit,*) 'd1 = [',diam1(1:n_part),'];'
1305  WRITE(mat_unit,*) 'd2 = [',diam2(1:n_part),'];'
1306  WRITE(mat_unit,*) 'rho1 = [',rho1(1:n_part),'];'
1307  WRITE(mat_unit,*) 'rho2 = [',rho2(1:n_part),'];'
1308 
1309 
1310  IF ( verbose_level .GE. 0 ) WRITE(*,*) 'write matlab file: done'
1311 
1312  CLOSE(mat_unit)
1313 
1314  END IF
1315 
1316  IF ( distribution .EQ. 'moments' ) THEN
1317 
1318  tend1 = .false.
1319 
1320  moments_search: DO
1321 
1322  READ(inp_unit , *, END = 300 ) card
1323 
1324  IF( trim(card) == 'MOMENTS' ) THEN
1325 
1326  EXIT moments_search
1327 
1328  END IF
1329 
1330  END DO moments_search
1331 
1332  READ(inp_unit,*) n_mom
1333 
1334  WRITE(*,*) 'input_moments'
1335 
1336  READ(inp_unit,*) solid_partial_mass_fraction(1:n_part)
1337 
1338  DO i = 0, n_mom-1
1339 
1340  READ(inp_unit,*) mom0(1:n_part,i)
1341 
1342  WRITE(*,*) mom0(1:n_part,i)
1343 
1344  END DO
1345 
1346  goto 310
1347 300 tend1 = .true.
1348 310 CONTINUE
1349 
1350  END IF
1351 
1352  ! Close input file
1353 
1354  CLOSE(inp_unit)
1355 
1356  ! Write a backup of the input file
1357 
1358  n_unit = n_unit + 1
1359 
1360  bak_unit = n_unit
1361 
1362  bak_file = trim(run_name)//'.bak'
1363 
1364  OPEN(bak_unit,file=bak_file,status='unknown')
1365 
1366  WRITE(bak_unit, control_parameters)
1367 
1368  IF ( inversion_flag ) WRITE(bak_unit, inversion_parameters)
1369 
1370  WRITE(bak_unit, plume_parameters)
1371 
1372  WRITE(bak_unit, atm_parameters)
1373 
1374  IF ( read_atm_profile .EQ. 'standard' ) THEN
1375 
1376  WRITE(bak_unit, std_atm_parameters )
1377 
1378  END IF
1379 
1380  IF ( read_atm_profile .EQ. 'table' ) THEN
1381 
1382  WRITE(bak_unit, table_atm_parameters )
1383 
1384  END IF
1385 
1386 
1387  WRITE(bak_unit, initial_values)
1388 
1389  WRITE(bak_unit, mixture_parameters)
1390 
1391  IF ( distribution .EQ. 'beta' ) THEN
1392 
1393  WRITE(bak_unit, beta_parameters)
1394 
1395  ELSEIF ( distribution .EQ. 'lognormal' ) THEN
1396 
1397  WRITE(bak_unit, lognormal_parameters)
1398 
1399  ELSEIF ( distribution .EQ. 'constant' ) THEN
1400 
1401  WRITE(bak_unit, constant_parameters)
1402 
1403  ELSEIF ( distribution .EQ. 'moments' ) THEN
1404 
1405  IF (( tend1 ) .OR. ( n_mom .EQ. 0 )) THEN
1406 
1407  WRITE(*,*) 'WARNING: input ', ' SAMPLING POINTS not found '
1408 
1409  ELSE
1410 
1411  WRITE(bak_unit,*) '''MOMENTS'''
1412  WRITE(bak_unit,*) n_mom
1413 
1414  DO i = 0, n_mom-1
1415 
1416  WRITE(bak_unit,*) mom0(1:n_part,i)
1417 
1418  END DO
1419 
1420  END IF
1421 
1422  END IF
1423 
1424  IF ( read_atm_profile .EQ. 'card' ) THEN
1425 
1426  WRITE(bak_unit,*) '''ATM_PROFILE'''
1427  WRITE(bak_unit,*) n_atm_profile
1428 
1429  DO i = 1, n_atm_profile
1430 
1431  WRITE(bak_unit,107) atm_profile0(1:7,i)
1432 
1433 107 FORMAT(7(1x,e14.7))
1434 
1435 
1436  END DO
1437 
1438  END IF
1439 
1440  CLOSE(bak_unit)
1441 
1442  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'end subroutine reainp'
1443 
1444  RETURN
1445 
1446  END SUBROUTINE read_inp
1447 
1448  !*****************************************************************************
1450  !
1455  !******************************************************************************
1456 
1457  SUBROUTINE open_file_units
1458 
1459  ! External variables
1460  USE particles_module, ONLY : n_part
1461  USE moments_module, ONLY : n_mom
1462  USE variables, ONLY : dakota_flag , hysplit_flag
1463 
1464  IMPLICIT NONE
1465 
1466 
1467  col_file = trim(run_name)//'.col'
1468  mom_file = trim(run_name)//'.mom'
1469  dak_file = trim(run_name)//'.dak'
1470  hy_file = trim(run_name)//'.hy'
1471 
1472  IF ( .NOT.dakota_flag ) THEN
1473 
1474  n_unit = n_unit + 1
1475  col_unit = n_unit
1476 
1477  OPEN(col_unit,file=col_file)
1478 
1479  n_unit = n_unit + 1
1480  mom_unit = n_unit
1481 
1482  OPEN(mom_unit,file=mom_file)
1483 
1484 
1485  WRITE(mom_unit,*) n_part
1486  WRITE(mom_unit,*) n_mom
1487 
1488  END IF
1489 
1490  IF ( hysplit_flag ) THEN
1491 
1492  n_unit = n_unit + 1
1493  hy_unit = n_unit
1494 
1495  OPEN(hy_unit,file=hy_file)
1496 
1497  END IF
1498 
1499  n_unit = n_unit + 1
1500  dak_unit = n_unit
1501 
1502  OPEN(dak_unit,file=dak_file)
1503 
1504  RETURN
1505 
1506  END SUBROUTINE open_file_units
1507 
1508  !*****************************************************************************
1510  !
1515  !******************************************************************************
1516 
1517  SUBROUTINE close_file_units
1518 
1519  USE variables, ONLY : dakota_flag , hysplit_flag
1520 
1521  IMPLICIT NONE
1522 
1523  IF ( .not.dakota_flag ) THEN
1524 
1525  CLOSE(col_unit)
1526  CLOSE(mom_unit)
1527 
1528  END IF
1529 
1530  IF ( hysplit_flag ) CLOSE ( hy_unit )
1531 
1532  CLOSE(dak_unit)
1533 
1534  RETURN
1535 
1536  END SUBROUTINE close_file_units
1537 
1538  !*****************************************************************************
1540  !
1546  !******************************************************************************
1547 
1548  SUBROUTINE write_column
1549 
1550  USE meteo_module, ONLY: rho_atm , ta, pa
1551 
1552  USE particles_module, ONLY: n_mom , n_part , solid_partial_mass_fraction , &
1553  mom , set_mom
1554 
1555  USE plume_module, ONLY: x , y , z , w , r , mag_u
1556  USE mixture_module, ONLY: rho_mix , tp , gas_mass_fraction , &
1557  atm_mass_fraction , wvapour_mass_fraction
1558 
1559  USE variables, ONLY: verbose_level
1560 
1561  IMPLICIT NONE
1562 
1563  REAL*8 :: mfr
1564 
1565  INTEGER :: i_part
1566 
1567  mfr = 3.14 * r**2 * rho_mix * mag_u
1568 
1569  ! WRITE(*,*) 'INPOUT: atm_mass_fraction',atm_mass_fraction
1570  ! READ(*,*)
1571 
1572  IF ( z .EQ. vent_height ) THEN
1573 
1574  WRITE(col_unit,97,advance="no")
1575 
1576  DO i_part=1,n_part
1577 
1578  WRITE(col_unit,98,advance="no")
1579 
1580  END DO
1581 
1582  WRITE(col_unit,99)
1583 
1584 97 FORMAT(1x,' z (m) ',1x,' r (m) ',1x,' x (m) ', &
1585  1x,' y (m) ',1x,'mix.dens(kg/m3)',1x,'temperature(C)', &
1586  1x,' vert vel (m/s)',1x,' mag vel (m/s) ',1x,' atm.mass fract', &
1587  1x,' wv mass fract ',1x)
1588 
1589 98 FORMAT(1x,'sol.mass fract ')
1590 
1591 99 FORMAT(1x,'atm.rho(kg/m3)',1x,' MFR (kg/s) ',1x,'atm.temp (K) ', &
1592  1x,' atm pres (Pa) ')
1593 
1594 
1595  END IF
1596 
1597  WRITE(col_unit,100) z , r , x , y , rho_mix , tp - 273.15d0 , w , mag_u,&
1598  atm_mass_fraction , wvapour_mass_fraction , &
1599  solid_partial_mass_fraction(1:n_part) , rho_atm , mfr , ta, pa
1600 
1601 !********* format for plume models comparison ********************************
1602 !
1603 ! WRITE(col_unit,100) z , r , x , y , rho_mix , tp - 273.15D0 , w , &
1604 ! atm_mass_fraction , wvapour_mass_fraction , &
1605 ! solid_partial_mass_fraction(1:n_part) * ( 1.D0 - gas_mass_fraction)&
1606 ! , rho_atm , mfr , ta
1607 
1608  WRITE(mom_unit,*) z , mom(1:n_part,0:n_mom-1),set_mom(1:n_part,0)
1609 
1610 100 FORMAT(33(1x,e15.8))
1611 
1612  IF ( verbose_level .GE. 1 ) THEN
1613 
1614  WRITE(*,*) '******************'
1615  WRITE(*,*) 'z',z
1616  WRITE(*,*) 'x',x
1617  WRITE(*,*) 'y',y
1618  WRITE(*,*) 'r',r
1619  WRITE(*,*) 'w',w
1620  WRITE(*,*) '******************'
1621  READ(*,*)
1622 
1623  END IF
1624 
1625  RETURN
1626 
1627  END SUBROUTINE write_column
1628 
1629  !*****************************************************************************
1631  !
1637  !******************************************************************************
1638 
1639  SUBROUTINE write_hysplit(last)
1640 
1641  USE meteo_module, ONLY: rho_atm , ta, pa
1642 
1643  USE particles_module, ONLY: n_mom , n_part , solid_partial_mass_fraction , &
1644  mom
1645 
1646  USE plume_module, ONLY: x , y , z , w , r , mag_u
1647  USE mixture_module, ONLY: rho_mix , tp , gas_mass_fraction , &
1648  atm_mass_fraction , wvapour_mass_fraction
1649 
1650 
1651  IMPLICIT NONE
1652 
1653  LOGICAL, INTENT(IN) :: last
1654 
1655  INTEGER :: i_part
1656 
1657  CHARACTER(len=8) :: x1 ! format descriptor
1658 
1659 
1660  IF ( z .EQ. vent_height ) THEN
1661 
1662  solid_mfr(1:n_part) = solid_partial_mass_fraction(1:n_part) * ( 1.d0 - &
1663  gas_mass_fraction) * pi_g * mag_u * r**2.d0 * rho_mix
1664 
1665  WRITE(hy_unit,107,advance="no")
1666 
1667  DO i_part=1,n_part
1668 
1669  WRITE(x1,'(I2.2)') i_part ! converting integer to string using a 'internal file'
1670 
1671  WRITE(hy_unit,108,advance="no") 'S mfr'//trim(x1)//' (kg/s)'
1672 
1673  END DO
1674 
1675  WRITE(hy_unit,*) ''
1676 
1677  ELSEIF ( z .GE. hy_z ) THEN
1678 
1679  solid_mfr_old(1:n_part) = solid_mfr(1:n_part)
1680 
1681  solid_mfr(1:n_part) = solid_partial_mass_fraction(1:n_part) * ( 1.d0 - &
1682  gas_mass_fraction) * pi_g * mag_u * r**2.d0 * rho_mix
1683 
1684  WRITE(hy_unit,110) 0.5d0 * ( hy_x + hy_x_old ) , 0.5d0 * ( hy_y + hy_y_old ) , &
1685  0.5d0 * ( hy_z + hy_z_old ) , solid_mfr_old(1:n_part) &
1686  - solid_mfr(1:n_part)
1687 
1688 ! WRITE(*,*) solid_partial_mass_fraction(1:n_part)
1689 ! WRITE(*,*) z, ( 1.D0 - gas_mass_fraction) * rho_mix * ( pi_g * r**2.D0 ) * mag_u
1690 
1691  hy_x_old = hy_x
1692  hy_y_old = hy_y
1693  hy_z_old = hy_z
1694 
1695  hy_x = x
1696  hy_y = y
1697  hy_z = hy_z + hy_deltaz
1698 
1699  END IF
1700 
1701  IF ( last ) THEN
1702 
1703  solid_mfr_old(1:n_part) = solid_mfr(1:n_part)
1704 
1705  solid_mfr(1:n_part) = solid_partial_mass_fraction(1:n_part) * ( 1.d0 - &
1706  gas_mass_fraction) * pi_g * mag_u * r**2.d0 * rho_mix
1707 
1708  hy_x = x
1709  hy_y = y
1710  hy_z = z
1711 
1712  WRITE(hy_unit,110) 0.5d0 * ( hy_x + hy_x_old ) , 0.5d0 * ( hy_y + &
1713  hy_y_old ) , 0.5d0 * ( hy_z + hy_z_old ) , solid_mfr_old(1:n_part) &
1714  - solid_mfr(1:n_part)
1715 
1716  END IF
1717 
1718 107 FORMAT(1x,' x (m) ',1x,' y (m) ', 1x,' z (m) ')
1719 
1720 108 FORMAT(2x,a)
1721 
1722 110 FORMAT(33(1x,e15.8))
1723 
1724  RETURN
1725 
1726  END SUBROUTINE write_hysplit
1727 
1728  !*****************************************************************************
1730  !
1738  !******************************************************************************
1739 
1740  SUBROUTINE write_dakota(description,value)
1741 
1742  USE variables, ONLY : verbose_level
1743 
1744  IMPLICIT NONE
1745 
1746  CHARACTER(20), INTENT(IN) :: description
1747  REAL*8, INTENT(IN) :: value
1748 
1749  WRITE(dak_unit,*) description,value
1750 
1751  IF ( verbose_level .GE. 2 ) THEN
1752 
1753  WRITE(*,*) description,value
1754 
1755  END IF
1756 
1757  RETURN
1758 
1759  END SUBROUTINE write_dakota
1760 
1761 
1762 
1763 END MODULE inpout
particles_module::particles_density
real *8 function particles_density(i_part, diam_in)
Particle density.
Definition: particles.f90:440
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
inpout::write_column
subroutine write_column
Write outputs.
Definition: inpout.f90:1548
inpout::initialize
subroutine initialize
Initialize variables.
Definition: inpout.f90:142
inpout::write_dakota
subroutine write_dakota(description, value)
Dakota outputs.
Definition: inpout.f90:1740
moments_module::coefficient
subroutine coefficient(nmax, c)
Binomial Coefficient.
Definition: moments.f90:414
mixture_module
Gas/particles mixture module.
Definition: mixture.f90:11
particles_module::allocate_particles
subroutine allocate_particles
Particles variables inizialization.
Definition: particles.f90:113
solver_module
Solver module.
Definition: solver_rise.f90:11
inpout::close_file_units
subroutine close_file_units
Close output units.
Definition: inpout.f90:1517
plume_module
Plume module.
Definition: plume.f90:11
inpout::open_file_units
subroutine open_file_units
Initialize output units.
Definition: inpout.f90:1457
moments_module
Method of Moments module.
Definition: moments.f90:11
inpout
Input/Output module.
Definition: inpout.f90:11
moments_module::beta_function
real *8 function beta_function(z, w)
Beta function.
Definition: moments.f90:316
particles_module
Particles module.
Definition: particles.f90:11
inpout::read_inp
subroutine read_inp
Read Input data.
Definition: inpout.f90:280
moments_module::wheeler_algorithm
subroutine wheeler_algorithm(m, xi, w)
Wheeler algorithm.
Definition: moments.f90:41
Generated on Wed Jun 1 2016 10:07:08 for PLUME-MoM by   doxygen 1.8.7