PLUME-MoM-TSM  1.0
VolcanicPlumeModel
inpout.f90
Go to the documentation of this file.
1 !********************************************************************
3 !
9 !********************************************************************
10 
11 MODULE inpout
12 
13  use, intrinsic :: iso_fortran_env
14  use, intrinsic :: ieee_arithmetic
15 
16  USE variables
17 
18  USE parameters_2d, ONLY : t_start , t_end , dt_output , wp , n_vars , c_d
19 
20  ! -- Variables for the namelist NUMERIC_PARAMETERS
21  USE parameters_2d, ONLY : rsource_cells , solver_scheme, dt0 , max_dt , cfl,&
22  limiter , theta, reconstr_coeff , interfaces_relaxation , n_rk
23 
24 
25  USE moments_module, ONLY : n_mom , n_nodes , n_sections
26 
27  USE plume_module, ONLY: vent_height, alpha_inp , beta_inp , particles_loss ,&
28  r0 , w0 , z , log10_mfr
29 
30  USE particles_module, ONLY: n_part , mom0 , mom
31 
32  USE particles_module, ONLY : solid_partial_mass_fraction , phi1 , rho1 , &
33  phi2 , rho2 , cp_part , settling_model , distribution , &
34  solid_mass_fraction0 , shape_factor , bin_partial_mass_fraction , &
35  solid_partial_mass_fraction0 , shape1 , shape2
36 
37  USE particles_module, ONLY : aggregation_model , particles_beta0 , &
38  phil , phir , m
39 
40  USE meteo_module, ONLY: h1 , h2 , rh , sphu_atm0 , visc_atm0 , &
41  rair , cpair , read_atm_profile , u_max , z_r , exp_wind , &
42  wind_mult_coeff , rwv , rel_hu
43 
44  USE solver_module, ONLY: dz0
45 
46  USE mixture_module, ONLY: t_mix0 , water_mass_fraction0, &
47  initial_neutral_density
48 
49  USE mixture_module, ONLY: n_gas , rvolcgas , cpvolcgas , rvolcgas_mix , &
50  volcgas_mass_fraction , volcgas_mix_mass_fraction , cpvolcgas_mix , &
51  rhovolcgas_mix , volcgas_mol_wt , rhovolcgas , volcgas_mass_fraction0, &
52  rho_lw, rho_ice, added_water_temp, added_water_mass_fraction
53 
54  IMPLICIT NONE
55 
56  REAL(wp) :: notset
57 
58 
60  INTEGER :: n_unit
61 
63  CHARACTER(LEN=30) :: inp_file
64 
66  CHARACTER(LEN=30) :: bak_file
67 
69  CHARACTER(LEN=30) :: run_name
70 
72  CHARACTER(LEN=30) :: sed_file
73 
75  CHARACTER(LEN=30) :: col_file
76 
78  CHARACTER(LEN=30) :: hy_file
79 
81  CHARACTER(LEN=30) :: hy_file_volcgas
82 
84  CHARACTER(LEN=30) :: mom_file
85 
87  CHARACTER(LEN=30) :: dak_file
88 
90  CHARACTER(LEN=30) :: inversion_file
91 
93  CHARACTER(LEN=50) :: atm_file
94 
96  INTEGER :: atm_unit
97 
98 
100  INTEGER :: bak_unit
101 
103  INTEGER :: inp_unit
104 
106  INTEGER :: col_unit
107 
109  INTEGER :: sed_unit
110 
112  INTEGER :: hy_unit
113 
114  INTEGER :: hy_lines
115 
116  INTEGER :: col_lines
117 
119 
120  INTEGER :: hy_unit_volcgas
121 
123  INTEGER :: temp_unit
124 
126  INTEGER :: mom_unit
127 
129  INTEGER :: dak_unit
130 
132  INTEGER :: inversion_unit
133 
134  REAL(wp) :: mfr0
135 
136  REAL(wp), ALLOCATABLE :: mu_lognormal(:) , sigma_lognormal(:)
137 
138  REAL(wp) :: month
139  REAL(wp) :: lat
140 
141  REAL(wp) :: phi_min , delta_phi
142 
143  REAL(wp) :: hy_deltaz , hy_z , hy_z_old , hy_x , hy_y , hy_x_old , hy_y_old
144 
145  REAL(wp), ALLOCATABLE :: solid_mfr(:) , solid_mfr_old(:), solid_mfr_init(:) , &
146  solid_mfr_oldold(:)
147 
148  REAL(wp) :: p_atm0 , t_atm0
149 
150  namelist / control_parameters / run_name , verbose_level , dakota_flag , &
151  inversion_flag , hysplit_flag , aggregation_flag, water_flag , &
152  umbrella_flag , entr_abv_nbl_flag
153 
154  namelist / mom_parameters / n_part , n_mom , n_nodes , n_sections
155 
156  namelist / particles_parameters / phi_min , delta_phi , distribution , &
157  solid_partial_mass_fraction , phi1 , rho1 , phi2 , rho2 , shape1 , &
158  shape2 , cp_part , shape_factor , particles_loss , settling_model
159 
160  namelist / inversion_parameters / height_obj , r_min , r_max , n_values , &
161  w_min , w_max , nbl_stop
162 
163  namelist / entrainment_parameters / alpha_inp , beta_inp
164 
165  namelist / water_parameters / rho_lw , rho_ice , added_water_temp , &
166 
167  added_water_mass_fraction
168 
169  namelist / atm_parameters / visc_atm0 , rair , cpair , wind_mult_coeff , &
170  read_atm_profile
171 
172  namelist / std_atm_parameters / sphu_atm0 , u_max , p_atm0 , t_atm0 , rel_hu
173 
174  namelist / table_atm_parameters / month , lat , u_max , z_r , exp_wind
175 
176  namelist / initial_values / r0 , w0 , log10_mfr , mfr0 , t_mix0 , &
177  initial_neutral_density , water_mass_fraction0 , vent_height , dz0 , &
178  n_gas
179 
180  namelist / aggregation_parameters / aggregation_model , particles_beta0
181 
182  namelist / hysplit_parameters / hy_deltaz , nbl_stop , n_cloud
183 
184  namelist / volcgas_parameters / rvolcgas , cpvolcgas , volcgas_mol_wt , &
185  volcgas_mass_fraction0
186 
187  namelist / lognormal_parameters / mu_lognormal , sigma_lognormal
188 
189  namelist / umbrella_run_parameters / t_end , dt_output , c_d , steady_flag
190 
191  namelist / numeric_parameters / rsource_cells , solver_scheme, dt0 , max_dt , &
192  cfl, limiter , theta , reconstr_coeff , interfaces_relaxation , n_rk
193 
194  SAVE
195 
196 CONTAINS
197 
198 
199  !*****************************************************************************
201  !
208  !******************************************************************************
209 
210  SUBROUTINE initialize
211 
212  ! External procedures
213 
214  USE particles_module, ONLY: allocate_particles
215 
216 
217  IMPLICIT NONE
218 
219  LOGICAL :: lexist
220 
221  REAL(wp) :: test
222 
223  notset = ieee_value(0.0_wp, ieee_quiet_nan)
224 
225 
226  !---------- default flags of the CONTROL_PARAMETERS namelist ----------------
227  dakota_flag = .false.
228  hysplit_flag = .false.
229  inversion_flag = .false.
230  aggregation_flag = .false.
231  water_flag = .false.
232  umbrella_flag = .false.
233  entr_abv_nbl_flag = .false.
234 
235  !------- default parameters of the ENTRAINMENT_PARAMETERS namelist ----------
236  alpha_inp = notset
237  beta_inp = notset
238 
239  !-------------- default of the INVERSION_PARAMETERS namelist ----------------
240  height_obj = notset
241  r_min = notset
242  r_max = notset
243  w_min = notset
244  w_max = notset
245  n_values = 0
246 
247  !-------------- default values of the INITIAL_VALUES namelist ---------------
248  initial_neutral_density = .false.
249  r0 = notset
250  w0 = notset
251  log10_mfr = notset
252  mfr0 = notset
253 
254  !-------------- default values of the AGGEGATION_PARAMETERS namelist --------
255  particles_beta0 = notset
256 
257  !-------------- default values of the STD_ATM_PARAMETERS namelist --------
258  sphu_atm0 = notset
259  rel_hu = notset
260  p_atm0 = notset
261  t_atm0 = notset
262  u_max = notset
263 
264  !------------ default values of the HYSPLIT_PARAMETERS namelist -------------
265  hy_deltaz = notset
266  nbl_stop = .true.
267  n_cloud = -1
268 
269  !---------- default values of the WATER_PARAMETERS namelist -----------------
270  rho_lw = notset
271  rho_ice = notset
272  added_water_temp = notset
273  added_water_mass_fraction = notset
274 
275 
276  !---------- default values of the UMBRELLA namelist -------------------------
277  t_end = 3600
278  dt_output = 600
279  c_d = notset
280  steady_flag = .false.
281 
282  !---------- default values of the ATM_PARAMETERS namelist -------------------
283  visc_atm0 = notset
284  rair = notset
285  cpair = notset
286  wind_mult_coeff = notset
287  read_atm_profile = ""
288  settling_model = "none"
289 
290  gi = 9.81_wp ! Gravity acceleration
291  pi_g = 4.0_wp * atan(1.0_wp)
292 
293  n_unit = 10
294 
295  inp_file = 'plume_model.inp'
296 
297  INQUIRE (file=inp_file,exist=lexist)
298 
299  IF (lexist .EQV. .false.) THEN
300 
301  !
302  !*** Initialization of variables readed in the input file (any version of
303  !*** the input file)
304  !
305 
306  !---------- parameters of the CONTROL_PARAMETERS namelist ----------------
307  run_name="default_run"
308  verbose_level = 0
309  dakota_flag = .false.
310  inversion_flag = .false.
311  hysplit_flag = .false.
312  water_flag = .false.
313  aggregation_flag = .false.
314 
315  !---------- parameters of the MOM_PARAMETERS namelist --------------------
316  n_part = 1
317  n_mom = 2
318  n_nodes = 5
319  n_sections = 11
320 
321  CALL allocate_particles
322 
323  !---------- parameters of the PARTICLES_PARAMETERS namelist --------------
324  phi_min = -4.0_wp
325  delta_phi = 1.0_wp
326  distribution = 'LOGNORMAL'
327  solid_partial_mass_fraction = 1.0_wp
328  phi1 = -1.0_wp
329  rho1 = 2000.0_wp
330  phi2 = 4.0_wp
331  rho2 = 2600.0_wp
332  shape1 = 1.0_wp
333  shape2 = 1.0_wp
334  shape_factor = 1.0_wp
335  cp_part = 1100.0_wp
336  particles_loss= .true.
337  settling_model="textor"
338 
339  !---------- parameters of the ENTRAINMENT_PARAMETERS namelist ------------
340  alpha_inp = 9.0e-2_wp
341  beta_inp = 0.60_wp
342 
343  !---------- parameters of the WATER_PARAMETERS namelist ------------------
344  rho_lw = 1000.0_wp
345  rho_ice = 920.0_wp
346  added_water_temp = 273.0_wp
347  added_water_mass_fraction = 0.0_wp
348 
349  !---------- parameters of the ATM_PARAMETERS namelist --------------------
350  visc_atm0 = 1.8e-5_wp
351  rair= 287.026_wp
352  cpair= 998.0_wp
353  wind_mult_coeff = 1.0_wp
354  read_atm_profile = "standard"
355  settling_model = "textor"
356 
357  !---------- parameters of the STD_ATM_PARAMETERS namelist ----------------
358  sphu_atm0 = 0.0_wp
359  u_max = 5.0_wp
360  z_r = 1000.0_wp
361  exp_wind = 0.0_wp
362 
363  !---------- parameters of the INITIAL_VALUES namelist --------------------
364  r0= 50.0_wp
365  w0 = 130.0_wp
366  t_mix0= 1373.0_wp
367  initial_neutral_density = .false.
368  water_mass_fraction0= 3.0e-2_wp
369  vent_height= 1500.0_wp
370  dz0= 5.0_wp
371  n_gas= 2
372 
373  ALLOCATE ( rvolcgas(n_gas) , cpvolcgas(n_gas) , volcgas_mol_wt(n_gas) , &
374  volcgas_mass_fraction(n_gas) , volcgas_mass_fraction0(n_gas) , &
375  rhovolcgas(n_gas) )
376 
377  !---------- parameters of the VOLCGAS_PARAMETERS namelist ----------------
378  rvolcgas(1) = 189.0_wp
379  rvolcgas(2) = 130.0_wp
380  cpvolcgas(1)= 844.0_wp
381  cpvolcgas(2)= 640.0_wp
382  volcgas_mol_wt(1) = 0.044_wp
383  volcgas_mol_wt(2) = 0.064_wp
384  volcgas_mass_fraction0(1) = 0.005_wp
385  volcgas_mass_fraction0(2) = 0.005_wp
386 
387  !---------- parameters of the LOGNORMAL_PARAMETERS namelist --------------
388  ALLOCATE( mu_lognormal(n_part) )
389  ALLOCATE( sigma_lognormal(n_part) )
390  mu_lognormal= 2.0_wp
391  sigma_lognormal= 1.6_wp
392 
393  inp_unit = n_unit
394 
395  OPEN(inp_unit,file=inp_file,status='NEW')
396 
397  WRITE(inp_unit, control_parameters )
398  WRITE(inp_unit, mom_parameters )
399  WRITE(inp_unit, particles_parameters )
400  WRITE(inp_unit, entrainment_parameters )
401  WRITE(inp_unit, atm_parameters )
402  WRITE(inp_unit, std_atm_parameters )
403  WRITE(inp_unit, initial_values )
404  WRITE(inp_unit, volcgas_parameters )
405  WRITE(inp_unit, lognormal_parameters )
406 
407  CLOSE(inp_unit)
408 
409  WRITE(*,*) 'Input file plume_model.inp not found'
410  WRITE(*,*) 'A new one with default values has been created'
411  stop
412 
413  END IF
414 
415  END SUBROUTINE initialize
416 
417  !******************************************************************************
419  !
425  !******************************************************************************
426 
427  SUBROUTINE read_inp
428 
429  ! External variables
430 
431  USE meteo_module, ONLY: rho_atm , pa , ta , atm_profile , n_atm_profile
432 
433  USE mixture_module, ONLY: water_volume_fraction0 , rgasmix , &
434  gas_mass_fraction, water_vapor_mass_fraction , &
435  liquid_water_mass_fraction , gas_volume_fraction, ice_mass_fraction
436 
437  USE particles_module, ONLY : m_quad , w_quad , f_quad , rho_quad
438 
439  ! External procedures
440 
441  USE meteo_module, ONLY : initialize_meteo
442 
443  USE meteo_module, ONLY : h_levels
444 
445  USE meteo_module, ONLY : rho_atm_month_lat , pres_atm_month_lat , &
446  temp_atm_month_lat
447 
448  ! USE mixture_module, ONLY : eval_wv
449 
450  USE particles_module, ONLY: particles_density , allocate_particles , &
451  deallocate_particles , eval_quad_values , init_quadrature_points , &
452  phifromm, particles_shape
453 
454  IMPLICIT NONE
455 
456  LOGICAL :: tend1
457  CHARACTER(LEN=80) :: card
458 
459  INTEGER :: ios
460 
461  INTEGER :: i , k , j
462 
463  REAL(wp), DIMENSION(max_n_part) :: solid_volume_fraction0
464  REAL(wp), ALLOCATABLE :: d_max(:)
465 
466  REAL(wp) :: solid_tot_volume_fraction0
467 
468  REAL(wp), DIMENSION(max_n_part) :: rho_solid_avg
469 
470  REAL(wp) :: rho_solid_tot_avg
471 
472  REAL(wp) :: diam
473 
474  REAL(wp) :: rhowv
475  REAL(wp) :: rho_gas
476  REAL(wp) :: rho_mix
477 
478  REAL(wp) :: alfa_s
479 
480 
481  REAL(wp), ALLOCATABLE :: atm_profile0(:,:)
482 
483  INTEGER :: i_part
484 
485  INTEGER, ALLOCATABLE :: coeff(:,:)
486 
487  REAL(wp), ALLOCATABLE :: rho_atm_month(:,:)
488 
489  REAL(wp) :: rho_atm_jan(100,13)
490  REAL(wp) :: rho_atm_apr(100,13)
491  REAL(wp) :: rho_atm_jul(100,13)
492  REAL(wp) :: rho_atm_oct(100,13)
493 
494  REAL(wp), ALLOCATABLE :: pres_atm_month(:,:)
495 
496  REAL(wp) :: pres_atm_jan(100,13)
497  REAL(wp) :: pres_atm_apr(100,13)
498  REAL(wp) :: pres_atm_jul(100,13)
499  REAL(wp) :: pres_atm_oct(100,13)
500 
501  REAL(wp), ALLOCATABLE :: temp_atm_month(:,:)
502 
503  REAL(wp) :: temp_atm_jan(100,13)
504  REAL(wp) :: temp_atm_apr(100,13)
505  REAL(wp) :: temp_atm_jul(100,13)
506  REAL(wp) :: temp_atm_oct(100,13)
507 
508  INTEGER :: atm_level
509 
510  INTEGER :: n_atm_levels
511 
512  REAL(wp) :: coeff_lat
513 
514  REAL(wp) :: Rrhovolcgas_mix
515 
516  INTEGER :: io
517 
518  INTEGER :: i_gas
519 
520  INTEGER :: i_sect
521 
522  INTEGER :: ip
523 
524  REAL(wp) :: rhop
525  REAL(wp) :: shapep
526 
527  namelist / bin_parameters / bin_partial_mass_fraction
528 
529  IF ( write_flag ) THEN
530 
531  WRITE(*,*)
532  WRITE(*,*) 'PlumeMoM (by M. de'' Michieli Vitturi)'
533  WRITE(*,*)
534  WRITE(*,*) '*** Starting the run ***'
535  WRITE(*,*)
536 
537  END IF
538 
539  n_unit = n_unit + 1
540 
541  inp_unit = n_unit
542 
543  inp_file = 'plume_model.inp'
544 
545  OPEN(inp_unit,file=inp_file,status='old')
546 
547  READ(inp_unit, control_parameters,iostat=io)
548 
549  IF ( io .EQ. 0 ) THEN
550 
551  n_unit = n_unit + 1
552  bak_unit = n_unit
553  bak_file = trim(run_name)//'.bak'
554 
555  OPEN(bak_unit,file=bak_file,status='unknown')
556  WRITE(bak_unit, control_parameters)
557  rewind(inp_unit)
558 
559  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read control_parameters: done'
560 
561  ELSE
562 
563  WRITE(*,*) 'Problem with namelist CONTROL_PARAMETERS'
564  stop
565 
566  END IF
567 
568  IF ( inversion_flag ) THEN
569 
570  READ(inp_unit, inversion_parameters,iostat=ios)
571 
572  IF ( ios .NE. 0 ) THEN
573 
574  WRITE(*,*) 'IOSTAT=',ios
575  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
576  WRITE(*,*) 'Please check the input file'
577  WRITE(*,inversion_parameters)
578  stop
579 
580  END IF
581 
582  READ(inp_unit, initial_values )
583 
584  IF ( ios .NE. 0 ) THEN
585 
586  WRITE(*,*) 'IOSTAT=',ios
587  WRITE(*,*) 'ERROR: problem with namelist INITIAL_VALUES'
588  WRITE(*,*) 'Please check the input file'
589  WRITE(*,initial_values)
590  stop
591 
592  END IF
593 
594  rewind(inp_unit)
595 
596  IF ( ( .NOT.isset(w0) ) .AND. ( .NOT.isset(r0) ) ) THEN
597 
598  IF ( umbrella_flag ) THEN
599 
600  WRITE(*,*) 'ERROR: problem with INVERSION'
601  WRITE(*,*) 'Search for both radius and velocity is not possible'
602  WRITE(*,*) 'with UMBRELLA_FLAG = T'
603  WRITE(*,*) 'Please check the input file'
604  stop
605 
606  END IF
607 
608  END IF
609 
610  IF ( ( .NOT. isset(height_obj) ) .OR. ( height_obj .LE. 0 ) ) THEN
611 
612  WRITE(*,*) ''
613  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
614  WRITE(*,*)
615  WRITE(*,inversion_parameters)
616  WRITE(*,*)
617  WRITE(*,*) 'Please check HEIGHT_OBJ value (>0 [m])'
618  WRITE(*,*) 'HEIGHT_OBJ =',height_obj
619  WRITE(*,*)
620  stop
621 
622  END IF
623 
624  IF ( verbose_level.GE.1 ) WRITE(*,*) 'read inversion_parameters: done'
625  WRITE(bak_unit, inversion_parameters)
626  write_flag = .false.
627  rewind(inp_unit)
628 
629  ELSE
630 
631  write_flag = .true.
632 
633  END IF
634 
635  !----------------- READ ENTRAINMENT_PARAMETERS namelist ---------------------
636  READ(inp_unit, entrainment_parameters,iostat=io)
637 
638  IF ( io .EQ. 0 ) THEN
639 
640  WRITE(bak_unit, entrainment_parameters)
641  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read entrainment_parameters: done'
642 
643  IF ( .NOT.isset(alpha_inp) ) THEN
644 
645  WRITE(*,*) 'ERROR: problem with namelist ENTRAINMENT_PARAMETERS'
646  WRITE(*,*) 'Please set alpha_inp (>0)'
647  WRITE(*,*)
648  WRITE(*,entrainment_parameters)
649  stop
650 
651  END IF
652 
653  IF ( .NOT.isset(beta_inp) ) THEN
654 
655  WRITE(*,*) 'ERROR: problem with namelist ENTRAINMENT_PARAMETERS'
656  WRITE(*,*) 'Please set beta_inp (>0)'
657  WRITE(*,*)
658  WRITE(*,entrainment_parameters)
659  stop
660 
661  END IF
662 
663  rewind(inp_unit)
664 
665  ELSE
666 
667  WRITE(*,*) 'IOSTAT=',ios
668  WRITE(*,*) 'ERROR: problem with namelist ENTRAINMENT_PARAMETERS'
669  WRITE(*,*) 'Please check the input file'
670  WRITE(*,entrainment_parameters)
671  rewind(inp_unit)
672  stop
673 
674  END IF
675 
676  !------------------- READ MoM_PARAMETERS namelist ---------------------------
677  READ(inp_unit, mom_parameters,iostat=io)
678 
679  IF ( io .EQ. 0 ) THEN
680 
681  WRITE(bak_unit, mom_parameters)
682 
683  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read MoM_parameters: done'
684 
685  CALL allocate_particles
686 
687  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'allocated particles parameters'
688 
689  rewind(inp_unit)
690 
691  ELSE
692 
693  WRITE(*,*) 'Problem with namelist MoM_PARAMETERS'
694  stop
695 
696  END IF
697 
698  !------------------- READ PARTICLES_PARAMETERS namelist ---------------------
699  READ(inp_unit, particles_parameters,iostat=io)
700 
701  IF ( io .EQ. 0 ) THEN
702 
703  WRITE(bak_unit, particles_parameters)
704  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read particles_parameters: done'
705  rewind(inp_unit)
706 
707  ELSE
708 
709  WRITE(*,*) 'Problem with namelist PARTICLES_PARAMETERS'
710  WRITE(*,*)
711  WRITE(*,particles_parameters)
712  stop
713 
714  END IF
715 
716  ! Compute the bins in phi-scale according to input parameters
717  DO i_sect=1,n_sections
718 
719  phil(i_sect) = phi_min + (n_sections-i_sect) * delta_phi
720  phir(i_sect) = phi_min + (n_sections-i_sect+1) * delta_phi
721 
722  END DO
723 
724  IF ( verbose_level .GE. 1 ) THEN
725 
726  WRITE(*,*) 'grain size sections:'
727  WRITE(*,"(100F6.2)") phil
728  WRITE(*,"(100F6.2)") phir
729 
730  END IF
731 
732  ! Compute the mass instervals for the different particles (1,n_part)
733  DO i_part = 1,n_part
734 
735  m(1,i_part) = 0.0_wp
736 
737  ! check on phi1 and phi2
738  IF ( isset(phi1(i_part)) .AND. isset(phi2(i_part)) ) THEN
739 
740  IF ( phi1(i_part) .GT. phi2(i_part) ) THEN
741 
742  WRITE(*,*) 'ERROR: problem with namelist PARTICLES_PARAMETERS'
743  WRITE(*,*) 'Please check the input file'
744  WRITE(*,*) 'phi1 MUST BE SMALLER THAN phi2'
745  WRITE(*,*) 'phi1',phi1
746  WRITE(*,*) 'phi2',phi2
747  stop
748 
749  END IF
750 
751  ELSE
752 
753  WRITE(*,*) 'ERROR: problem with namelist PARTICLES_PARAMETERS'
754  WRITE(*,*) 'Please check the input file'
755  WRITE(*,*) 'phi1',phi1
756  WRITE(*,*) 'phi2',phi2
757  stop
758 
759  END IF
760 
761  IF ( isset(shape_factor(i_part)) ) THEN
762 
763  IF ( isset(shape1(i_part)) .AND. isset(shape2(i_part)) ) THEN
764 
765  WRITE(*,*) 'ERROR: problem with namelist PARTICLES_PARAMETERS'
766  WRITE(*,*) 'Please check the input file'
767  WRITE(*,*) 'Shape factor',shape_factor
768  WRITE(*,*) 'Shape1',shape1
769  WRITE(*,*) 'Shape2',shape2
770  WRITE(*,*) 'IF SHAPE_FACTOR is defined, SHAPE1 and SHAPE2'
771  WRITE(*,*) 'are not used'
772  stop
773 
774  ELSE
775 
776  shape1(i_part) = shape_factor(i_part)
777  shape2(i_part) = shape_factor(i_part)
778 
779  END IF
780 
781  ELSE
782 
783  IF ( isset(shape1(i_part)) .AND. isset(shape2(i_part)) ) THEN
784 
785 
786  ELSE
787 
788  WRITE(*,*) 'ERROR: problem with namelist PARTICLES_PARAMETERS'
789  WRITE(*,*) 'Please check the input file'
790  WRITE(*,*) 'Shape factor',shape_factor
791  WRITE(*,*) 'Shape1',shape1
792  WRITE(*,*) 'Shape2',shape2
793  rewind(inp_unit)
794  stop
795 
796  END IF
797 
798  END IF
799 
800  DO i_sect = 1,n_sections
801 
802  diam = 1.e-3_wp * 2.0_wp**( - phir(i_sect) )
803  rhop = particles_density( i_part,phir(i_sect) )
804  shapep = particles_shape( i_part,phir(i_sect) )
805 
806  m(i_sect+1,i_part) = rhop * (shapep * diam**3)
807 
808  END DO
809 
810  END DO
811 
812  CALL init_quadrature_points
813 
814  distribution = lower(distribution)
815 
816  IF ( distribution .EQ. 'lognormal' ) THEN
817 
818  ALLOCATE( mu_lognormal(n_part) )
819  ALLOCATE( sigma_lognormal(n_part) )
820 
821  READ(inp_unit, lognormal_parameters,iostat=io)
822 
823  IF ( io .EQ. 0 ) THEN
824 
825  WRITE(bak_unit, lognormal_parameters)
826  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read lognormal_parameters: done'
827  rewind(inp_unit)
828 
829  ELSE
830 
831  WRITE(*,*) 'Problem with namelist LOGNORMAL_PARAMETERS'
832  stop
833 
834  END IF
835 
836  ELSEIF ( distribution .EQ. 'bin' ) THEN
837 
838  READ(inp_unit, bin_parameters,iostat=io)
839 
840  IF ( io .EQ. 0 ) THEN
841 
842  WRITE(bak_unit, bin_parameters)
843  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read bin_parameters: done'
844  rewind(inp_unit)
845 
846  DO i_part=1,n_part
847  bin_partial_mass_fraction(1:n_sections,i_part) = &
848  bin_partial_mass_fraction(1:n_sections,i_part) / &
849  sum(bin_partial_mass_fraction(1:n_sections,i_part))
850 
851  END DO
852 
853  ELSE
854 
855  WRITE(*,*) 'Problem reading namelist BIN_PARAMETERS'
856  WRITE(*,*) 'Please check the input file'
857  WRITE(*,*)
858  WRITE(*,bin_parameters)
859  stop
860 
861  END IF
862 
863  ELSE
864 
865  WRITE(*,*) 'Error in namelist PARTICLES_PARAMETERS'
866  WRITE(*,*) 'Please check the values of distribution: ',distribution
867  WRITE(*,particles_parameters)
868  stop
869 
870  END IF
871 
872  !--------------------- READ WATER_PARAMETERS namelist -----------------------
873  IF (water_flag) THEN
874 
875  READ(inp_unit, water_parameters,iostat=io)
876 
877  IF ( io .EQ. 0 ) THEN
878 
879  IF ( .NOT.isset(rho_lw) ) THEN
880 
881  WRITE(*,*) 'Namelist WATER_PARAMETERS'
882  WRITE(*,*) 'Plase define RHO_LW (kg/m3)'
883  stop
884 
885  END IF
886 
887  IF ( .NOT.isset(rho_ice) ) THEN
888 
889  WRITE(*,*) 'Namelist WATER_PARAMETERS'
890  WRITE(*,*) 'Plase define RHO_ICE (kg/m3)'
891  stop
892 
893  END IF
894 
895  IF ( .NOT.isset(added_water_mass_fraction) ) THEN
896 
897  WRITE(*,*) ''
898  WRITE(*,*) 'ERROR: problem with namelist WATER_PARAMETERS'
899  WRITE(*,*)
900  WRITE(*,water_parameters)
901  WRITE(*,*)
902  WRITE(*,*) 'Please check ADDED_WATER_MASS_FRACTION value [0;1]'
903  WRITE(*,*) 'ADDED_WATER_MASS_FRACTION =',added_water_mass_fraction
904  WRITE(*,*)
905  stop
906 
907  ELSE
908 
909  IF ( ( added_water_mass_fraction .LT. 0.0_wp ) .OR. &
910  ( added_water_mass_fraction .GE. 1.0_wp ) ) THEN
911 
912  WRITE(*,*) 'Namelist WATER_PARAMETERS'
913  WRITE(*,*) 'added_water_mass_fraction should be >=0 and <1'
914  WRITE(*,*) 'actual value:',added_water_mass_fraction
915  stop
916 
917  END IF
918 
919  END IF
920 
921  IF ( .NOT.isset(added_water_temp) ) THEN
922 
923  WRITE(*,*) ''
924  WRITE(*,*) 'ERROR: problem with namelist WATER_PARAMETERS'
925  WRITE(*,*)
926  WRITE(*,water_parameters)
927  WRITE(*,*)
928  WRITE(*,*) 'Please check ADDED_WATER_TEMP value [K]'
929  WRITE(*,*) 'ADDED_WATER_TEMP =',added_water_temp
930  WRITE(*,*)
931  stop
932 
933  END IF
934 
935  WRITE(bak_unit, water_parameters)
936 
937  rewind(inp_unit)
938 
939  ELSE
940 
941  WRITE(*,*) 'Problem reading namelist WATER_PARAMETERS'
942  WRITE(*,*) 'Please check the input file'
943  WRITE(*,*)
944  WRITE(*,water_parameters)
945  stop
946 
947  END IF
948 
949  ELSE
950 
951  rho_ice = 920.0_wp
952  rho_lw = 1000.0_wp
953  added_water_mass_fraction = 0.0_wp
954  added_water_temp = 273.0_wp
955 
956  END IF
957 
958  READ(inp_unit, atm_parameters,iostat=io)
959 
960  IF ( io .EQ. 0 ) THEN
961 
962  IF ( .NOT.isset(added_water_temp) ) THEN
963 
964  WRITE(*,*) ''
965  WRITE(*,*) 'ERROR: problem with namelist ATM_PARAMETERS'
966  WRITE(*,*)
967  WRITE(*,atm_parameters)
968  WRITE(*,*)
969  WRITE(*,*) 'Please check VISC_ATM0 value [Pa s]'
970  WRITE(*,*) 'VISC_ATM0 =',visc_atm0
971  WRITE(*,*)
972  stop
973 
974  END IF
975 
976  IF ( .NOT.isset(rair) ) THEN
977 
978  WRITE(*,*) ''
979  WRITE(*,*) 'ERROR: problem with namelist ATM_PARAMETERS'
980  WRITE(*,*)
981  WRITE(*,atm_parameters)
982  WRITE(*,*)
983  WRITE(*,*) 'Please check RAIR value [J K-1 kg-1]'
984  WRITE(*,*) 'RAIR =',rair
985  WRITE(*,*)
986  stop
987 
988  END IF
989 
990  IF ( .NOT.isset(cpair) ) THEN
991 
992  WRITE(*,*) ''
993  WRITE(*,*) 'ERROR: problem with namelist ATM_PARAMETERS'
994  WRITE(*,*)
995  WRITE(*,atm_parameters)
996  WRITE(*,*)
997  WRITE(*,*) 'Please check CPAIR value [J kg-1 K-1]'
998  WRITE(*,*) 'CPAIR =',cpair
999  WRITE(*,*)
1000  stop
1001 
1002  END IF
1003 
1004  IF ( .NOT.isset(wind_mult_coeff) ) THEN
1005 
1006  WRITE(*,*) ''
1007  WRITE(*,*) 'ERROR: problem with namelist ATM_PARAMETERS'
1008  WRITE(*,*)
1009  WRITE(*,atm_parameters)
1010  WRITE(*,*)
1011  WRITE(*,*) 'Please check WIND_MULT_COEFF value'
1012  WRITE(*,*) 'WIND_MULT_COEFF =',wind_mult_coeff
1013  WRITE(*,*)
1014  stop
1015 
1016  END IF
1017 
1018  WRITE(bak_unit, atm_parameters)
1019  rewind(inp_unit)
1020 
1021  ELSE
1022 
1023  WRITE(*,*) 'Problem with namelist ATM_PARAMETERS'
1024  stop
1025 
1026  END IF
1027 
1028 
1029  IF ( read_atm_profile .EQ. 'table' ) THEN
1030 
1031  n_atm_levels = 0
1032 
1033  READ( inp_unit, table_atm_parameters )
1034  WRITE( bak_unit, table_atm_parameters )
1035 
1036  n_unit = n_unit + 1
1037 
1038  atm_unit = n_unit
1039 
1040  atm_file = '../AtmProfile_info/Density_April.txt'
1041 
1042  OPEN(atm_unit,file=atm_file,status='old')
1043 
1044  READ(atm_unit,*)
1045 
1046  atm_level = 0
1047 
1048  atm_read_levels_apr: DO
1049 
1050  atm_level = atm_level + 1
1051 
1052  READ(atm_unit,*,iostat=io ) rho_atm_apr(atm_level,1:8)
1053 
1054  IF ( io > 0 ) EXIT
1055 
1056  END DO atm_read_levels_apr
1057 
1058  CLOSE(atm_unit)
1059 
1060  n_unit = n_unit + 1
1061 
1062  atm_unit = n_unit
1063 
1064  atm_file = '../AtmProfile_info/Density_Jan.txt'
1065 
1066  OPEN(atm_unit,file=atm_file,status='old')
1067 
1068  READ(atm_unit,*)
1069 
1070  atm_level = 0
1071 
1072  atm_read_levels_jan: DO
1073 
1074  atm_level = atm_level + 1
1075 
1076  READ(atm_unit,*,iostat=io) rho_atm_jan(atm_level,1:8)
1077 
1078  IF ( io > 0 ) EXIT
1079 
1080  END DO atm_read_levels_jan
1081 
1082  CLOSE(atm_unit)
1083 
1084  n_unit = n_unit + 1
1085 
1086  atm_unit = n_unit
1087 
1088  atm_file = '../AtmProfile_info/Density_July.txt'
1089 
1090  OPEN(atm_unit,file=atm_file,status='old')
1091 
1092  READ(atm_unit,*)
1093 
1094  atm_level = 0
1095 
1096  atm_read_levels_jul: DO
1097 
1098  atm_level = atm_level + 1
1099 
1100  READ(atm_unit,*,iostat=io) rho_atm_jul(atm_level,1:8)
1101 
1102  IF ( io > 0 ) EXIT
1103 
1104  END DO atm_read_levels_jul
1105 
1106  CLOSE(atm_unit)
1107 
1108  n_unit = n_unit + 1
1109 
1110  atm_unit = n_unit
1111 
1112  atm_file = '../AtmProfile_info/Density_Oct.txt'
1113 
1114  OPEN(atm_unit,file=atm_file,status='old')
1115 
1116  READ(atm_unit,*)
1117 
1118  atm_level = 0
1119 
1120  atm_read_levels_oct: DO
1121 
1122  atm_level = atm_level + 1
1123 
1124  READ(atm_unit,*,iostat=io) rho_atm_oct(atm_level,1:8)
1125 
1126  IF ( io > 0 ) EXIT
1127 
1128  n_atm_levels = atm_level
1129 
1130  END DO atm_read_levels_oct
1131 
1132  CLOSE(atm_unit)
1133 
1134  ! ----- READ PRESSURES -------
1135 
1136  n_unit = n_unit + 1
1137 
1138  atm_unit = n_unit
1139 
1140  atm_file = '../AtmProfile_info/Pressure_April.txt'
1141 
1142  OPEN(atm_unit,file=atm_file,status='old')
1143 
1144  READ(atm_unit,*)
1145 
1146  atm_level = 0
1147 
1148  pres_read_levels_apr: DO
1149 
1150  atm_level = atm_level + 1
1151 
1152  READ(atm_unit,*,iostat=io) pres_atm_apr(atm_level,1:8)
1153 
1154  IF ( io > 0 ) EXIT
1155 
1156  END DO pres_read_levels_apr
1157 
1158  CLOSE(atm_unit)
1159 
1160  n_unit = n_unit + 1
1161 
1162  atm_unit = n_unit
1163 
1164  atm_file = '../AtmProfile_info/Pressure_Jan.txt'
1165 
1166  OPEN(atm_unit,file=atm_file,status='old')
1167 
1168  READ(atm_unit,*)
1169 
1170  atm_level = 0
1171 
1172  pres_read_levels_jan: DO
1173 
1174  atm_level = atm_level + 1
1175 
1176  READ(atm_unit,*,iostat=io) pres_atm_jan(atm_level,1:8)
1177 
1178  IF ( io > 0 ) EXIT
1179 
1180  END DO pres_read_levels_jan
1181 
1182  CLOSE(atm_unit)
1183 
1184  n_unit = n_unit + 1
1185 
1186  atm_unit = n_unit
1187 
1188  atm_file = '../AtmProfile_info/Pressure_July.txt'
1189 
1190  OPEN(atm_unit,file=atm_file,status='old')
1191 
1192  READ(atm_unit,*)
1193 
1194  atm_level = 0
1195 
1196  pres_read_levels_jul: DO
1197 
1198  atm_level = atm_level + 1
1199 
1200  READ(atm_unit,*,iostat=io) pres_atm_jul(atm_level,1:8)
1201 
1202  IF ( io > 0 ) EXIT
1203 
1204  END DO pres_read_levels_jul
1205 
1206  CLOSE(atm_unit)
1207 
1208 
1209  n_unit = n_unit + 1
1210 
1211  atm_unit = n_unit
1212 
1213  atm_file = '../AtmProfile_info/Pressure_Oct.txt'
1214 
1215  OPEN(atm_unit,file=atm_file,status='old')
1216 
1217  READ(atm_unit,*)
1218 
1219  atm_level = 0
1220 
1221  pres_read_levels_oct: DO
1222 
1223  atm_level = atm_level + 1
1224 
1225  READ(atm_unit,*,iostat=io) pres_atm_oct(atm_level,1:8)
1226 
1227  IF ( io > 0 ) EXIT
1228 
1229  n_atm_levels = atm_level
1230 
1231  END DO pres_read_levels_oct
1232 
1233  CLOSE(atm_unit)
1234 
1235 
1236  ! ----- READ TEMPERATURES -------
1237 
1238  n_unit = n_unit + 1
1239 
1240  atm_unit = n_unit
1241 
1242  atm_file = '../AtmProfile_info/Temp_April.txt'
1243 
1244  OPEN(atm_unit,file=atm_file,status='old')
1245 
1246  READ(atm_unit,*)
1247 
1248  atm_level = 0
1249 
1250  temp_read_levels_apr: DO
1251 
1252  atm_level = atm_level + 1
1253 
1254  READ(atm_unit,*,iostat=io) temp_atm_apr(atm_level,1:8)
1255 
1256  IF ( io > 0 ) EXIT
1257 
1258  END DO temp_read_levels_apr
1259 
1260  CLOSE(atm_unit)
1261 
1262  n_unit = n_unit + 1
1263 
1264  atm_unit = n_unit
1265 
1266  atm_file = '../AtmProfile_info/Temp_Jan.txt'
1267 
1268  OPEN(atm_unit,file=atm_file,status='old')
1269 
1270  READ(atm_unit,*)
1271 
1272  atm_level = 0
1273 
1274  temp_read_levels_jan: DO
1275 
1276  atm_level = atm_level + 1
1277 
1278  READ(atm_unit,*,iostat=io) temp_atm_jan(atm_level,1:8)
1279 
1280  IF ( io > 0 ) EXIT
1281 
1282  END DO temp_read_levels_jan
1283 
1284  CLOSE(atm_unit)
1285 
1286  n_unit = n_unit + 1
1287 
1288  atm_unit = n_unit
1289 
1290  atm_file = '../AtmProfile_info/Temp_July.txt'
1291 
1292  OPEN(atm_unit,file=atm_file,status='old')
1293 
1294  READ(atm_unit,*)
1295 
1296  atm_level = 0
1297 
1298  temp_read_levels_jul: DO
1299 
1300  atm_level = atm_level + 1
1301 
1302  READ(atm_unit,*,iostat=io) temp_atm_jul(atm_level,1:8)
1303 
1304  IF ( io > 0 ) EXIT
1305 
1306  END DO temp_read_levels_jul
1307 
1308  CLOSE(atm_unit)
1309 
1310 
1311  n_unit = n_unit + 1
1312 
1313  atm_unit = n_unit
1314 
1315  atm_file = '../AtmProfile_info/Temp_Oct.txt'
1316 
1317  OPEN(atm_unit,file=atm_file,status='old')
1318 
1319  READ(atm_unit,*)
1320 
1321  atm_level = 0
1322 
1323  temp_read_levels_oct: DO
1324 
1325  atm_level = atm_level + 1
1326 
1327  READ(atm_unit,*,iostat=io) temp_atm_oct(atm_level,1:8)
1328 
1329  IF ( io > 0 ) EXIT
1330 
1331  n_atm_levels = atm_level
1332 
1333  END DO temp_read_levels_oct
1334 
1335  CLOSE(atm_unit)
1336 
1337  ALLOCATE( h_levels(n_atm_levels) )
1338 
1339  ALLOCATE(rho_atm_month_lat(n_atm_levels),rho_atm_month(n_atm_levels,8))
1340  ALLOCATE(pres_atm_month_lat(n_atm_levels),pres_atm_month(n_atm_levels,8))
1341  ALLOCATE(temp_atm_month_lat(n_atm_levels),temp_atm_month(n_atm_levels,8))
1342 
1343  IF ((month .GE. 0.0_wp) .and. (month .LE. 1.0_wp)) THEN
1344  WRITE(*,*) 'winter'
1345  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_jan(1:n_atm_levels,1:8)
1346  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_jan(1:n_atm_levels,1:8)
1347  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_jan(1:n_atm_levels,1:8)
1348 
1349  ELSEIF ((month .GT. 1.0_wp) .and. (month .LE. 2.0_wp)) THEN
1350  WRITE(*,*) 'spring'
1351  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_apr(1:n_atm_levels,1:8)
1352  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_apr(1:n_atm_levels,1:8)
1353  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_apr(1:n_atm_levels,1:8)
1354 
1355  ELSEIF ((month .GT. 2.0_wp) .and. (month .LE. 3.0_wp)) THEN
1356  WRITE(*,*) 'summer'
1357  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_jul(1:n_atm_levels,1:8)
1358  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_jul(1:n_atm_levels,1:8)
1359  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_jul(1:n_atm_levels,1:8)
1360 
1361  ELSEIF ((month .GT. 3.0_wp) .and. (month .LE. 4.0_wp)) THEN
1362  WRITE(*,*) 'autumn'
1363  rho_atm_month(1:n_atm_levels,1:8) = rho_atm_apr(1:n_atm_levels,1:8)
1364  pres_atm_month(1:n_atm_levels,1:8) = pres_atm_apr(1:n_atm_levels,1:8)
1365  temp_atm_month(1:n_atm_levels,1:8) = temp_atm_apr(1:n_atm_levels,1:8)
1366 
1367  END IF
1368 
1369  IF ( ( lat .GE. 0.0_wp ) .AND. ( lat .LE. 15.0_wp ) ) THEN
1370 
1371  coeff_lat = 1.0_wp - ( lat - 0.0_wp ) / ( 15.0_wp - 0.0_wp )
1372 
1373  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1374  rho_atm_month(1:n_atm_levels,2) + ( 1.0_wp - coeff_lat ) * &
1375  rho_atm_month(1:n_atm_levels,3)
1376 
1377  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1378  pres_atm_month(1:n_atm_levels,2) + ( 1.0_wp - coeff_lat ) * &
1379  pres_atm_month(1:n_atm_levels,3)
1380 
1381  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1382  temp_atm_month(1:n_atm_levels,2) + ( 1.0_wp - coeff_lat ) * &
1383  temp_atm_month(1:n_atm_levels,3)
1384 
1385  ELSEIF ( ( lat .GT. 15.0_wp ) .AND. ( lat .LE. 30.0_wp ) ) THEN
1386 
1387  coeff_lat = 1.0_wp - ( lat - 15.0_wp ) / ( 30.0_wp - 15.0_wp )
1388 
1389  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1390  rho_atm_month(1:n_atm_levels,3) + ( 1.0_wp - coeff_lat ) * &
1391  rho_atm_month(1:n_atm_levels,4)
1392 
1393  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1394  pres_atm_month(1:n_atm_levels,3) + ( 1.0_wp - coeff_lat ) * &
1395  pres_atm_month(1:n_atm_levels,5)
1396 
1397  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1398  temp_atm_month(1:n_atm_levels,3) + ( 1.0_wp - coeff_lat ) * &
1399  temp_atm_month(1:n_atm_levels,5)
1400 
1401  ELSEIF ( ( lat .GT. 30.0_wp ) .AND. ( lat .LE. 45.0_wp ) ) THEN
1402 
1403  coeff_lat = 1.0_wp - ( lat - 30.0_wp ) / ( 45.0_wp - 30.0_wp )
1404 
1405  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1406  rho_atm_month(1:n_atm_levels,4) + ( 1.0_wp - coeff_lat ) * &
1407  rho_atm_month(1:n_atm_levels,5)
1408 
1409  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1410  pres_atm_month(1:n_atm_levels,4) + ( 1.0_wp - coeff_lat ) * &
1411  pres_atm_month(1:n_atm_levels,5)
1412 
1413  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1414  temp_atm_month(1:n_atm_levels,4) + ( 1.0_wp - coeff_lat ) * &
1415  temp_atm_month(1:n_atm_levels,5)
1416 
1417  ELSEIF ( ( lat .GT. 45.0_wp ) .AND. ( lat .LE. 60.0_wp ) ) THEN
1418 
1419  coeff_lat = 1.0_wp - ( lat - 45.0_wp ) / ( 60.0_wp - 45.0_wp )
1420 
1421  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1422  rho_atm_month(1:n_atm_levels,5) + ( 1.0_wp - coeff_lat ) * &
1423  rho_atm_month(1:n_atm_levels,6)
1424 
1425  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1426  pres_atm_month(1:n_atm_levels,5) + ( 1.0_wp - coeff_lat ) * &
1427  pres_atm_month(1:n_atm_levels,6)
1428 
1429  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1430  temp_atm_month(1:n_atm_levels,5) + ( 1.0_wp - coeff_lat ) * &
1431  temp_atm_month(1:n_atm_levels,6)
1432 
1433  ELSEIF ( ( lat .GT. 60.0_wp ) .AND. ( lat .LE. 75.0_wp ) ) THEN
1434 
1435  coeff_lat = 1.0_wp - ( lat - 60.0_wp ) / ( 75.0_wp - 60.0_wp )
1436 
1437  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1438  rho_atm_month(1:n_atm_levels,6) + ( 1.0_wp - coeff_lat ) * &
1439  rho_atm_month(1:n_atm_levels,7)
1440 
1441  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1442  pres_atm_month(1:n_atm_levels,6) + ( 1.0_wp - coeff_lat ) * &
1443  pres_atm_month(1:n_atm_levels,7)
1444 
1445  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1446  temp_atm_month(1:n_atm_levels,6) + ( 1.0_wp - coeff_lat ) * &
1447  temp_atm_month(1:n_atm_levels,7)
1448 
1449  ELSEIF ( ( lat .GT. 75.0_wp ) .AND. ( lat .LE. 90.0_wp ) ) THEN
1450 
1451  coeff_lat = 1.0_wp - ( lat - 75.0_wp ) / ( 90.0_wp - 75.0_wp )
1452 
1453  rho_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1454  rho_atm_month(1:n_atm_levels,7) &
1455  + ( 1.0_wp - coeff_lat ) * rho_atm_month(1:n_atm_levels,8)
1456 
1457  pres_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1458  pres_atm_month(1:n_atm_levels,7) &
1459  + ( 1.0_wp - coeff_lat ) * pres_atm_month(1:n_atm_levels,8)
1460 
1461  temp_atm_month_lat(1:n_atm_levels) = coeff_lat * &
1462  temp_atm_month(1:n_atm_levels,7) &
1463  + ( 1.0_wp - coeff_lat ) * temp_atm_month(1:n_atm_levels,8)
1464 
1465  END IF
1466 
1467  pres_atm_month_lat(1:n_atm_levels) = &
1468  100.0_wp * pres_atm_month_lat(1:n_atm_levels)
1469 
1470  h_levels(1:n_atm_levels) = 1000.0_wp * temp_atm_month(1:n_atm_levels,1)
1471 
1472  ELSEIF ( read_atm_profile .EQ. 'card' ) THEN
1473 
1474  tend1 = .false.
1475 
1476  WRITE(*,*) 'search atm_profile'
1477 
1478  atm_profile_search: DO
1479 
1480  READ(inp_unit,*, end = 200 ) card
1481 
1482  IF( trim(card) == 'ATM_PROFILE' ) THEN
1483 
1484  EXIT atm_profile_search
1485 
1486  END IF
1487 
1488  END DO atm_profile_search
1489 
1490  READ(inp_unit,*) n_atm_profile
1491 
1492  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'n_atm_profile',n_atm_profile
1493 
1494  ALLOCATE( atm_profile(7,n_atm_profile) )
1495  ALLOCATE( atm_profile0(7,n_atm_profile) )
1496 
1497  DO i = 1, n_atm_profile
1498 
1499  READ(inp_unit,*) atm_profile0(1:7,i)
1500 
1501  atm_profile(1:7,i) = atm_profile0(1:7,i)
1502  ! convert from km to meters
1503  atm_profile(1,i) = atm_profile(1,i) * 1000.0_wp
1504 
1505  ! convert from hPa to Pa
1506  atm_profile(3,i) = atm_profile(3,i) * 100.0_wp
1507 
1508  atm_profile(6,i) = atm_profile(6,i) * wind_mult_coeff
1509  atm_profile(7,i) = atm_profile(7,i) * wind_mult_coeff
1510 
1511  IF ( verbose_level .GE. 1 ) WRITE(*,*) i,atm_profile(1,i)
1512 
1513  END DO
1514 
1515  GOTO 210
1516 200 tend1 = .true.
1517 210 CONTINUE
1518 
1519  rewind(inp_unit)
1520 
1521  ELSEIF ( read_atm_profile .EQ. 'standard' ) THEN
1522 
1523 
1524  READ( inp_unit,std_atm_parameters,iostat=ios )
1525 
1526  IF ( ios .NE. 0 ) THEN
1527 
1528  WRITE(*,*) 'IOSTAT=',ios
1529  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1530  WRITE(*,*) 'Please check the input file'
1531  WRITE(*,std_atm_parameters)
1532  stop
1533 
1534  ELSE
1535 
1536  IF ( .NOT. isset(u_max) ) THEN
1537 
1538  WRITE(*,*) ''
1539  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1540  WRITE(*,*)
1541  WRITE(*,std_atm_parameters)
1542  WRITE(*,*)
1543  WRITE(*,*) 'Please check u_max value (>0 [m/s])'
1544  WRITE(*,*) 'u_max =',u_max
1545  WRITE(*,*)
1546  stop
1547 
1548  END IF
1549 
1550  IF ( .NOT. isset(sphu_atm0) ) THEN
1551 
1552  IF ( .NOT. isset(rel_hu) ) THEN
1553 
1554  WRITE(*,*) ''
1555  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1556  WRITE(*,*)
1557  WRITE(*,std_atm_parameters)
1558  WRITE(*,*)
1559  WRITE(*,*) 'Please assign sphu_atm0 or rel_hu'
1560  WRITE(*,*)
1561  stop
1562 
1563  ELSEIF ( ( rel_hu .LT. 0.0_wp ) .OR. ( rel_hu .GT. 1.0_wp ) ) THEN
1564 
1565  WRITE(*,*) ''
1566  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1567  WRITE(*,*)
1568  WRITE(*,std_atm_parameters)
1569  WRITE(*,*)
1570  WRITE(*,*) 'Please check rel_hu value (0<=rel_hu<=1)'
1571  WRITE(*,*) 'rel_hu =',rel_hu
1572  WRITE(*,*)
1573  stop
1574 
1575  END IF
1576 
1577  ELSE
1578 
1579  IF ( isset(rel_hu) ) THEN
1580 
1581  WRITE(*,*) ''
1582  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1583  WRITE(*,*)
1584  WRITE(*,std_atm_parameters)
1585  WRITE(*,*)
1586  WRITE(*,*) 'Please assign only sphu_atm0 or rel_hu'
1587  WRITE(*,*)
1588  stop
1589 
1590  ELSE
1591 
1592  IF ( sphu_atm0 .LE. 2.0e-6_wp ) THEN
1593 
1594  WRITE(*,*) 'WARNING: sphu_atm0 value at sea level'
1595  WRITE(*,*) 'should be higher than value at tropopause'
1596  WRITE(*,*) 'base (2.0E-6 kg/kg)'
1597  WRITE(*,*) 'Value changed to 2.01E-6'
1598  WRITE(*,*)
1599  sphu_atm0 = 2.01e-6_wp
1600 
1601  END IF
1602 
1603  END IF
1604 
1605  END IF
1606 
1607  IF ( .NOT. isset(p_atm0) ) THEN
1608 
1609  WRITE(*,*) ''
1610  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1611  WRITE(*,*)
1612  WRITE(*,std_atm_parameters)
1613  WRITE(*,*)
1614  WRITE(*,*) 'Please check p_atm0 value (Pa)'
1615  WRITE(*,*) 'p_atm0 =',p_atm0
1616  WRITE(*,*)
1617  stop
1618 
1619  ELSE
1620 
1621  pa = p_atm0
1622 
1623  END IF
1624 
1625  IF ( .NOT. isset(t_atm0) ) THEN
1626 
1627  WRITE(*,*) ''
1628  WRITE(*,*) 'ERROR: problem with namelist STD_ATM_PARAMETERS'
1629  WRITE(*,*)
1630  WRITE(*,std_atm_parameters)
1631  WRITE(*,*)
1632  WRITE(*,*) 'Please check t_atm0 value (K)'
1633  WRITE(*,*) 't_atm0 =',t_atm0
1634  WRITE(*,*)
1635  stop
1636 
1637  ELSE
1638 
1639  ta = t_atm0
1640 
1641  END IF
1642 
1643 
1644  WRITE(bak_unit, std_atm_parameters)
1645  rewind(inp_unit)
1646 
1647  END IF
1648 
1649  END IF
1650 
1651  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read atm_parameters: done'
1652 
1653  READ(inp_unit,initial_values,iostat=ios)
1654 
1655  IF ( ios .NE. 0 ) THEN
1656 
1657  WRITE(*,*) 'IOSTAT=',ios
1658  WRITE(*,*) 'ERROR: problem with namelist INITIAL_VALUES'
1659  WRITE(*,*) 'Please check the input file'
1660  WRITE(*,initial_values)
1661  stop
1662 
1663  ELSE
1664 
1665  ALLOCATE ( rvolcgas(n_gas) , cpvolcgas(n_gas) , &
1666  volcgas_mass_fraction(n_gas) , volcgas_mol_wt(n_gas) , &
1667  rhovolcgas(n_gas) , volcgas_mass_fraction0(n_gas) )
1668 
1669  WRITE(bak_unit, initial_values)
1670 
1671  rewind(inp_unit)
1672 
1673  END IF
1674 
1675  IF ( inversion_flag ) THEN
1676 
1677  IF ( isset(mfr0) ) THEN
1678 
1679  WRITE(*,*) 'WARNING: you should not assign mfr when inversion is true'
1680  WRITE(*,*) 'in the input file: mfr0',mfr0
1681  stop
1682 
1683  END IF
1684 
1685  IF ( isset(log10_mfr) ) THEN
1686 
1687  WRITE(*,*) 'WARNING: you should not assign mfr when inversion is true'
1688  WRITE(*,*) 'in the input file: log10_mfr',log10_mfr
1689  stop
1690 
1691  END IF
1692 
1693  IF ( isset(r0) .AND. isset(w0) ) THEN
1694 
1695  WRITE(*,*) 'WARNING: you should not assign R0 and W0 with inversion=true'
1696  WRITE(*,*) 'R0',r0
1697  WRITE(*,*) 'W0',w0
1698  stop
1699 
1700  END IF
1701 
1702  IF ( isset(w0) ) THEN
1703 
1704  IF ( isset(w_min) .OR. isset(w_max) ) THEN
1705 
1706  WRITE(*,*) 'WARNING: you should not assign W0,W_MIN and W_MAX with inversion=true'
1707  WRITE(*,*) 'W0',w0
1708  WRITE(*,*) 'W_MIN',w_min
1709  WRITE(*,*) 'W_MAX',w_max
1710  stop
1711 
1712  END IF
1713 
1714  ELSE
1715 
1716  IF ( ( .NOT. isset(w_min) ) .OR. ( w_min .LE. 0 ) ) THEN
1717 
1718  WRITE(*,*) ''
1719  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
1720  WRITE(*,*)
1721  WRITE(*,inversion_parameters)
1722  WRITE(*,*)
1723  WRITE(*,*) 'Please check W_MIN value (>0 [m/s])'
1724  WRITE(*,*) 'W_MIN =',w_min
1725  WRITE(*,*)
1726  stop
1727 
1728  END IF
1729 
1730 
1731  IF ( ( .NOT. isset(w_max) ) .OR. ( w_max .LE. w_min ) ) THEN
1732 
1733  WRITE(*,*) ''
1734  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
1735  WRITE(*,*)
1736  WRITE(*,inversion_parameters)
1737  WRITE(*,*)
1738  WRITE(*,*) 'Please check W_MAX value (>W_MIN [m/s])'
1739  WRITE(*,*) 'W_MAX =',w_max
1740  WRITE(*,*)
1741  stop
1742 
1743  END IF
1744 
1745  END IF
1746 
1747 
1748  IF ( isset(r0) ) THEN
1749 
1750  IF ( isset(r_min) .OR. isset(r_max) ) THEN
1751 
1752  WRITE(*,*) 'WARNING: you should not assign R0,R_MIN and R_MAX with inversion=true'
1753  WRITE(*,*) 'R0',r0
1754  WRITE(*,*) 'R_MIN',r_min
1755  WRITE(*,*) 'R_MAX',r_max
1756  stop
1757 
1758  END IF
1759 
1760  ELSE
1761 
1762  IF ( ( .NOT. isset(r_min) ) .OR. ( r_min .LE. 0 ) ) THEN
1763 
1764  WRITE(*,*) ''
1765  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
1766  WRITE(*,*)
1767  WRITE(*,inversion_parameters)
1768  WRITE(*,*)
1769  WRITE(*,*) 'Please check R_MIN value (>0 [m])'
1770  WRITE(*,*) 'R_MIN =',r_min
1771  WRITE(*,*)
1772  stop
1773 
1774  END IF
1775 
1776  IF ( ( .NOT. isset(r_max) ) .OR. ( r_max .LE. r_min ) ) THEN
1777 
1778  WRITE(*,*) ''
1779  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
1780  WRITE(*,*)
1781  WRITE(*,inversion_parameters)
1782  WRITE(*,*)
1783  WRITE(*,*) 'Please check R_MAX value (>R_MIN [m])'
1784  WRITE(*,*) 'R_MAX =',r_max
1785  WRITE(*,*)
1786  stop
1787 
1788  END IF
1789 
1790  END IF
1791 
1792  IF ( isset(r_min) .AND. isset(w_min) ) THEN
1793 
1794  IF ( n_values .LE. 0 ) THEN
1795 
1796  WRITE(*,*) ''
1797  WRITE(*,*) 'ERROR: problem with namelist INVERSION_PARAMETERS'
1798  WRITE(*,*)
1799  WRITE(*,inversion_parameters)
1800  WRITE(*,*)
1801  WRITE(*,*) 'Please check N_VALUES value (>0 [integer])'
1802  WRITE(*,*) 'N_VALUES =',n_values
1803  WRITE(*,*)
1804  stop
1805 
1806  END IF
1807 
1808  END IF
1809 
1810 
1811  ELSEIF ( isset(mfr0) ) THEN
1812 
1813  IF ( isset(log10_mfr) ) THEN
1814 
1815  WRITE(*,*) 'WARNING: only one of these parameters can be assigned in'
1816  WRITE(*,*) 'the input file: log10_mfr,mfr0',log10_mfr,mfr0
1817  stop
1818 
1819  ELSE
1820 
1821  IF ( .NOT.isset(mfr0) ) THEN
1822 
1823  WRITE(*,*) ''
1824  WRITE(*,*) 'ERROR: problem with namelist INITIAL_VALUES'
1825  WRITE(*,*)
1826  WRITE(*,initial_values)
1827  WRITE(*,*)
1828  WRITE(*,*) 'Please check MFR0 value (>0 [kg/s])'
1829  WRITE(*,*) 'MFR0 =',mfr0
1830  WRITE(*,*)
1831  stop
1832 
1833  ELSE
1834 
1835  log10_mfr = log10(mfr0)
1836 
1837  IF ( write_flag ) WRITE(*,*) 'LOG10 mass eruption rate =',log10_mfr
1838 
1839  END IF
1840 
1841  END IF
1842 
1843  ELSEIF ( isset(w0) ) THEN
1844 
1845  IF ( isset(log10_mfr) .AND. ( .NOT.isset(r0) ) ) THEN
1846 
1847  IF ( write_flag ) WRITE(*,*) &
1848  'WARNING: initial radius calculated from MER and velocity'
1849 
1850  END IF
1851 
1852  IF ( .NOT.isset(log10_mfr) .AND. ( .NOT.isset(r0) ) ) THEN
1853 
1854  WRITE(*,*) ''
1855  WRITE(*,*) 'ERROR: problem with namelist INITIAL_VALUES'
1856  WRITE(*,*)
1857  WRITE(*,initial_values)
1858  WRITE(*,*)
1859  WRITE(*,*) 'Not enough input parameters assigned in INITIAL_VALUES'
1860  WRITE(*,*) 'MFR0',mfr0
1861  WRITE(*,*) 'LOG10_MFR',log10_mfr
1862  WRITE(*,*) 'W0',w0
1863  WRITE(*,*) 'R0',r0
1864  WRITE(*,*)
1865  stop
1866 
1867  END IF
1868 
1869  ELSEIF ( isset(r0) ) THEN
1870 
1871  IF ( isset(log10_mfr) .AND. ( .NOT.isset(w0) ) ) THEN
1872 
1873  IF ( write_flag ) WRITE(*,*) &
1874  'WARNING: initial radius calculated from MER and radius'
1875 
1876  END IF
1877 
1878  IF ( .NOT.isset(log10_mfr) .AND. ( .NOT.isset(w0) ) ) THEN
1879 
1880  WRITE(*,*) ''
1881  WRITE(*,*) 'ERROR: problem with namelist INITIAL_VALUES'
1882  WRITE(*,*)
1883  WRITE(*,initial_values)
1884  WRITE(*,*)
1885  WRITE(*,*) 'Not enough input parameters assigned in INITIAL_VALUES'
1886  WRITE(*,*) 'MFR0',mfr0
1887  WRITE(*,*) 'LOG10_MFR',log10_mfr
1888  WRITE(*,*) 'W0',w0
1889  WRITE(*,*) 'R0',r0
1890  WRITE(*,*)
1891  stop
1892 
1893  END IF
1894 
1895  ELSEIF ( ( .NOT.isset(log10_mfr) ).AND. ( .NOT.isset(r0) ) .AND. &
1896  ( .NOT.isset(w0) ) ) THEN
1897 
1898  WRITE(*,*) ''
1899  WRITE(*,*) 'ERROR: problem with namelist INITIAL_VALUES'
1900  WRITE(*,*)
1901  WRITE(*,initial_values)
1902  WRITE(*,*)
1903  WRITE(*,*) 'Not enough input parameters assigned in INITIAL_VALUES'
1904  WRITE(*,*) 'MFR0',mfr0
1905  WRITE(*,*) 'LOG10_MFR',log10_mfr
1906  WRITE(*,*) 'W0',w0
1907  WRITE(*,*) 'R0',r0
1908  WRITE(*,*)
1909  stop
1910 
1911  END IF
1912 
1913  IF ( isset(log10_mfr) .AND. isset(w0) .AND. isset(r0) ) THEN
1914 
1915  WRITE(*,*) 'ERROR: too many input parameters: input log10_mfr or w0 and r0'
1916  stop
1917 
1918  END IF
1919 
1920  z = vent_height
1921 
1922 
1923  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'read initial_parameters: done'
1924 
1925  ! ----- AGGREGATION
1926  IF ( aggregation_flag ) THEN
1927 
1928  READ(inp_unit, aggregation_parameters,iostat=ios)
1929 
1930  IF ( ios .NE. 0 ) THEN
1931 
1932  WRITE(*,*) 'IOSTAT=',ios
1933  WRITE(*,*) 'ERROR: problem with namelist AGGREGATION_PARAMETERS'
1934  WRITE(*,*) 'Please check the input file'
1935  WRITE(*,aggregation_parameters)
1936  stop
1937 
1938  ELSE
1939 
1940  rewind(inp_unit)
1941 
1942  aggregation_model = lower(aggregation_model)
1943 
1944  IF ( aggregation_model.EQ.'costa') THEN
1945 
1946  IF ( .not.water_flag ) THEN
1947 
1948  WRITE(*,*) ''
1949  WRITE(*,*) 'ERROR: only wet aggregation is possible'
1950  WRITE(*,*) 'with ''Costa'' aggregation model'
1951  WRITE(*,*) 'WATER FLAG =',water_flag
1952 
1953  stop
1954 
1955  END IF
1956 
1957  ELSEIF ( aggregation_model.EQ.'constant') THEN
1958 
1959  IF ( .not.isset(particles_beta0) ) THEN
1960 
1961  WRITE(*,*) ''
1962  WRITE(*,*) 'ERROR: particles_beta0 requested'
1963  WRITE(*,*) 'with ''constant'' aggregation model'
1964  WRITE(*,*) 'PARTICLES_BETA0 = ',particles_beta0
1965 
1966  stop
1967 
1968  END IF
1969 
1970  ELSEIF ( aggregation_model.NE.'brownian') THEN
1971 
1972  WRITE(*,*) 'ERROR: problem with namelist AGGREGATION_PARAMETERS'
1973  WRITE(*,*) 'Please check aggregation_model:',aggregation_model
1974  stop
1975 
1976  END IF
1977 
1978  WRITE(bak_unit, aggregation_parameters)
1979 
1980  END IF
1981 
1982  END IF
1983 
1984  IF ( hysplit_flag ) THEN
1985 
1986  READ(inp_unit,hysplit_parameters,iostat=ios)
1987 
1988  IF ( ios .NE. 0 ) THEN
1989 
1990  WRITE(*,*) 'IOSTAT=',ios
1991  WRITE(*,*) 'ERROR: problem with namelist HYSPLIT_PARAMETERS'
1992  WRITE(*,*) 'Please check the input file'
1993  WRITE(*,hysplit_parameters)
1994  stop
1995 
1996  ELSE
1997 
1998  IF ( .NOT. isset(hy_deltaz) ) THEN
1999 
2000  WRITE(*,*) ''
2001  WRITE(*,*) 'ERROR: problem with namelist HYSPLIT_PARAMETERS'
2002  WRITE(*,*)
2003  WRITE(*,hysplit_parameters)
2004  WRITE(*,*)
2005  WRITE(*,*) 'Please check hy_deltaz value (>0 [m])'
2006  WRITE(*,*) 'hy_deltaz =',hy_deltaz
2007  WRITE(*,*)
2008 
2009  stop
2010 
2011  END IF
2012 
2013  IF ( n_cloud .EQ. -1 ) THEN
2014 
2015  WRITE(*,*) ''
2016  WRITE(*,*) 'ERROR: problem with namelist HYSPLIT_PARAMETERS'
2017  WRITE(*,*) 'Please check n_cloud value (>0 [integer])'
2018  WRITE(*,*) 'n_cloud =',n_cloud
2019 
2020  stop
2021 
2022  END IF
2023 
2024 
2025  rewind(inp_unit)
2026  WRITE(bak_unit, hysplit_parameters)
2027 
2028  ALLOCATE( solid_mfr(n_part) , solid_mfr_old(n_part) )
2029 
2030  hy_z = vent_height + hy_deltaz
2031  hy_z_old = vent_height
2032  hy_x_old = 0.0_wp
2033  hy_y_old = 0.0_wp
2034 
2035  END IF
2036 
2037 
2038  END IF
2039 
2040  ! ---------
2041 
2042  rvolcgas(1:n_gas) = -1.0_wp
2043  cpvolcgas(1:n_gas) = -1.0_wp
2044  volcgas_mol_wt(1:n_gas) = -1.0_wp
2045  volcgas_mass_fraction0(1:n_gas) = -1.0_wp
2046 
2047  IF ( n_gas .GT. 0.0_wp ) THEN
2048 
2049  READ(inp_unit, volcgas_parameters)
2050 
2051  IF ( any( rvolcgas(1:n_gas) ==-1.0_wp ) ) THEN
2052 
2053  WRITE(*,*) 'Error in namelist VOLCGAS PARAMETERS'
2054  WRITE(*,*) 'Please check the values of rvolcgas',rvolcgas(1:n_gas)
2055  stop
2056 
2057  END IF
2058 
2059  IF ( any( cpvolcgas(1:n_gas) .EQ. -1.0_wp ) ) THEN
2060 
2061  WRITE(*,*) 'Error in namelist VOLCGAS PARAMETERS'
2062  WRITE(*,*) 'Please check the values of cpvolcgas',cpvolcgas(1:n_gas)
2063  stop
2064 
2065  END IF
2066 
2067  IF ( any( volcgas_mol_wt(1:n_gas) .EQ. -1.0_wp ) ) THEN
2068 
2069  WRITE(*,*) 'Error in namelist VOLCGAS PARAMETERS'
2070  WRITE(*,*) 'Please check the values of rvolcgas' , &
2071  volcgas_mol_wt(1:n_gas)
2072  stop
2073 
2074  END IF
2075 
2076  IF ( any( volcgas_mass_fraction0(1:n_gas) .EQ. -1.0_wp ) ) THEN
2077 
2078  WRITE(*,*) 'Error in namelist VOLCGAS PARAMETERS'
2079  WRITE(*,*) 'Please check the values of rvolcgas', &
2080  volcgas_mass_fraction0(1:n_gas)
2081  stop
2082 
2083  END IF
2084 
2085 
2086  IF ( ( sum( volcgas_mass_fraction0(1:n_gas) ) + water_mass_fraction0 ) &
2087  .GE. 1.0_wp ) THEN
2088 
2089  WRITE(*,*) 'WARNING: Sum of gas mass fractions :', &
2090  sum( volcgas_mass_fraction0(1:n_part) + water_mass_fraction0 )
2091 
2092  !READ(*,*)
2093 
2094  END IF
2095 
2096  END IF
2097 
2098  rvolcgas_mix = 0.0_wp
2099  cpvolcgas_mix = 0.0_wp
2100  rrhovolcgas_mix = 0.0_wp
2101 
2102  CALL initialize_meteo
2103 
2104  CALL phifromm
2105 
2106  IF ( n_gas .GT. 0 ) THEN
2107 
2108  DO i_gas = 1,n_gas
2109 
2110  rvolcgas_mix = rvolcgas_mix + volcgas_mass_fraction0(i_gas) &
2111  * rvolcgas(i_gas)
2112 
2113  cpvolcgas_mix = cpvolcgas_mix + volcgas_mass_fraction0(i_gas) &
2114  * cpvolcgas(i_gas)
2115 
2116  rrhovolcgas_mix = rrhovolcgas_mix + volcgas_mass_fraction0(i_gas) &
2117  / ( pa / ( rvolcgas(i_gas) * t_mix0 ) )
2118 
2119  END DO
2120 
2121  rvolcgas_mix = rvolcgas_mix / sum( volcgas_mass_fraction0(1:n_gas) )
2122 
2123  cpvolcgas_mix = cpvolcgas_mix / sum( volcgas_mass_fraction0(1:n_gas) )
2124 
2125  rhovolcgas_mix = sum(volcgas_mass_fraction0(1:n_gas)) / rrhovolcgas_mix
2126 
2127  volcgas_mix_mass_fraction = sum(volcgas_mass_fraction0(1:n_gas))
2128 
2129  ELSE
2130 
2131  rvolcgas_mix = 0.0_wp
2132 
2133  cpvolcgas_mix = 0.0_wp
2134 
2135  rhovolcgas_mix = 0.0_wp
2136 
2137  volcgas_mix_mass_fraction = 0.0_wp
2138 
2139  END IF
2140 
2141  IF ( verbose_level .GE. 1 ) THEN
2142 
2143  WRITE(*,*) 'volcgas_mix_mass_fraction',volcgas_mix_mass_fraction
2144 
2145  END IF
2146 
2147  rhowv = pa / ( rwv * t_mix0 )
2148 
2149  ! ---- We assume all volcanic H2O at the vent is water vapor
2150  water_vapor_mass_fraction = water_mass_fraction0
2151 
2152  liquid_water_mass_fraction = 0.0_wp
2153 
2154  gas_mass_fraction = water_vapor_mass_fraction + volcgas_mix_mass_fraction
2155 
2156  IF ( n_gas .GT. 0 ) THEN
2157 
2158  rho_gas = gas_mass_fraction / ( water_vapor_mass_fraction / rhowv &
2159  + volcgas_mix_mass_fraction / rhovolcgas_mix )
2160 
2161  ELSE
2162 
2163  rho_gas = rhowv
2164 
2165  END IF
2166 
2167  IF ( verbose_level .GE. 1 ) THEN
2168 
2169  WRITE(*,*) 'rvolcgas_mix :', rvolcgas_mix
2170  WRITE(*,*) 'cpvolcgas_mix :', cpvolcgas_mix
2171  WRITE(*,*) 'rhovolcgas_mix :', rhovolcgas_mix
2172  WRITE(*,*) 'rhowv :', rhowv
2173  WRITE(*,*) 'rho_gas :', rho_gas
2174  !READ(*,*)
2175 
2176  END IF
2177 
2178  IF ( n_gas .GT. 0 ) WRITE(bak_unit, volcgas_parameters)
2179 
2180  IF ( sum( solid_partial_mass_fraction(1:n_part) ) .NE. 1.0_wp ) THEN
2181 
2182  WRITE(*,*) 'WARNING: Sum of solid mass fractions :', &
2183  sum( solid_partial_mass_fraction(1:n_part) )
2184 
2185  solid_partial_mass_fraction(1:n_part) = &
2186  solid_partial_mass_fraction(1:n_part) &
2187  / sum( solid_partial_mass_fraction(1:n_part) )
2188 
2189  IF ( verbose_level .GE. 1 ) THEN
2190 
2191  WRITE(*,*) ' Modified solid mass fractions :', &
2192  solid_partial_mass_fraction(1:n_part)
2193 
2194  END IF
2195 
2196  END IF
2197 
2198  solid_partial_mass_fraction0 = solid_partial_mass_fraction
2199 
2200  ! solid mass fractions in the mixture
2201  solid_mass_fraction0(1:n_part) = ( 1.0_wp - water_mass_fraction0 &
2202  - volcgas_mix_mass_fraction ) * solid_partial_mass_fraction(1:n_part)
2203 
2204  WRITE(*,*) '---------- INITIALIZATION ------------'
2205  WRITE(*,*)
2206  WRITE(*,*) 'SOLID PARTIAL MASS DISTRIBUTOINS'
2207 
2208  ! loop to initialize the moments of order 1. These values are updated later
2209  ! with the bulk densities of the different particles famlies
2210  DO i_part = 1,n_part
2211 
2212  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'i_part',i_part
2213 
2214  DO i_sect = 1,n_sections
2215 
2216  IF ( distribution .EQ. 'lognormal' ) THEN
2217 
2218  ! evaluate the moments of order 1 (mass) from the parameters of the
2219  ! lognormal distribution
2220  mom0(1,i_sect,i_part) = normpdf(phir(i_sect), mu_lognormal(i_part),&
2221  sigma_lognormal(i_part) )
2222 
2223  ELSEIF ( distribution .EQ. 'bin' ) THEN
2224 
2225  ! assign the moments of order 1 (mass) from the values read in
2226  ! input
2227  mom0(1,i_sect,i_part) = bin_partial_mass_fraction(n_sections-i_sect+1,i_part)
2228 
2229  END IF
2230 
2231  END DO
2232 
2233  mom0(1,:,i_part) = mom0(1,:,i_part) / sum( mom0(1,:,i_part) )
2234 
2235  IF ( verbose_level .GE. 0 ) THEN
2236 
2237  WRITE(*,*) 'Particle phase:',i_part
2238  WRITE(*,"(30F8.2)") phil(n_sections:1:-1)
2239  WRITE(*,"(30F8.2)") phir(n_sections:1:-1)
2240  WRITE(*,"(30ES8.1)") mom0(1,n_sections:1:-1,i_part)
2241  WRITE(*,*)
2242  !READ(*,*)
2243 
2244  END IF
2245 
2246  ! compute the moments of order 0 (number of particles) from the
2247  ! moments of order 1 (mass of particles)
2248  mom0(0,1:n_sections,i_part) = numberfrommass(m(1:n_sections,i_part), &
2249  m(2:n_sections+1,i_part) , mom0(1,1:n_sections,i_part) )
2250 
2251  END DO
2252 
2253  mom = mom0
2254 
2255  CALL eval_quad_values
2256 
2257  DO i_part = 1,n_part
2258 
2259  ! the density of the particles phases are evaluated here. It is
2260  ! independent from the mass fraction of the particles phases, so
2261  ! it is possible to evaluate them with the "uncorrected" moments
2262  rho_solid_avg(i_part) = 1.0_wp / ( sum( f_quad(:,:,i_part) * &
2263  w_quad(:,:,i_part) * m_quad(:,:,i_part) / rho_quad(:,:,i_part) ) / &
2264  sum(f_quad(:,:,i_part) * w_quad(:,:,i_part) * m_quad(:,:,i_part) ) )
2265 
2266  IF ( verbose_level .GE. 1 ) THEN
2267 
2268  WRITE(*,*) 'rho avg',rho_solid_avg(i_part)
2269  READ(*,*)
2270 
2271  END IF
2272 
2273  END DO
2274 
2275  ! the average solid density is evaluated through the mass fractions and
2276  ! the densities of the particles phases
2277  rho_solid_tot_avg = 1.0_wp / sum( solid_partial_mass_fraction(1:n_part) / &
2278  rho_solid_avg(1:n_part) )
2279 
2280  IF ( verbose_level .GE. 1 ) THEN
2281 
2282  WRITE(*,*)
2283  WRITE(*,*) '******* CHECK ON MASS AND VOLUME FRACTIONS *******'
2284  WRITE(*,*) 'rho solid avg', rho_solid_tot_avg
2285 
2286  END IF
2287 
2288  IF ( initial_neutral_density ) THEN
2289 
2290  ! CHECK AND CORRECT
2291 
2292  rho_mix = rho_atm
2293 
2294  solid_tot_volume_fraction0 = ( rho_mix - rho_gas ) / &
2295  ( rho_solid_tot_avg - rho_gas )
2296 
2297  gas_volume_fraction = 1.0_wp - solid_tot_volume_fraction0
2298 
2299  ELSE
2300 
2301  gas_volume_fraction = rho_solid_tot_avg / ( rho_gas * ( 1.0_wp / &
2302  gas_mass_fraction - 1.0_wp ) + rho_solid_tot_avg )
2303 
2304  solid_tot_volume_fraction0 = 1.0_wp - gas_volume_fraction
2305 
2306  rho_mix = gas_volume_fraction * rho_gas + solid_tot_volume_fraction0 &
2307  * rho_solid_tot_avg
2308 
2309  END IF
2310 
2311  IF ( verbose_level .GE. 1 ) THEN
2312 
2313  WRITE(*,*) 'gas_volume_fraction',gas_volume_fraction
2314  WRITE(*,*) 'solid_tot_volume_fraction0',solid_tot_volume_fraction0
2315  WRITE(*,*) 'rho_gas',rho_gas
2316  WRITE(*,*) 'rho_mix',rho_mix
2317 
2318  WRITE(*,*) 'gas_mass_fraction',gas_mass_fraction
2319  WRITE(*,*) 'solid_mass_fractions',solid_mass_fraction0(1:n_part)
2320 
2321  END IF
2322 
2323  DO i_part = 1,n_part
2324 
2325  ! the volume fraction of the particle phases ( with respect to the
2326  ! solid phase only) is evaluated
2327  alfa_s = solid_partial_mass_fraction(i_part) * rho_solid_tot_avg / &
2328  rho_solid_avg(i_part)
2329 
2330  ! this is the volume fraction of the particles phases in the mixture
2331  solid_volume_fraction0(i_part) = solid_tot_volume_fraction0 * alfa_s
2332 
2333  IF ( verbose_level .GE. 1 ) THEN
2334 
2335  WRITE(*,*) 'i_part =',i_part
2336  WRITE(*,*) 'alfa_s',i_part,alfa_s
2337  WRITE(*,*) 'solid_volume_fraction0',solid_volume_fraction0(i_part)
2338  WRITE(*,*) 'solid_partial_mass_fract', &
2339  solid_partial_mass_fraction(i_part)
2340  WRITE(*,*) 'solid_mass_fract', solid_mass_fraction0(i_part)
2341  WRITE(*,*)
2342 
2343  END IF
2344 
2345  END DO
2346 
2347  ! The values of the linear reconstructions at the quadrature points are
2348  ! computed with the corrected values of the moments
2349  CALL eval_quad_values
2350 
2351  IF ( verbose_level .GE. 1 ) THEN
2352 
2353  WRITE(*,*) 'gas volume fraction', gas_volume_fraction
2354  WRITE(*,*) 'gas mass fraction', gas_mass_fraction
2355 
2356  END IF
2357 
2358  IF ( umbrella_flag ) THEN
2359 
2360  nbl_stop = .true.
2361  WRITE(*,*) 'Plume equations integrated up to neutral buoyance level'
2362 
2363  ! ------- READ run_parameters NAMELIST -----------------------------------
2364  READ(inp_unit, umbrella_run_parameters,iostat=ios )
2365 
2366  IF ( ios .NE. 0 ) THEN
2367 
2368  WRITE(*,*) 'IOSTAT=',ios
2369  WRITE(*,*) 'ERROR: problem with namelist UMBRELLA_RUN_PARAMETERS'
2370  WRITE(*,*) 'Please check the input file'
2371  WRITE(*,umbrella_run_parameters)
2372  stop
2373 
2374  ELSE
2375 
2376  IF ( ( .NOT.isset(c_d) ) .OR. ( c_d .LT. 0.0_wp ) ) THEN
2377 
2378  WRITE(*,*) ''
2379  WRITE(*,*) 'ERROR: problem with namelist UMBRELLA_RUN_PARAMETERS'
2380  WRITE(*,*)
2381  WRITE(*,umbrella_run_parameters)
2382  WRITE(*,*)
2383  WRITE(*,*) 'Please check C_D value'
2384  WRITE(*,*) 'C_D =',c_d
2385  WRITE(*,*)
2386  stop
2387 
2388  END IF
2389 
2390  rewind(inp_unit)
2391  WRITE(bak_unit, umbrella_run_parameters)
2392 
2393  END IF
2394 
2395  ! ------- READ numeric_parameters NAMELIST ----------------------------------
2396 
2397  READ(inp_unit,numeric_parameters,iostat=ios )
2398 
2399  IF ( ios .NE. 0 ) THEN
2400 
2401  WRITE(*,*) 'IOSTAT=',ios
2402  WRITE(*,*) 'ERROR: problem with namelist NUMERIC_PARAMETERS'
2403  WRITE(*,*) 'Please check the input file'
2404  WRITE(*,numeric_parameters)
2405  stop
2406 
2407  ELSE
2408 
2409  rewind(inp_unit)
2410  WRITE(bak_unit, numeric_parameters)
2411 
2412  END IF
2413 
2414  IF ( ( solver_scheme .NE. 'LxF' ) .AND. ( solver_scheme .NE. 'KT' ) &
2415  .AND. ( solver_scheme .NE. 'GFORCE' ) &
2416  .AND. ( solver_scheme .NE. 'UP' ) ) THEN
2417 
2418  WRITE(*,*) 'WARNING: no correct solver scheme selected',solver_scheme
2419  WRITE(*,*) 'Chose between: LxF, GFORCE or KT'
2420  stop
2421 
2422  END IF
2423 
2424  IF ( ( cfl .GT. 0.25 ) .OR. ( cfl .LT. 0.0_wp ) ) THEN
2425 
2426  WRITE(*,*) 'WARNING: wrong value of cfl ',cfl
2427  WRITE(*,*) 'Choose a value between 0.0 and 0.25'
2428  READ(*,*)
2429 
2430  END IF
2431 
2432  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'Limiters',limiter(1:n_vars)
2433 
2434  limiter(n_vars+1) = limiter(2)
2435  limiter(n_vars+2) = limiter(3)
2436 
2437  IF ( ( maxval(limiter(1:n_vars)) .GT. 3 ) .OR. &
2438  ( minval(limiter(1:n_vars)) .LT. 0 ) ) THEN
2439 
2440  WRITE(*,*) 'WARNING: wrong limiter ',limiter(1:n_vars)
2441  WRITE(*,*) 'Choose among: none, minmod,superbee,van_leer'
2442  stop
2443 
2444  END IF
2445 
2446  IF ( verbose_level .GE. 0 ) THEN
2447 
2448  WRITE(*,*) 'Linear reconstruction and b. c. applied to variables:'
2449  WRITE(*,*) 'h,hu,hv'
2450 
2451  END IF
2452 
2453  IF ( ( reconstr_coeff .GT. 1.0_wp ) .OR. ( reconstr_coeff .LT. 0.0_wp ) )&
2454  THEN
2455 
2456  WRITE(*,*) 'WARNING: wrong value of reconstr_coeff ',reconstr_coeff
2457  WRITE(*,*) 'Change the value between 0.0 and 1.0 in the input file'
2458  READ(*,*)
2459 
2460  END IF
2461 
2462  END IF
2463  ! ---------------------------------------------------------------------------
2464 
2465 
2466  ! Close input file
2467 
2468  CLOSE(inp_unit)
2469 
2470  IF ( read_atm_profile .EQ. 'card' ) THEN
2471 
2472  WRITE(bak_unit,*) '''ATM_PROFILE'''
2473  WRITE(bak_unit,*) n_atm_profile
2474 
2475  DO i = 1, n_atm_profile
2476 
2477  WRITE(bak_unit,107) atm_profile0(1:7,i)
2478 
2479 107 FORMAT(7(1x,es14.7))
2480 
2481 
2482  END DO
2483 
2484  END IF
2485 
2486  CLOSE(bak_unit)
2487 
2488  IF ( verbose_level .GE. 1 ) WRITE(*,*) 'end subroutine reainp'
2489 
2490  RETURN
2491 
2492  END SUBROUTINE read_inp
2493 
2494  !******************************************************************************
2496  !
2501  !******************************************************************************
2502 
2503  SUBROUTINE open_file_units
2504 
2505  ! External variables
2506  USE particles_module, ONLY : n_part
2507  USE moments_module, ONLY : n_mom
2508  USE variables, ONLY : dakota_flag , hysplit_flag
2509 
2510  IMPLICIT NONE
2511 
2512  col_file = trim(run_name)//'.col'
2513  sed_file = trim(run_name)//'.sed'
2514  mom_file = trim(run_name)//'.mom'
2515  dak_file = trim(run_name)//'.dak'
2516  hy_file = trim(run_name)//'.hy'
2517  hy_file_volcgas = trim(run_name)//'_volcgas.hy'
2518  inversion_file = trim(run_name)//'.inv'
2519 
2520 
2521  n_unit = n_unit + 1
2522  col_unit = n_unit
2523 
2524  OPEN(col_unit,file=col_file)
2525 
2526  n_unit = n_unit + 1
2527  sed_unit = n_unit
2528 
2529  OPEN(sed_unit,file=sed_file)
2530 
2531  n_unit = n_unit + 1
2532  mom_unit = n_unit
2533 
2534  OPEN(mom_unit,file=mom_file)
2535 
2536  WRITE(mom_unit,*) n_part
2537  WRITE(mom_unit,*) n_mom
2538 
2539 
2540  n_unit = n_unit + 1
2541  hy_unit = n_unit
2542 
2543  IF ( hysplit_flag ) THEN
2544 
2545  OPEN(hy_unit,file=hy_file)
2546 
2547  END IF
2548 
2549  n_unit = n_unit + 1
2550  dak_unit = n_unit
2551 
2552  OPEN(dak_unit,file=dak_file)
2553 
2554  IF ( inversion_flag ) THEN
2555 
2556  n_unit = n_unit + 1
2557  inversion_unit = n_unit
2558 
2559  OPEN(inversion_unit,file=inversion_file)
2560  WRITE(inversion_unit,187)
2561 187 FORMAT(1x,' radius (m) ',1x,' velocity (m/s) ',1x, &
2562  'MER (kg/s) ', 1x,'plume height (m)',1x, &
2563  ' inversion ',1x,'column regime')
2564 
2565  END IF
2566 
2567  RETURN
2568 
2569  END SUBROUTINE open_file_units
2570 
2571  !******************************************************************************
2573  !
2578  !******************************************************************************
2579 
2580  SUBROUTINE close_file_units
2581 
2582  USE variables, ONLY : dakota_flag , hysplit_flag
2583 
2584  IMPLICIT NONE
2585 
2586  CLOSE(col_unit)
2587  CLOSE(sed_unit)
2588  CLOSE(mom_unit)
2589 
2590  IF ( hysplit_flag ) CLOSE ( hy_unit )
2591 
2592  IF ( .NOT. umbrella_flag ) CLOSE(dak_unit)
2593 
2594  IF ( inversion_flag ) CLOSE ( inversion_unit )
2595 
2596  RETURN
2597 
2598  END SUBROUTINE close_file_units
2599 
2600  !******************************************************************************
2602  !
2608  !******************************************************************************
2609 
2610  SUBROUTINE write_column
2611 
2612  USE meteo_module, ONLY: rho_atm , ta, pa
2613 
2614  USE particles_module, ONLY: n_mom , n_part , solid_partial_mass_fraction , &
2615  mom , cum_particle_loss_rate
2616 
2617  USE plume_module, ONLY: x , y , z , w , r , mag_u
2618 
2619  USE mixture_module, ONLY: rho_mix , t_mix , atm_mass_fraction , &
2620  volcgas_mix_mass_fraction , volcgas_mass_fraction, &
2621  dry_air_mass_fraction , water_vapor_mass_fraction , &
2622  liquid_water_mass_fraction, ice_mass_fraction
2623 
2624  USE variables, ONLY: verbose_level
2625 
2626  IMPLICIT NONE
2627 
2628  REAL(wp) :: mfr
2629 
2630  INTEGER :: i_part , j_part , i_mom , i_sect
2631 
2632  INTEGER :: i_gas
2633 
2634  CHARACTER(15) :: mom_str
2635  CHARACTER(LEN=2) :: i_part_string , i_sect_string
2636 
2637  mfr = 3.14 * r**2 * rho_mix * w
2638 
2639  ! WRITE(*,*) 'INPOUT: atm_mass_fraction',atm_mass_fraction
2640  ! READ(*,*)
2641 
2642  IF ( z .EQ. vent_height ) THEN
2643 
2644  col_lines = 0
2645 
2646  WRITE(col_unit,97,advance="no")
2647  WRITE(sed_unit,96,advance="no")
2648 
2649  DO i_part=1,n_part
2650 
2651  DO i_sect=1,n_sections
2652 
2653  i_part_string = lettera(i_part)
2654  i_sect_string = lettera(i_sect)
2655  mom_str = ' rhoB'//i_part_string//'_'//i_sect_string
2656 
2657  WRITE(col_unit,98,advance="no") mom_str
2658  WRITE(sed_unit,98,advance="no") mom_str
2659 
2660  END DO
2661 
2662  END DO
2663 
2664  DO i_gas=1,n_gas
2665 
2666  WRITE(col_unit,99,advance="no")
2667 
2668  END DO
2669 
2670  WRITE(col_unit,100)
2671 
2672  WRITE(sed_unit,*) ''
2673 
2674 96 FORMAT(1x,' z(m) ',1x,' r(m) ',1x,' x(m) ', &
2675  1x,' y(m) ')
2676 
2677 97 FORMAT(1x,' z(m) ',1x,' r(m) ',1x,' x(m) ', &
2678  1x,' y(m) ',1x,'mix.dens(kg/m3)',1x,'temperature(C)', &
2679  1x,' vert.vel.(m/s)',1x,' mag.vel.(m/s) ',1x,' d.a.massfract ', &
2680  1x,' w.v.massfract ',1x,' l.w.massfract ',1x' i.massfract ',1x)
2681 
2682 
2683 
2684 98 FORMAT(1x,a15)
2685 198 FORMAT(1x,'agr.massfract ')
2686 
2687 99 FORMAT(1x,' volgas.massf ')
2688 
2689 100 FORMAT(1x,' volgasmix.massf',1x,'atm.rho(kg/m3)',1x,' MFR(kg/s) ', &
2690  1x,'atm.temp(K) ', 1x,' atm.pres.(Pa) ')
2691 
2692 
2693  END IF
2694 
2695  col_lines = col_lines + 1
2696 
2697  WRITE(col_unit,101,advance="no") z , r , x , y , rho_mix , t_mix-273.15_wp ,&
2698  w , mag_u, dry_air_mass_fraction , water_vapor_mass_fraction , &
2699  liquid_water_mass_fraction , ice_mass_fraction
2700 
2701  WRITE(sed_unit,141,advance="no") z , r , x , y
2702 
2703 101 FORMAT(12(1x,es16.9))
2704 141 FORMAT(4(1x,es16.9))
2705 
2706  DO i_part=1,n_part
2707 
2708  DO i_sect=1,n_sections
2709 
2710  WRITE(col_unit,102,advance="no") mom(1,i_sect,i_part)
2711  WRITE(sed_unit,102,advance="no") abs(cum_particle_loss_rate(i_part,i_sect))
2712 
2713  END DO
2714 
2715  END DO
2716 
2717 102 FORMAT(1(1x,es16.9))
2718 
2719 
2720  WRITE(col_unit,103) volcgas_mass_fraction(1:n_gas) , &
2721  volcgas_mix_mass_fraction , rho_atm , mfr , ta, pa
2722 
2723  WRITE(sed_unit,*) ''
2724 
2725 
2726 103 FORMAT(20(1x,es16.9))
2727 
2728  WRITE(mom_unit,104,advance="no") z
2729 
2730 104 FORMAT(1(1x,es15.8))
2731 
2732  DO i_mom=0,n_mom-1
2733 
2734  DO i_sect=1,n_sections
2735 
2736  DO i_part=1,n_part
2737 
2738  WRITE(mom_unit,105,advance="no") mom(i_mom,i_sect,i_part)
2739 
2740  END DO
2741 
2742  END DO
2743 
2744  END DO
2745 
2746 105 FORMAT(1(1x,es16.9))
2747 
2748  WRITE(mom_unit,*) " "
2749 
2750  IF ( verbose_level .GE. 1 ) THEN
2751 
2752  WRITE(*,*) '******************'
2753  WRITE(*,*) 'z',z
2754  WRITE(*,*) 'x',x
2755  WRITE(*,*) 'y',y
2756  WRITE(*,*) 'r',r
2757  WRITE(*,*) 'w',w
2758  WRITE(*,*) '******************'
2759 
2760  END IF
2761 
2762  RETURN
2763 
2764  END SUBROUTINE write_column
2765 
2766  !******************************************************************************
2768  !
2776  !******************************************************************************
2777 
2778  SUBROUTINE write_dakota(description,value)
2779 
2780  USE variables, ONLY : verbose_level
2781 
2782  IMPLICIT NONE
2783 
2784  CHARACTER(20), INTENT(IN) :: description
2785 
2786  REAL(wp), INTENT(IN) :: value
2787 
2788  WRITE(dak_unit,*) description,value
2789 
2790  IF ( verbose_level .GE. 2 ) THEN
2791 
2792  WRITE(*,*) description,value
2793 
2794  END IF
2795 
2796  RETURN
2797 
2798  END SUBROUTINE write_dakota
2799 
2800  !******************************************************************************
2802  !
2805  !
2815  !******************************************************************************
2816 
2817  SUBROUTINE write_inversion(r0,w_opt,opt_mfr,opt_height,search_flag, &
2818  opt_regime)
2819 
2820  REAL(wp),INTENT(IN) :: r0,w_opt,opt_mfr,opt_height
2821  LOGICAL,INTENT(IN) :: search_flag
2822  INTEGER,INTENT(IN) :: opt_regime
2823 
2824  WRITE(inversion_unit,181) r0,w_opt,opt_mfr,opt_height,search_flag,opt_regime
2825 
2826 181 FORMAT(2(2x,f15.8),1(1x,es15.8),1(1x,f15.6)4x,l,7x,i4)
2827 
2828  RETURN
2829 
2830  END SUBROUTINE write_inversion
2831 
2832  !******************************************************************************
2834  !
2840  !******************************************************************************
2841 
2842  SUBROUTINE write_zero_hysplit
2843 
2844  USE particles_module, ONLY: n_part
2845 
2846  IMPLICIT NONE
2847 
2848  CHARACTER(len=8) :: x1 , x2 ! format descriptor
2849 
2850  INTEGER :: i_part , i_sect
2851 
2852  REAL(wp), ALLOCATABLE :: delta_solid(:)
2853 
2854  INTEGER :: n_tot
2855 
2856  n_tot = n_part * n_sections
2857 
2858  OPEN(hy_unit,file=hy_file)
2859 
2860  WRITE(hy_unit,107,advance="no")
2861 
2862  DO i_part=1,n_part
2863 
2864  WRITE(x1,'(I2.2)') i_part ! converting integer to string
2865 
2866  DO i_sect=1,n_sections
2867 
2868  WRITE(x2,'(I2.2)') i_sect ! converting integer to string
2869 
2870  WRITE(hy_unit,108,advance="no")'S_'//trim(x1)//'_'//trim(x2)//' (kg/s)'
2871 
2872  END DO
2873 
2874  END DO
2875 
2876  WRITE(hy_unit,*) ''
2877 
2878  ALLOCATE( delta_solid(n_tot) )
2879 
2880  delta_solid(1:n_part) = 0.0_wp
2881 
2882  WRITE(hy_unit,110) 0.0_wp , 0.0_wp , vent_height + 0.50_wp * hy_deltaz , &
2883  delta_solid(1:n_part)
2884 
2885  DEALLOCATE( delta_solid )
2886 
2887  CLOSE(hy_unit)
2888 
2889 107 FORMAT(1x,' x (m) ',1x,' y (m) ', 1x,' z (m) ')
2890 
2891 108 FORMAT(2x,a)
2892 
2893 110 FORMAT(50(1x,e15.8))
2894 
2895  END SUBROUTINE write_zero_hysplit
2896 
2897  !******************************************************************************
2899  !
2905  !******************************************************************************
2906 
2907  SUBROUTINE check_hysplit
2908 
2909  USE meteo_module, ONLY: rho_atm , ta, pa , interp_1d_scalar
2910  USE meteo_module, ONLY : cos_theta , sin_theta , u_atm , zmet
2911 
2912  USE particles_module, ONLY: n_mom , n_part , solid_partial_mass_fraction , &
2913  mom
2914 
2915  USE plume_module, ONLY: x , y , z , w , r , mag_u
2916 
2917  USE mixture_module, ONLY: rho_mix , t_mix , atm_mass_fraction , &
2918  volcgas_mix_mass_fraction , volcgas_mass_fraction, &
2919  dry_air_mass_fraction , water_vapor_mass_fraction , &
2920  liquid_water_mass_fraction, ice_mass_fraction
2921 
2922  USE variables, ONLY : height_nbl
2923 
2924  IMPLICIT NONE
2925 
2926  CHARACTER(len=8) :: x1 , x2 ! format descriptor
2927 
2928  INTEGER :: i , j , n_hy
2929 
2930  REAL(wp) :: temp_k,mfr
2931  REAL(wp) :: da_mf,wv_mf,lw_mf, ice_mf, volcgas_tot_mf
2932  REAL(wp), ALLOCATABLE :: x_col(:) , y_col(:) , z_col(:) , r_col(:)
2933  REAL(wp), ALLOCATABLE :: mom_col(:,:) , mfr_col(:)
2934  REAL(wp), ALLOCATABLE :: volcgas_mf(:,:)
2935  REAL(wp), ALLOCATABLE :: solid_mass_flux(:,:) , solid_mass_loss_cum(:,:)
2936  REAL(wp), ALLOCATABLE :: volcgas_mass_flux(:,:)
2937  REAL(wp) :: z_min , z_max , z_bot , z_top , x_top , x_bot , y_bot , y_top
2938  REAL(wp) :: r_bot , r_top
2939  REAL(wp) :: solid_bot , solid_top
2940  REAL(wp) :: solid_loss_bot , solid_loss_top
2941  REAL(wp) :: gas_top
2942  REAL(wp), ALLOCATABLE :: delta_solid(:) , cloud_solid(:)
2943  REAL(wp), ALLOCATABLE :: cloud_gas(:)
2944  REAL(wp), ALLOCATABLE :: solid_tot(:)
2945 
2946 
2947  REAL(wp) :: angle_release , start_angle
2948  REAL(wp) :: delta_angle
2949  REAL(wp) :: dx , dy , dz , dv(3)
2950 
2951  REAL(wp) :: vect(3) , vect0(3) , v(3) , c , s
2952  REAL(wp) :: mat_v(3,3) , mat_R(3,3)
2953 
2954  INTEGER :: i_part , i_sect
2955  INTEGER :: n_tot
2956 
2957  INTEGER :: read_col_unit , read_sed_unit
2958 
2959 
2960  n_tot = n_part * n_sections
2961 
2962  ALLOCATE( x_col(col_lines) , y_col(col_lines) , z_col(col_lines) )
2963  ALLOCATE( r_col(col_lines) )
2964  ALLOCATE( mom_col(n_tot,col_lines) )
2965  ALLOCATE( mfr_col(col_lines) )
2966  ALLOCATE( volcgas_mf(n_gas,col_lines) )
2967  ALLOCATE( solid_mass_flux(n_tot,col_lines) )
2968  ALLOCATE( solid_mass_loss_cum(n_tot,col_lines) )
2969  ALLOCATE( volcgas_mass_flux(n_gas,col_lines) )
2970  ALLOCATE( delta_solid(n_tot) , cloud_solid(n_tot) )
2971  ALLOCATE( cloud_gas(n_gas) )
2972  ALLOCATE( solid_tot(n_tot) )
2973 
2974  n_unit = n_unit + 1
2975  read_col_unit = n_unit
2976 
2977  OPEN(read_col_unit,file=col_file)
2978 
2979  READ(read_col_unit,*)
2980 
2981  DO i = 1,col_lines
2982 
2983  READ(read_col_unit,111) z_col(i) , r_col(i) , x_col(i) , y_col(i) , &
2984  rho_mix , temp_k , w , mag_u, da_mf , wv_mf , lw_mf , ice_mf , &
2985  mom_col(1:n_tot,i) , volcgas_mf(1:n_gas,i) , volcgas_tot_mf , &
2986  rho_atm , mfr_col(i) , ta, pa
2987 
2988  solid_mass_flux(1:n_tot,i) = mom_col(1:n_tot,i) * pi_g * r_col(i)**2 &
2989  * w
2990 
2991  volcgas_mass_flux(1:n_gas,i) = volcgas_mf(1:n_gas,i) &
2992  * rho_mix * pi_g * r_col(i)**2 * w
2993 
2994  !WRITE(*,*) 'Solid mass flux (kg/s): ',solid_mass_flux(1:n_tot,i)
2995  !WRITE(*,*) 'Total solid mass flux (kg/s): ',SUM(solid_mass_flux(1:n_tot,i))
2996  !WRITE(*,*) 'solid_pmf: ',solid_pmf(1:n_tot,i)
2997  !WRITE(*,*) 'Sum solid mass fractions: ',SUM(solid_pmf(1:n_tot,i))
2998  !WRITE(*,*) z_col(i) , solid_mass_loss_cum(1:n_tot,i)
2999  !READ(*,*)
3000  !WRITE(*,*) 'volcgas_mass_flux ',volcgas_mass_flux(1:n_gas,i), z_col(i)
3001  !READ(*,*)
3002 
3003  END DO
3004 
3005 111 FORMAT(90(1x,es16.9))
3006 
3007  CLOSE(read_col_unit)
3008 
3009  n_unit = n_unit + 1
3010  read_sed_unit = n_unit
3011 
3012  OPEN(read_sed_unit,file=sed_file)
3013 
3014  READ(read_sed_unit,*)
3015 
3016  DO i = 1,col_lines
3017 
3018  READ(read_sed_unit,112) z_col(i) , r_col(i) , x_col(i) , y_col(i) , &
3019  solid_mass_loss_cum(1:n_tot,i)
3020 
3021  END DO
3022 
3023 112 FORMAT(200(1x,es16.9))
3024 
3025  CLOSE(read_sed_unit)
3026 
3027  OPEN(hy_unit,file=hy_file)
3028 
3029  WRITE(hy_unit,107,advance="no")
3030 
3031  DO i_part=1,n_part
3032 
3033  WRITE(x1,'(I2.2)') i_part ! converting integer to string
3034 
3035  DO i_sect=1,n_sections
3036 
3037  WRITE(x2,'(I2.2)') i_sect ! converting integer to string
3038 
3039  WRITE(hy_unit,108,advance="no")'S_'//trim(x1)//'_'//trim(x2)//' (kg/s)'
3040 
3041  END DO
3042 
3043  END DO
3044 
3045  WRITE(hy_unit,*) ''
3046 
3047  z_min = z_col(1)
3048 
3049  IF ( nbl_stop ) THEN
3050 
3051  z_max = height_nbl + z_min
3052 
3053  ELSE
3054 
3055  z_max = z_col(col_lines)
3056 
3057  END IF
3058 
3059  n_hy = floor( ( z_max - z_min ) / hy_deltaz )
3060 
3061  solid_tot(1:n_tot) = 0.0_wp
3062 
3063  DO i = 1,n_hy
3064 
3065  z_bot = z_min + (i-1) * hy_deltaz
3066  z_top = z_min + i * hy_deltaz
3067 
3068  z = z_bot
3069 
3070  DO j = 1,n_tot
3071 
3072  CALL interp_1d_scalar(z_col, x_col, z_bot, x_bot)
3073  CALL interp_1d_scalar(z_col, x_col, z_top, x_top)
3074 
3075  CALL interp_1d_scalar(z_col, x_col, z_bot, y_bot)
3076  CALL interp_1d_scalar(z_col, x_col, z_top, y_top)
3077 
3078  CALL interp_1d_scalar(z_col, r_col, z_bot, r_bot)
3079  CALL interp_1d_scalar(z_col, r_col, z_top, r_top)
3080 
3081  ! CALL interp_1d_scalar(z_col, solid_mass_flux(j,:), z_bot, solid_bot)
3082  ! CALL interp_1d_scalar(z_col, solid_mass_flux(j,:), z_top, solid_top)
3083  ! delta_solid(j) = solid_bot - solid_top
3084 
3085  CALL interp_1d_scalar(z_col, solid_mass_loss_cum(j,:), z_bot, solid_loss_bot)
3086  CALL interp_1d_scalar(z_col, solid_mass_loss_cum(j,:), z_top, solid_loss_top)
3087 
3088  delta_solid(j) = abs(solid_loss_top - solid_loss_bot)
3089 
3090  END DO
3091 
3092  IF ( n_cloud .EQ. 1 ) THEN
3093 
3094  IF ( verbose_level .GE. 1 ) THEN
3095 
3096  WRITE(*,110) 0.5_wp * ( x_top + x_bot ) , 0.5_wp * ( y_top+y_bot ) , &
3097  0.5_wp * ( z_top + z_bot ) , delta_solid(1:n_part)
3098 
3099  !READ(*,*)
3100 
3101  END IF
3102 
3103  WRITE(hy_unit,110) 0.5_wp * ( x_top+x_bot ) , 0.5_wp * ( y_top+y_bot ) ,&
3104  0.5_wp * ( z_top + z_bot ) , delta_solid(1:n_tot)
3105 
3106  ELSE
3107 
3108  CALL zmet
3109 
3110  IF ( u_atm .LT. 1.0d+3 ) THEN
3111 
3112  delta_angle = 2.0_wp*pi_g/n_cloud
3113 
3114  ELSE
3115 
3116  delta_angle = pi_g / ( n_cloud - 1.0_wp )
3117 
3118  END IF
3119 
3120 
3121  DO j=1,n_cloud
3122 
3123  start_angle = atan2(sin_theta,cos_theta)
3124  angle_release = (j-1) * delta_angle - 0.5_wp*pi_g
3125 
3126  dx = 0.5_wp* ( r_bot + r_top ) * cos(start_angle + angle_release)
3127  dy = 0.5_wp* ( r_bot + r_top ) * sin(start_angle + angle_release)
3128  dz = 0.0_wp
3129 
3130  !WRITE(*,*) "dx,dy ",dx,dy
3131  !WRITE(*,*) "start_angle ",start_angle
3132  !WRITE(*,*) "angle_release ",angle_release
3133  !WRITE(*,*) "delta_angle ",delta_angle
3134  !READ(*,*)
3135 
3136 
3137  IF ( verbose_level .GE. 1 ) THEN
3138 
3139  WRITE(*,110) 0.5_wp * ( x_top + x_bot ) + dx , &
3140  0.5_wp * ( y_top + y_bot ) + dy , &
3141  0.5_wp * ( z_top + z_bot ) + dz , &
3142  delta_solid(1:n_tot)/n_cloud
3143 
3144  END IF
3145 
3146  WRITE(hy_unit,110) 0.5_wp * ( x_top + x_bot ) + dx , &
3147  0.5_wp * ( y_top + y_bot ) + dy , &
3148  0.5_wp * ( z_top + z_bot ) + dz , &
3149  delta_solid(1:n_tot)/n_cloud
3150 
3151  END DO
3152 
3153  END IF
3154 
3155  solid_tot(1:n_tot) = solid_tot(1:n_tot) + delta_solid(1:n_tot)
3156 
3157  END DO
3158 
3159  ! WRITE THE RELEASE FROM THE MIDDLE OF LAST INTERVAL
3160 
3161  z_bot = z_min + n_hy * hy_deltaz
3162  z_top = z_max
3163 
3164  DO j = 1,n_tot
3165 
3166 
3167  CALL interp_1d_scalar(z_col, x_col, z_bot, x_bot)
3168  CALL interp_1d_scalar(z_col, x_col, z_top, x_top)
3169 
3170  CALL interp_1d_scalar(z_col, x_col, z_bot, y_bot)
3171  CALL interp_1d_scalar(z_col, x_col, z_top, y_top)
3172 
3173  CALL interp_1d_scalar(z_col, r_col, z_bot, r_bot)
3174  CALL interp_1d_scalar(z_col, r_col, z_top, r_top)
3175 
3176  ! CALL interp_1d_scalar(z_col, solid_mass_flux(j,:), z_bot, solid_bot)
3177  CALL interp_1d_scalar(z_col, solid_mass_flux(j,:), z_top, solid_top)
3178  ! delta_solid(j) = solid_bot - solid_top
3179 
3180  cloud_solid(j) = solid_top
3181 
3182  CALL interp_1d_scalar(z_col, solid_mass_loss_cum(j,:), z_bot, solid_loss_bot)
3183  CALL interp_1d_scalar(z_col, solid_mass_loss_cum(j,:), z_top, solid_loss_top)
3184 
3185  delta_solid(j) = solid_loss_top - solid_loss_bot
3186 
3187  END DO
3188 
3189  solid_tot(1:n_tot) = solid_tot(1:n_tot) + delta_solid(1:n_tot)
3190  solid_tot(1:n_tot) = solid_tot(1:n_tot) + cloud_solid(1:n_tot)
3191 
3192  IF ( n_cloud .EQ. 1 ) THEN
3193 
3194  IF ( verbose_level .GE. 1 ) THEN
3195 
3196  WRITE(*,110) 0.5_wp * ( x_top + x_bot ) , 0.5_wp * ( y_top + y_bot ) , &
3197  0.5_wp * ( z_top + z_bot ) , delta_solid(1:n_tot)
3198 
3199  END IF
3200 
3201  WRITE(hy_unit,110) 0.5_wp * ( x_top + x_bot ) , 0.5_wp * ( y_top+y_bot ) , &
3202  0.5_wp * ( z_top + z_bot ) , delta_solid(1:n_tot)
3203 
3204  ELSE
3205 
3206  IF ( u_atm .LT. 1.0d+3 ) THEN
3207 
3208  delta_angle = 2.0_wp*pi_g/n_cloud
3209 
3210  ELSE
3211 
3212  delta_angle = pi_g / ( n_cloud - 1.0_wp )
3213 
3214  END IF
3215 
3216 
3217  DO i=1,n_cloud
3218 
3219  start_angle = atan2(sin_theta,cos_theta)
3220  angle_release = (i-1) * delta_angle - 0.5_wp*pi_g
3221 
3222  dx = 0.5* ( r_bot + r_top ) * cos(start_angle + angle_release)
3223  dy = 0.5* ( r_bot + r_top ) * sin(start_angle + angle_release)
3224 
3225  dz = 0.0_wp
3226 
3227  IF ( verbose_level .GE. 1 ) THEN
3228 
3229  WRITE(*,110) 0.5_wp * ( x_top + x_bot ) + dx , &
3230  0.5_wp * ( y_top + y_bot ) + dy , &
3231  0.5_wp * ( z_top + z_bot ) + dz , &
3232  delta_solid(1:n_tot)/n_cloud
3233 
3234  END IF
3235 
3236  WRITE(hy_unit,110) 0.5_wp * ( x_top + x_bot ) + dx , &
3237  0.5_wp * ( y_top + y_bot ) + dy , &
3238  0.5_wp * ( z_top + z_bot ) + dz , &
3239  delta_solid(1:n_tot)/n_cloud
3240 
3241  END DO
3242 
3243  END IF
3244 
3245  ! WRITE THE RELEASE AT THE TOP OF THE COLUMN (OR NBL.)
3246 
3247  IF ( umbrella_flag ) THEN
3248 
3249  IF ( verbose_level .GE. 1 ) THEN
3250 
3251  WRITE(*,110) x_top+dx , y_top+dy , z_top+dz , &
3252  cloud_solid(1:n_tot)
3253 
3254  END IF
3255 
3256  WRITE(hy_unit,110) x_top+dx , y_top+dy , z_top+dz , &
3257  cloud_solid(1:n_tot)
3258 
3259  ELSE
3260 
3261  IF ( n_cloud .EQ. 1 ) THEN
3262 
3263  IF ( verbose_level .GE. 1 ) THEN
3264 
3265  WRITE(*,110) x_top , y_top , z_top , cloud_solid(1:n_tot)
3266 
3267  END IF
3268 
3269  WRITE(hy_unit,110) x_top , y_top , z_top , cloud_solid(1:n_tot)
3270  ! Write data for umbrella cloud initialization
3271  !OPEN(112, file = 'NBL.hy', status = 'replace')
3272  !WRITE(112,*) x_top , y_top , z_top , cloud_solid(1:n_tot)
3273  !CLOSE(112)
3274 
3275  ELSE
3276 
3277  IF ( u_atm .LT. 1.0d+3 ) THEN
3278 
3279  delta_angle = 2.0_wp*pi_g/n_cloud
3280 
3281  ELSE
3282 
3283  delta_angle = pi_g / ( n_cloud - 1.0_wp )
3284 
3285  END IF
3286 
3287  DO i=1,n_cloud
3288 
3289  start_angle = atan2(sin_theta,cos_theta)
3290  angle_release = (i-1) * delta_angle - 0.5_wp*pi_g
3291 
3292  dx = 0.5* ( r_bot + r_top ) * cos(start_angle + angle_release)
3293  dy = 0.5* ( r_bot + r_top ) * sin(start_angle + angle_release)
3294  dz = 0.0_wp
3295 
3296  IF ( verbose_level .GE. 1 ) THEN
3297 
3298  WRITE(*,110) x_top+dx , y_top+dy , z_top+dz , &
3299  cloud_solid(1:n_tot)/n_cloud
3300 
3301  END IF
3302 
3303  WRITE(hy_unit,110) x_top+dx , y_top+dy , z_top+dz , &
3304  cloud_solid(1:n_tot)/n_cloud
3305 
3306  END DO
3307 
3308  END IF
3309 
3310  END IF
3311 
3312  ! WRITE(*,*) 'z_max',z_max
3313  WRITE(*,*) 'Solid mass released in the atmosphere (kg/s): ',sum(solid_tot)
3314 
3315 107 FORMAT(1x,' x (m) ',1x,' y (m) ', 1x,' z (m) ')
3316 
3317 108 FORMAT(2x,a)
3318 
3319 110 FORMAT(90(1x,e15.8))
3320 
3321 
3322  ! Write hysplit file for volcanig gas only
3323 
3324  OPEN(hy_unit_volcgas,file=hy_file_volcgas)
3325 
3326  WRITE(hy_unit_volcgas,207,advance="no")
3327 
3328  DO i=1,n_gas
3329 
3330  WRITE(x1,'(I2.2)') i ! converting integer to string using a 'internal file'
3331 
3332  WRITE(hy_unit_volcgas,208,advance="no") 'VG fr '//trim(x1)//' (kg/s)'
3333 
3334  END DO
3335 
3336 
3337  WRITE(hy_unit_volcgas,*) ''
3338 
3339  z_min = z_col(1)
3340 
3341  IF ( nbl_stop ) THEN
3342 
3343  z_max = height_nbl + z_min
3344 
3345  ELSE
3346 
3347  z_max = z_col(col_lines)
3348 
3349  END IF
3350 
3351 
3352  ! WRITE(*,*) 'z_min',z_min
3353 
3354  n_hy = floor( ( z_max - z_min ) / hy_deltaz )
3355 
3356  z_bot = z_min + n_hy * hy_deltaz
3357  z_top = z_max
3358 
3359  !WRITE(*,*) 'volcgas_mass_flux : ',volcgas_mass_flux(n_gas,:)
3360  DO j = 1,n_gas
3361 
3362 
3363  CALL interp_1d_scalar(z_col, volcgas_mass_flux(j,:), z_top, gas_top)
3364 
3365  CALL interp_1d_scalar(z_col, x_col, z_bot, x_bot)
3366  CALL interp_1d_scalar(z_col, x_col, z_top, x_top)
3367 
3368  CALL interp_1d_scalar(z_col, x_col, z_bot, y_bot)
3369  CALL interp_1d_scalar(z_col, x_col, z_top, y_top)
3370 
3371  CALL interp_1d_scalar(z_col, r_col, z_bot, r_bot)
3372  CALL interp_1d_scalar(z_col, r_col, z_top, r_top)
3373 
3374 
3375  cloud_gas(j) = gas_top
3376 
3377  END DO
3378 
3379  !WRITE(*,*) 'cloud_gas(j) : ',gas_top
3380  !WRITE(*,*) 'cloud_gas(1:n_gas) : ',cloud_gas(1:n_gas)
3381 
3382  IF ( umbrella_flag ) THEN
3383 
3384  IF ( verbose_level .GE. 1 ) THEN
3385 
3386  WRITE(*,210) x_top+dx , y_top+dy , z_top+dz , &
3387  cloud_gas(1:n_gas)
3388 
3389  END IF
3390 
3391  WRITE(hy_unit_volcgas,210) x_top+dx , y_top+dy , z_top+dz , &
3392  cloud_gas(1:n_gas)
3393  ELSE
3394 
3395  IF ( n_cloud .EQ. 1 ) THEN
3396 
3397  IF ( verbose_level .GE. 1 ) THEN
3398 
3399  WRITE(*,210) x_top , y_top , z_top , cloud_gas(1:n_gas)
3400 
3401  END IF
3402 
3403  WRITE(hy_unit_volcgas,210) x_top , y_top , z_top , cloud_gas(1:n_gas)
3404 
3405  ELSE
3406 
3407  IF ( u_atm .LT. 1.0d+3 ) THEN
3408 
3409  delta_angle = 2.0_wp*pi_g/n_cloud
3410 
3411  ELSE
3412 
3413  delta_angle = pi_g / ( n_cloud - 1.0_wp )
3414 
3415  END IF
3416 
3417  DO i=1,n_cloud
3418 
3419  start_angle = atan2(sin_theta,cos_theta)
3420  angle_release = (i-1) * delta_angle - 0.5_wp*pi_g
3421 
3422  dx = 0.5* ( r_bot + r_top ) * cos(start_angle + angle_release)
3423  dy = 0.5* ( r_bot + r_top ) * sin(start_angle + angle_release)
3424 
3425 
3426  IF ( verbose_level .GE. 1 ) THEN
3427 
3428  WRITE(*,210) x_top+dx , y_top+dy , z_top , cloud_gas(1:n_gas) &
3429  / n_cloud
3430 
3431  END IF
3432 
3433  WRITE(hy_unit_volcgas,210) x_top+dx , y_top+dy , z_top , &
3434  cloud_gas(1:n_gas)/n_cloud
3435 
3436  END DO
3437 
3438  END IF
3439 
3440  END IF
3441 
3442 207 FORMAT(1x,' x (m) ',1x,' y (m) ', 1x,' z (m) ')
3443 
3444 208 FORMAT(2x,a)
3445 
3446 210 FORMAT(33(1x,e15.8))
3447 
3448  RETURN
3449 
3450  END SUBROUTINE check_hysplit
3451 
3452  !------------------------------------------------------------------------------
3456  !
3461  !------------------------------------------------------------------------------
3462 
3463  FUNCTION cross(a, b)
3464  REAL(wp), DIMENSION(3) :: cross
3465  REAL(wp), DIMENSION(3), INTENT(IN) :: a, b
3466 
3467  cross(1) = a(2) * b(3) - a(3) * b(2)
3468  cross(2) = a(3) * b(1) - a(1) * b(3)
3469  cross(3) = a(1) * b(2) - a(2) * b(1)
3470 
3471  END FUNCTION cross
3472 
3473  !------------------------------------------------------------------------------
3477  !
3484  !------------------------------------------------------------------------------
3485 
3486  FUNCTION normpdf(phi, mu,sigma )
3487  REAL(wp) :: normpdf
3488  REAL(wp), INTENT(IN) :: phi,mu,sigma
3489 
3490  normpdf = 1.0_wp / ( sigma * sqrt( 2.0_wp * pi_g ) ) * exp( -0.5_wp * &
3491  ( ( phi - mu ) / sigma )**2 )
3492 
3493  END FUNCTION normpdf
3494 
3495  !------------------------------------------------------------------------------
3499  !
3503  !------------------------------------------------------------------------------
3504 
3505  FUNCTION lower( string ) result (new)
3506  character(len=*) :: string
3507 
3508  character(len=len(string)) :: new
3509 
3510  integer :: i
3511  integer :: k
3512  INTEGER :: length
3513 
3514  length = len(string)
3515  new = string
3516  do i = 1,len(string)
3517  k = iachar(string(i:i))
3518  if ( k >= iachar('A') .and. k <= iachar('Z') ) then
3519  k = k + iachar('a') - iachar('A')
3520  new(i:i) = achar(k)
3521  endif
3522  enddo
3523  end function lower
3524 
3525  !------------------------------------------------------------------------------
3529  !
3540  !------------------------------------------------------------------------------
3541 
3542  FUNCTION numberfrommass(MassL,MassR,Mass) result (Number)
3543 
3544  REAL(wp), INTENT(IN) :: MassL(:)
3545  REAL(wp), INTENT(IN) :: MassR(:)
3546  REAL(wp), INTENT(IN) :: Mass(:)
3547 
3548  REAL(wp) :: Number(size(mass))
3549 
3550  REAL(wp) :: a , b
3551  REAL(wp) :: x1 , x2
3552 
3553  INTEGER :: j
3554 
3555  b = 0.0_wp
3556 
3557  DO j=n_sections,1,-1
3558 
3559  x1 = massl(j)
3560  x2 = massr(j)
3561  a = ( 6.0_wp * mass(j) / ( x2-x1 ) - b * ( x1+2.0_wp*x2 ) ) / ( 2.0_wp*x1+x2 )
3562 
3563  IF ( a .GT. 0.0_wp ) THEN
3564 
3565  ELSE
3566 
3567  a = 2.0_wp * mass(j) / ( x2**2-x1**2 )
3568  b = a
3569 
3570  END IF
3571 
3572  number(j) = 0.5_wp*(a+b)*(x2-x1)
3573 
3574  b = a
3575 
3576  END DO
3577 
3578  IF ( verbose_level .GE. 1 ) THEN
3579 
3580  WRITE(*,*) 'MassL'
3581  WRITE(*,"(30ES8.1)") massl
3582  WRITE(*,*) 'MassR'
3583  WRITE(*,"(30ES8.1)") massr
3584  WRITE(*,*) 'Mass'
3585  WRITE(*,"(30ES8.1)") mass
3586  WRITE(*,*) 'Number'
3587  WRITE(*,"(30ES8.1)") number
3588 
3589 
3590  READ(*,*)
3591 
3592  END IF
3593 
3594  RETURN
3595 
3596  END FUNCTION numberfrommass
3597 
3598  !------------------------------------------------------------------------------
3602  !
3608  !------------------------------------------------------------------------------
3609 
3610  CHARACTER*2 FUNCTION lettera(k)
3611  IMPLICIT NONE
3612  CHARACTER ones,tens,hund,thou
3613  !
3614  INTEGER :: k
3615  !
3616  INTEGER :: iten, ione, ihund, ithou
3617  !
3618  ithou=int(k/1000)
3619  ihund=int((k-(ithou*1000))/100)
3620  iten=int((k-(ithou*1000)-(ihund*100))/10)
3621  ione=k-ithou*1000-ihund*100-iten*10
3622  ones=char(ione+48)
3623  tens=char(iten+48)
3624  hund=char(ihund+48)
3625  thou=char(ithou+48)
3626  lettera=tens//ones
3627  !
3628  RETURN
3629  END FUNCTION lettera
3630 
3631 
3632 
3633 END MODULE inpout
inpout::write_zero_hysplit
subroutine write_zero_hysplit
Hysplit output initialization.
Definition: inpout.f90:2843
particles_module::phifromm
subroutine phifromm
Compute size from mass.
Definition: particles.f90:1291
inpout::close_file_units
subroutine close_file_units
Close output units.
Definition: inpout.f90:2581
particles_module::phi1
real(wp), dimension(:), allocatable phi1
First diameter for the density function.
Definition: particles.f90:74
inpout::hy_y
real(wp) hy_y
Definition: inpout.f90:143
moments_module::n_mom
integer n_mom
Number of moments.
Definition: moments.f90:24
inpout::col_file
character(len=30) col_file
Name of output file for values along the profile.
Definition: inpout.f90:75
particles_module::particles_shape
real(wp) function particles_shape(i_part, phi)
Particle shape factor.
Definition: particles.f90:662
variables::dakota_flag
logical dakota_flag
Flag for dakota run (less files on output)
Definition: variables.f90:36
particles_module::n_part
integer n_part
number of particle phases
Definition: particles.f90:23
inpout::n_unit
integer n_unit
Counter for unit files.
Definition: inpout.f90:60
parameters_2d::n_rk
integer n_rk
Runge-Kutta order.
Definition: parameters_2d.f90:55
inpout::delta_phi
real(wp) delta_phi
Definition: inpout.f90:141
inpout::write_dakota
subroutine write_dakota(description, value)
Dakota outputs.
Definition: inpout.f90:2779
inpout::bak_file
character(len=30) bak_file
Name of output file for backup of input parameters.
Definition: inpout.f90:66
meteo_module::h2
real(wp) h2
Definition: meteo.f90:20
meteo_module::pa
real(wp) pa
Atmospheric pressure.
Definition: meteo.f90:79
inpout::initialize
subroutine initialize
Initialize variables.
Definition: inpout.f90:211
parameters_2d::rsource_cells
integer rsource_cells
Definition: parameters_2d.f90:34
inpout::lower
character(len=len(string)) function lower(string)
Convert from capital to small.
Definition: inpout.f90:3506
meteo_module::rho_atm
real(wp) rho_atm
Atmospheric density.
Definition: meteo.f90:60
variables::inversion_flag
logical inversion_flag
Definition: variables.f90:67
mixture_module::ice_mass_fraction
real(wp) ice_mass_fraction
mass fraction of ice in the mixture
Definition: mixture.f90:118
particles_module::distribution
character(len=20) distribution
Ditribution of the particles: .
Definition: particles.f90:110
mixture_module::liquid_water_mass_fraction
real(wp) liquid_water_mass_fraction
mass fraction of liquid water in the mixture
Definition: mixture.f90:106
mixture_module::volcgas_mass_fraction
real(wp), dimension(:), allocatable volcgas_mass_fraction
mass fractions of volcanic gases
Definition: mixture.f90:82
inpout::hy_deltaz
real(wp) hy_deltaz
Definition: inpout.f90:143
parameters_2d::dt0
real(wp) dt0
Initial time step.
Definition: parameters_2d.f90:18
meteo_module::cpair
real(wp) cpair
specific heat capacity for dry air
Definition: meteo.f90:91
plume_module::w0
real(wp) w0
initial vertical velocity of the plume
Definition: plume.f90:36
meteo_module::rh
real(wp) rh
Relative humidity for standard atmosphere.
Definition: meteo.f90:45
particles_module::m_quad
real(wp), dimension(:,:,:), allocatable m_quad
Abscissa of quadrature formulas.
Definition: particles.f90:131
particles_module::aggregation_model
character(len=20) aggregation_model
Aggregation kernel model: .
Definition: particles.f90:125
parameters_2d
Parameters.
Definition: parameters_2d.f90:7
meteo_module::u_max
real(wp) u_max
Definition: meteo.f90:136
particles_module::shape1
real(wp), dimension(:), allocatable shape1
Shape factor at phi=phi1.
Definition: particles.f90:80
meteo_module::z_r
real(wp) z_r
Definition: meteo.f90:136
particles_module::cum_particle_loss_rate
real(wp), dimension(:,:), allocatable cum_particle_loss_rate
cumulative rate of particles lost up to the integration height ( kg s^-1)
Definition: particles.f90:50
particles_module::init_quadrature_points
subroutine init_quadrature_points
Quadrature initialization.
Definition: particles.f90:1227
inpout::atm_unit
integer atm_unit
Atmosphere input unit.
Definition: inpout.f90:96
inpout::hy_file_volcgas
character(len=30) hy_file_volcgas
Name of output file for hysplit volcanic gas.
Definition: inpout.f90:81
solver_module::dz0
real(wp) dz0
Initial integration step.
Definition: solver_rise.f90:52
inpout::atm_file
character(len=50) atm_file
Name of file for the parameters of the atmosphere.
Definition: inpout.f90:93
mixture_module::volcgas_mol_wt
real(wp), dimension(:), allocatable volcgas_mol_wt
molecular weight of additional volcanic gases
Definition: mixture.f90:76
inpout::mom_file
character(len=30) mom_file
Name of output file for backup of input parameters.
Definition: inpout.f90:84
mixture_module::gas_volume_fraction
real(wp) gas_volume_fraction
gas vlume fraction in the mixture
Definition: mixture.f90:22
mixture_module::water_vapor_mass_fraction
real(wp) water_vapor_mass_fraction
mass fraction of water vapor in the mixture
Definition: mixture.f90:112
variables::n_cloud
integer n_cloud
Definition: variables.f90:58
particles_module::allocate_particles
subroutine allocate_particles
Particles variables allocation.
Definition: particles.f90:198
moments_module::n_sections
integer n_sections
Definition: moments.f90:26
mixture_module::rho_ice
real(wp) rho_ice
Density of ice in the mixture.
Definition: mixture.f90:133
variables::water_flag
logical water_flag
Flag for water.
Definition: variables.f90:70
inpout::hy_z_old
real(wp) hy_z_old
Definition: inpout.f90:143
variables::pi_g
real *8 pi_g
Greek pi.
Definition: variables.f90:30
inpout::inp_file
character(len=30) inp_file
Name of input file.
Definition: inpout.f90:63
mixture_module::water_volume_fraction0
real(wp) water_volume_fraction0
initial water volume fraction
Definition: mixture.f90:52
meteo_module::ta
real(wp) ta
Atmospheric temperature.
Definition: meteo.f90:76
mixture_module::dry_air_mass_fraction
real(wp) dry_air_mass_fraction
mass fraction of dry air in the mixture
Definition: mixture.f90:100
moments_module::n_nodes
integer n_nodes
Number of nodes for the quadrature.
Definition: moments.f90:18
variables::write_flag
logical write_flag
Definition: variables.f90:74
particles_module
Particles module.
Definition: particles.f90:11
plume_module::x
real(wp) x
plume location (downwind)
Definition: plume.f90:19
inpout::sed_file
character(len=30) sed_file
Name of output file for particle loss rate.
Definition: inpout.f90:72
mixture_module::rvolcgas_mix
real(wp) rvolcgas_mix
gas constant of volcanic gases mixture ( J/(kg K) )
Definition: mixture.f90:88
particles_module::deallocate_particles
subroutine deallocate_particles
Definition: particles.f90:282
particles_module::rho1
real(wp), dimension(:), allocatable rho1
Density at phi=phi1.
Definition: particles.f90:77
mixture_module::volcgas_mix_mass_fraction
real(wp) volcgas_mix_mass_fraction
mass fraction of the volcanic gas in the mixture
Definition: mixture.f90:97
meteo_module::u_atm
real(wp) u_atm
Horizonal wind speed.
Definition: meteo.f90:54
meteo_module::wind_mult_coeff
real(wp) wind_mult_coeff
Definition: meteo.f90:143
variables::aggregation_flag
logical aggregation_flag
Definition: variables.f90:72
inpout
Input/Output module.
Definition: inpout.f90:11
parameters_2d::n_vars
integer n_vars
Number of conservative variables.
Definition: parameters_2d.f90:50
inpout::solid_mfr
real(wp), dimension(:), allocatable solid_mfr
Definition: inpout.f90:145
inpout::sed_unit
integer sed_unit
Particle loss values along the column data unit.
Definition: inpout.f90:109
inpout::dak_unit
integer dak_unit
Dakota variables data unit.
Definition: inpout.f90:129
inpout::hy_unit_volcgas
integer hy_unit_volcgas
hysplit volcanic gas data unit
Definition: inpout.f90:120
mixture_module::added_water_temp
real(wp) added_water_temp
Definition: mixture.f90:135
particles_module::f_quad
real(wp), dimension(:,:,:), allocatable f_quad
Values of linear reconstructions at quadrature points.
Definition: particles.f90:158
plume_module::beta_inp
real(wp) beta_inp
entrainment coefficient (normal direction)
Definition: plume.f90:30
variables::hysplit_flag
logical hysplit_flag
Flag for hysplit run.
Definition: variables.f90:39
variables::nbl_stop
logical nbl_stop
Flag for hysplit output .
Definition: variables.f90:45
variables::r_min
real(wp) r_min
Definition: variables.f90:77
inpout::hy_y_old
real(wp) hy_y_old
Definition: inpout.f90:143
meteo_module::exp_wind
real(wp) exp_wind
Definition: meteo.f90:136
parameters_2d::dt_output
real(wp) dt_output
time interval for the output of the solution
Definition: parameters_2d.f90:81
mixture_module::added_water_mass_fraction
real(wp) added_water_mass_fraction
Definition: mixture.f90:137
plume_module::alpha_inp
real(wp) alpha_inp
entrainment coefficient (parallel direction)
Definition: plume.f90:29
particles_module::rho_quad
real(wp), dimension(:,:,:), allocatable rho_quad
Particle densities at quadrature points.
Definition: particles.f90:146
inpout::write_inversion
subroutine write_inversion(r0, w_opt, opt_mfr, opt_height, search_flag, opt_regime)
Write inversion file.
Definition: inpout.f90:2819
mixture_module::cpvolcgas_mix
real(wp) cpvolcgas_mix
specific heat of volcanic gases mixture
Definition: mixture.f90:91
particles_module::cp_part
real(wp), dimension(:), allocatable cp_part
Heat capacity of particle phases.
Definition: particles.f90:92
inpout::inversion_file
character(len=30) inversion_file
Name of output file for the inversion variables.
Definition: inpout.f90:90
plume_module::y
real(wp) y
plume location (crosswind)
Definition: plume.f90:20
meteo_module::sphu_atm0
real(wp) sphu_atm0
Atmospheric specific humidity at sea level (kg/kg)
Definition: meteo.f90:66
inpout::p_atm0
real(wp) p_atm0
Definition: inpout.f90:148
particles_module::bin_partial_mass_fraction
real(wp), dimension(:,:), allocatable bin_partial_mass_fraction
mass fraction of the bins of particle with respect to the total solid
Definition: particles.f90:44
particles_module::solid_partial_mass_fraction
real(wp), dimension(:), allocatable solid_partial_mass_fraction
mass fraction of the particle phases with respect to the total solid
Definition: particles.f90:26
inpout::notset
real(wp) notset
Definition: inpout.f90:56
parameters_2d::c_d
real(wp) c_d
Definition: parameters_2d.f90:36
inpout::sigma_lognormal
real(wp), dimension(:), allocatable sigma_lognormal
Definition: inpout.f90:136
parameters_2d::cfl
real(wp) cfl
Courant-Friedrichs-Lewy parameter.
Definition: parameters_2d.f90:22
meteo_module
Meteo module.
Definition: meteo.f90:11
plume_module::r0
real(wp) r0
initial radius of the plume
Definition: plume.f90:37
mixture_module::rgasmix
real(wp) rgasmix
universal constant for the mixture
Definition: mixture.f90:28
meteo_module::visc_atm0
real(wp) visc_atm0
Atmospheric kinematic viscosity at sea level.
Definition: meteo.f90:73
inpout::solid_mfr_init
real(wp), dimension(:), allocatable solid_mfr_init
Definition: inpout.f90:145
variables::w_min
real(wp) w_min
Definition: variables.f90:79
meteo_module::interp_1d_scalar
subroutine interp_1d_scalar(x1, f1, x2, f2)
Scalar interpolation.
Definition: meteo.f90:511
mixture_module::atm_mass_fraction
real(wp) atm_mass_fraction
mass fraction of the entrained air in the mixture
Definition: mixture.f90:94
inpout::normpdf
real(wp) function normpdf(phi, mu, sigma)
Normal probability density function.
Definition: inpout.f90:3487
variables::steady_flag
logical steady_flag
Flag to stop the umbrella solver when a steady upwind spreading is reached.
Definition: variables.f90:56
inpout::inp_unit
integer inp_unit
Input data unit.
Definition: inpout.f90:103
meteo_module::pres_atm_month_lat
real(wp), dimension(:), allocatable pres_atm_month_lat
Definition: meteo.f90:138
parameters_2d::solver_scheme
character(len=20) solver_scheme
Finite volume method: .
Definition: parameters_2d.f90:75
mixture_module
Gas/particles mixture module.
Definition: mixture.f90:11
mixture_module::rho_mix
real(wp) rho_mix
mixture density
Definition: mixture.f90:31
particles_module::m
real(wp), dimension(:,:), allocatable m
boundaries of the sections in mass scale (kg)
Definition: particles.f90:170
inpout::hy_file
character(len=30) hy_file
Name of output file for hysplit.
Definition: inpout.f90:78
plume_module
Plume module.
Definition: plume.f90:11
inpout::t_atm0
real(wp) t_atm0
Definition: inpout.f90:148
inpout::mom_unit
integer mom_unit
Moments values along the column data unit.
Definition: inpout.f90:126
inpout::open_file_units
subroutine open_file_units
Initialize output units.
Definition: inpout.f90:2504
meteo_module::cos_theta
real(wp) cos_theta
Wind angle.
Definition: meteo.f90:48
inpout::write_column
subroutine write_column
Write outputs.
Definition: inpout.f90:2611
inpout::mfr0
real(wp) mfr0
Definition: inpout.f90:134
variables::r_max
real(wp) r_max
Definition: variables.f90:78
meteo_module::h_levels
real(wp), dimension(:), allocatable h_levels
Definition: meteo.f90:141
inpout::col_lines
integer col_lines
Definition: inpout.f90:116
meteo_module::zmet
subroutine zmet
Meteo parameters.
Definition: meteo.f90:229
inpout::hy_z
real(wp) hy_z
Definition: inpout.f90:143
inpout::solid_mfr_oldold
real(wp), dimension(:), allocatable solid_mfr_oldold
Definition: inpout.f90:145
meteo_module::rair
real(wp) rair
perfect gas constant for dry air ( J/(kg K) )
Definition: meteo.f90:88
moments_module
Method of Moments module.
Definition: moments.f90:11
inpout::phi_min
real(wp) phi_min
Definition: inpout.f90:141
plume_module::particles_loss
logical particles_loss
logical defining if we loose particles
Definition: plume.f90:32
mixture_module::rvolcgas
real(wp), dimension(:), allocatable rvolcgas
gas constants for volcanic gases
Definition: mixture.f90:70
variables::isset
logical function isset(var)
Input variable check.
Definition: variables.f90:104
variables::height_nbl
real(wp) height_nbl
Definition: variables.f90:60
particles_module::mom
real(wp), dimension(:,:,:), allocatable mom
Moments of the particles diameter.
Definition: particles.f90:53
particles_module::phil
real(wp), dimension(:), allocatable phil
left boundaries of the sections in phi-scale
Definition: particles.f90:164
parameters_2d::t_end
real(wp) t_end
end time for the run
Definition: parameters_2d.f90:79
inpout::check_hysplit
subroutine check_hysplit
Hysplit outputs.
Definition: inpout.f90:2908
inpout::hy_x
real(wp) hy_x
Definition: inpout.f90:143
inpout::dak_file
character(len=30) dak_file
Name of output file for the variables used by dakota.
Definition: inpout.f90:87
variables::entr_abv_nbl_flag
logical entr_abv_nbl_flag
Flag to allow for entrainment above neutral buoyancy level.
Definition: variables.f90:50
plume_module::mag_u
real(wp) mag_u
velocity magnitude along the centerline
Definition: plume.f90:26
particles_module::solid_partial_mass_fraction0
real(wp), dimension(:), allocatable solid_partial_mass_fraction0
init mass fraction of the particle phases with respect to the total solid
Definition: particles.f90:29
inpout::read_inp
subroutine read_inp
Read Input data.
Definition: inpout.f90:428
meteo_module::atm_profile
real(wp), dimension(:,:), allocatable atm_profile
atmospheric profile above the vent. It is an array with n_atm_profile rows and 7 columns: ...
Definition: meteo.f90:132
particles_module::w_quad
real(wp), dimension(:,:,:), allocatable w_quad
Weights of quadrature formulas.
Definition: particles.f90:134
solver_module
Solver module.
Definition: solver_rise.f90:11
meteo_module::rho_atm_month_lat
real(wp), dimension(:), allocatable rho_atm_month_lat
Definition: meteo.f90:138
mixture_module::initial_neutral_density
logical initial_neutral_density
logical defining if the plume has neutral density at the base
Definition: mixture.f90:34
particles_module::shape_factor
real(wp), dimension(:), allocatable shape_factor
shape factor for settling velocity (Pfeiffer)
Definition: particles.f90:71
plume_module::vent_height
real(wp) vent_height
height of the base of the plume
Definition: plume.f90:35
plume_module::z
real(wp) z
plume vertical coordinate
Definition: plume.f90:21
variables::umbrella_flag
logical umbrella_flag
Flag to solve the model for the umbrella spreading.
Definition: variables.f90:53
inpout::hy_lines
integer hy_lines
Definition: inpout.f90:114
particles_module::particles_beta0
real(wp) particles_beta0
Definition: particles.f90:98
parameters_2d::limiter
integer, dimension(10) limiter
Limiter for the slope in the linear reconstruction: .
Definition: parameters_2d.f90:68
mixture_module::cpvolcgas
real(wp), dimension(:), allocatable cpvolcgas
specific heat capacity for volcanic gases
Definition: mixture.f90:73
variables::wp
integer, parameter wp
working precision
Definition: variables.f90:21
variables::n_values
integer n_values
Definition: variables.f90:81
particles_module::solid_mass_fraction0
real(wp), dimension(:), allocatable solid_mass_fraction0
initial mass fraction of the particle phases with respect to the mixture
Definition: particles.f90:38
inpout::cross
real(wp) function, dimension(3) cross(a, b)
Cross product.
Definition: inpout.f90:3464
inpout::numberfrommass
real(wp) function, dimension(size(mass)) numberfrommass(MassL, MassR, Mass)
Particle number from mass in each bin.
Definition: inpout.f90:3543
particles_module::phir
real(wp), dimension(:), allocatable phir
right boundaries of the sections in phi-scale
Definition: particles.f90:167
mixture_module::t_mix
real(wp) t_mix
mixture temperature
Definition: mixture.f90:37
parameters_2d::max_dt
real(wp) max_dt
Largest time step allowed.
Definition: parameters_2d.f90:20
meteo_module::temp_atm_month_lat
real(wp), dimension(:), allocatable temp_atm_month_lat
Definition: meteo.f90:138
meteo_module::initialize_meteo
subroutine initialize_meteo
Meteo parameters initialization.
Definition: meteo.f90:160
inpout::run_name
character(len=30) run_name
Name of the run (used for the output and backup files)
Definition: inpout.f90:69
mixture_module::rhovolcgas_mix
real(wp) rhovolcgas_mix
volcanic gases mixture density
Definition: mixture.f90:85
parameters_2d::t_start
real(wp) t_start
initial time for the run
Definition: parameters_2d.f90:78
particles_module::particles_density
real(wp) function particles_density(i_part, phi)
Particle density.
Definition: particles.f90:607
variables::gi
real *8 gi
Gravity acceleration.
Definition: variables.f90:24
inpout::mu_lognormal
real(wp), dimension(:), allocatable mu_lognormal
Definition: inpout.f90:136
inpout::lat
real(wp) lat
Definition: inpout.f90:139
mixture_module::volcgas_mass_fraction0
real(wp), dimension(:), allocatable volcgas_mass_fraction0
initial mass fractions of volcanic gases
Definition: mixture.f90:79
mixture_module::rhovolcgas
real(wp), dimension(:), allocatable rhovolcgas
volcanic gases densities
Definition: mixture.f90:67
plume_module::r
real(wp) r
plume radius
Definition: plume.f90:22
meteo_module::rwv
real(wp), parameter rwv
gas constant for water vapor ( J/(kg K) )
Definition: meteo.f90:118
meteo_module::rel_hu
real(wp) rel_hu
Definition: meteo.f90:36
mixture_module::t_mix0
real(wp) t_mix0
initial temperature
Definition: mixture.f90:49
inpout::month
real(wp) month
Definition: inpout.f90:138
inpout::lettera
character *2 function lettera(k)
Numeric to String conversion.
Definition: inpout.f90:3611
mixture_module::rho_lw
real(wp) rho_lw
Density of liquid water in the mixture.
Definition: mixture.f90:130
particles_module::phi2
real(wp), dimension(:), allocatable phi2
Second diameter for the density function.
Definition: particles.f90:83
parameters_2d::theta
real(wp) theta
Van Leer limiter parameter.
Definition: parameters_2d.f90:77
particles_module::eval_quad_values
subroutine eval_quad_values
Quadrature values computation.
Definition: particles.f90:1179
inpout::hy_x_old
real(wp) hy_x_old
Definition: inpout.f90:143
plume_module::log10_mfr
real(wp) log10_mfr
Definition: plume.f90:38
particles_module::mom0
real(wp), dimension(:,:,:), allocatable mom0
Initial moments of the particles diameter.
Definition: particles.f90:56
particles_module::rho2
real(wp), dimension(:), allocatable rho2
Density at phi=phi2.
Definition: particles.f90:86
inpout::col_unit
integer col_unit
Output values along the column data unit.
Definition: inpout.f90:106
variables::height_obj
real(wp) height_obj
Definition: variables.f90:76
mixture_module::water_mass_fraction0
real(wp) water_mass_fraction0
initial water mass fraction
Definition: mixture.f90:55
inpout::solid_mfr_old
real(wp), dimension(:), allocatable solid_mfr_old
Definition: inpout.f90:145
mixture_module::n_gas
integer n_gas
volcanic gas species number
Definition: mixture.f90:64
parameters_2d::reconstr_coeff
real(wp) reconstr_coeff
Slope coefficient in the linear reconstruction.
Definition: parameters_2d.f90:26
inpout::hy_unit
integer hy_unit
hysplit data unit
Definition: inpout.f90:112
meteo_module::n_atm_profile
integer n_atm_profile
Definition: meteo.f90:120
variables
Global variables.
Definition: variables.f90:10
particles_module::shape2
real(wp), dimension(:), allocatable shape2
Shape factor at phi=phi2.
Definition: particles.f90:89
mixture_module::gas_mass_fraction
real(wp) gas_mass_fraction
gas mass fraction in the mixture
Definition: mixture.f90:19
plume_module::w
real(wp) w
plume vertical velocity
Definition: plume.f90:25
particles_module::settling_model
character *10 settling_model
Settling model: .
Definition: particles.f90:104
parameters_2d::interfaces_relaxation
logical interfaces_relaxation
Flag to add the relaxation terms after the linear reconstruction: .
Definition: parameters_2d.f90:32
meteo_module::read_atm_profile
character *10 read_atm_profile
Definition: meteo.f90:134
inpout::inversion_unit
integer inversion_unit
Inversion variables data unit.
Definition: inpout.f90:132
meteo_module::h1
real(wp) h1
Definition: meteo.f90:19
meteo_module::sin_theta
real(wp) sin_theta
Definition: meteo.f90:48
inpout::bak_unit
integer bak_unit
Backup input unit.
Definition: inpout.f90:100
variables::verbose_level
integer verbose_level
Level of verbose output (0 = minimal output on screen)
Definition: variables.f90:33
inpout::temp_unit
integer temp_unit
hysplit scratch unit
Definition: inpout.f90:123
variables::w_max
real(wp) w_max
Definition: variables.f90:80
Generated on Tue Jun 23 2020 12:45:55 for PLUME-MoM-TSM by   doxygen 1.8.11