forked from firemodels/fds
-
Notifications
You must be signed in to change notification settings - Fork 0
/
part.f90
4184 lines (3395 loc) · 158 KB
/
part.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!> \brief Variables and subroutines pertaining to Lagrangian particles
MODULE PART
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE TRAN
USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE,PACK_PARTICLE
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
USE SOOT_ROUTINES, ONLY: DROPLET_SCRUBBING
IMPLICIT NONE (TYPE,EXTERNAL)
PRIVATE
PUBLIC INSERT_ALL_PARTICLES,MOVE_PARTICLES,PARTICLE_MASS_ENERGY_TRANSFER,REMOVE_PARTICLES,&
GENERATE_PARTICLE_DISTRIBUTIONS,PARTICLE_MOMENTUM_TRANSFER
CONTAINS
!> \brief Generate size distribution statistics for all particles classes
SUBROUTINE GENERATE_PARTICLE_DISTRIBUTIONS
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE PHYSICAL_FUNCTIONS, ONLY : PARTICLE_SIZE_DISTRIBUTION
REAL(EB) :: LL,UL,BIN_SIZE,TNOW,DD,DI
INTEGER :: I,J,IL,IU,ILPC
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC=>NULL()
TYPE (RAMPS_TYPE), POINTER :: RM=>NULL()
IF (N_LAGRANGIAN_CLASSES==0) RETURN ! Don't waste time if no particles
TNOW=CURRENT_TIME()
PART_CLASS_LOOP: DO ILPC=1,N_LAGRANGIAN_CLASSES
LPC=>LAGRANGIAN_PARTICLE_CLASS(ILPC)
IF_SIZE_DISTRIBUTION: IF (.NOT.LPC%MONODISPERSE .AND. (LPC%DIAMETER>0._EB) .OR. LPC%CNF_RAMP_INDEX>=0) THEN
IF (LPC%CNF_RAMP_INDEX<0) THEN
CALL PARTICLE_SIZE_DISTRIBUTION(LPC%DIAMETER,LPC%R_CNF(:),LPC%CNF(:),LPC%CVF(:),NDC,LPC%GAMMA,LPC%SIGMA,LPC%DISTRIBUTION)
ELSE
RM=>RAMPS(LPC%CNF_RAMP_INDEX)
LPC%DIAMETER = RM%SPAN*0.5_EB ! Not used for MONODISPERSE
DD = RM%SPAN/NDC
LPC%R_CNF(0) = RM%T_MIN
LPC%CNF(0) = 0._EB
LPC%CVF(0) = 0._EB
DO I=1,NDC
DI = RM%T_MIN + I*DD
LPC%R_CNF(I) = DI
LPC%CNF(I) = EVALUATE_RAMP(DI,LPC%CNF_RAMP_INDEX)
LPC%CVF(I) = LPC%CVF(I-1) + (DI-0.5*DD)**3*(LPC%CNF(I)-LPC%CNF(I-1))
ENDDO
LPC%R_CNF = 1.E-6_EB*0.5_EB*LPC%R_CNF ! Convert diameter in microns to radius in meters.
LPC%CNF = LPC%CNF/LPC%CNF(NDC)
LPC%CVF = LPC%CVF/LPC%CVF(NDC)
ENDIF
BIN_SIZE = (LPC%R_CNF(NDC)-LPC%R_CNF(0))/REAL(LPC%N_STRATA,EB)
STRATIFY: DO I=1,LPC%N_STRATA
LL = LPC%R_CNF(0) + (I-1)*BIN_SIZE
UL = LPC%R_CNF(0) + I *BIN_SIZE
IL = 1
LL_LOOP: DO J=1,NDC
IF (LPC%R_CNF(J)>LL) THEN
IL = J-1
LPC%STRATUM_INDEX_LOWER(I) = J-1
EXIT LL_LOOP
ENDIF
ENDDO LL_LOOP
IU = NDC
UL_LOOP: DO J=NDC,1,-1
IF (LPC%R_CNF(J)<=UL) THEN
IU = J
LPC%STRATUM_INDEX_UPPER(I) = J
EXIT UL_LOOP
ENDIF
ENDDO UL_LOOP
LPC%W_CNF(I) = LPC%CNF(IU) - LPC%CNF(IL)
ENDDO STRATIFY
! Compute mean droplet volume for distribution
LPC%MEAN_DROPLET_VOLUME = 0._EB
DO I=1,NDC
LPC%MEAN_DROPLET_VOLUME = LPC%MEAN_DROPLET_VOLUME + ( LPC%CNF(I) - LPC%CNF(I-1) ) * FOTHPI*LPC%R_CNF(I-1)**3
ENDDO
ENDIF IF_SIZE_DISTRIBUTION
! If pacticles/PARTICLEs can break up, compute normalized (median = 1) size distribution for child PARTICLEs
IF (LPC%BREAKUP .AND. .NOT.LPC%MONODISPERSE) THEN
IF(LPC%BREAKUP_CNF_RAMP_INDEX<0) THEN
CALL PARTICLE_SIZE_DISTRIBUTION(1._EB,LPC%BREAKUP_R_CNF(:),LPC%BREAKUP_CNF(:),LPC%BREAKUP_CVF(:),NDC, &
LPC%BREAKUP_GAMMA,LPC%BREAKUP_SIGMA,LPC%BREAKUP_DISTRIBUTION)
ELSE
RM=>RAMPS(LPC%BREAKUP_CNF_RAMP_INDEX)
DD=RM%SPAN/NDC
DO I=1,NDC
DI=RM%T_MIN+(I-0.5_EB)*DD
LPC%BREAKUP_R_CNF(I) = 1E-6_EB*0.5_EB*DI
LPC%BREAKUP_CNF(I) = EVALUATE_RAMP(DI,LPC%BREAKUP_CNF_RAMP_INDEX)
ENDDO
LPC%BREAKUP_CNF=LPC%BREAKUP_CNF/LPC%BREAKUP_CNF(NDC)
ENDIF
ENDIF
ENDDO PART_CLASS_LOOP
T_USED(8)=T_USED(8)+CURRENT_TIME()-TNOW
END SUBROUTINE GENERATE_PARTICLE_DISTRIBUTIONS
!> \brief Insert Lagrangian particles into the domain every time step
SUBROUTINE INSERT_ALL_PARTICLES(T,NM)
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,RANDOM_CHOICE
USE GEOMETRY_FUNCTIONS, ONLY: RANDOM_RECTANGLE,RANDOM_CONE,RANDOM_RING,CONE_MESH_INTERSECTION_VOLUME,UNIFORM_RING,&
RING_MESH_INTERSECTION_ARC
USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER, SHUTDOWN
USE TRAN, ONLY: GET_IJK
USE DEVICE_VARIABLES
USE CONTROL_VARIABLES
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
REAL :: RN,RN2,RN3,INIT_CALLS
REAL(EB) :: PHI_RN,FLOW_RATE,THETA_RN,SPHI,CPHI,MASS_SUM,D_PRES_FACTOR, &
STHETA,CTHETA,PWT0,PARTICLE_SPEED,SHIFT1,SHIFT2,XTMP,YTMP,ZTMP,VLEN, &
TRIGT1,TRIGT2,TNOW,TSI,PIPE_PRESSURE,X1,X2,Y1,Y2,Z1,Z2, &
ETA,ETA_MAX,ETA_MIN,XI,YJ,ZK,X_OFFSET,Y_OFFSET,Z_OFFSET, &
VXMIN,VXMAX,VYMIN,VYMAX,VZMIN,VZMAX,VDX,VDY,VDZ
REAL(EB), PARAMETER :: VENT_OFFSET=0.5
INTEGER :: IP,KS,II,JJ,KK,IC,IL,IU,ILPC,DROP_SUM,IIG,JJG,KKG,IW,IOR,STRATUM,INIT_INDEX,ICF,ITS,ITERATIONS
INTEGER :: N,N_INSERT,ILAT,NEW_LP_INDEX,III,IOS,LU_VEG_IN,NVOX,VXI
INTEGER, ALLOCATABLE, DIMENSION(:) :: LP_INDEX_LOOKUP
LOGICAL :: INSERT_ANOTHER_BATCH
TYPE (PROPERTY_TYPE), POINTER :: PY=>NULL()
TYPE (TABLES_TYPE), POINTER :: TA=>NULL()
TYPE (DEVICE_TYPE), POINTER :: DV=>NULL()
TYPE (SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP=>NULL()
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC=>NULL()
TYPE (INITIALIZATION_TYPE), POINTER :: IN=>NULL()
TYPE (WALL_TYPE), POINTER :: WC
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
IF (N_LAGRANGIAN_CLASSES==0) RETURN ! Don't waste time if no particles
TNOW=CURRENT_TIME()
CALL POINT_TO_MESH(NM)
! Insert particles at spray nozzles (INSERT_SPRAY_PARTICLES), VENT surfaces (INSERT_VENT_PARTICLES), specified
! volumes (INSERT_VOLUMETRIC_PARTICLES), or HVAC ducts (INSERT_DUCT_PARTICLES). If the insertions are to be made according
! to a user-specified clock, allow the possibility via the OVERALL_INSERT_LOOP of inserting multiple sets of particles.
OVERALL_INSERT_LOOP: DO
INSERT_ANOTHER_BATCH = .FALSE.
! Insert spray particles/droplets
CALL INSERT_SPRAY_PARTICLES
! Insert particles at a solid boundary or vent.
CALL INSERT_VENT_PARTICLES
! Insert particles within a volume specified by an INIT line.
DO INIT_INDEX=1,N_INIT
IN => INITIALIZATION(INIT_INDEX)
IF (IN%INVOKED_BY_SURF) CYCLE
X_OFFSET = 0._EB
Y_OFFSET = 0._EB
Z_OFFSET = 0._EB
! Check if INIT line links to a binary file of bulk density data
IF (IN%BULK_DENSITY_FILE=='null') THEN
CALL INSERT_VOLUMETRIC_PARTICLES
ELSE
IF (.NOT. IN%ALREADY_INSERTED(NM)) THEN
! Build out of blocks
IN%SHAPE = 'BLOCK'
LU_VEG_IN = GET_FILE_NUMBER()
OPEN(UNIT=LU_VEG_IN,FILE=TRIM(IN%BULK_DENSITY_FILE),&
STATUS='OLD',FORM='UNFORMATTED',ACTION='READ',IOSTAT=IOS)
IF (IOS==0) THEN
READ(LU_VEG_IN) VXMIN,VXMAX,VYMIN,VYMAX,VZMIN,VZMAX
! Skip if volume containing vegetation is entirely outside the current mesh
IF (VXMIN>XF .OR. VXMAX<XS .OR. VYMIN>YF .OR. VYMAX<YS .OR. VZMIN>ZF .OR. VZMAX<ZS) THEN
IN%ALREADY_INSERTED(NM)=.TRUE.
CYCLE
ENDIF
! Voxel resolution
READ(LU_VEG_IN) VDX,VDY,VDZ
! Number of vegetation containing voxels
READ(LU_VEG_IN) NVOX
! Create a pseudo-INIT block for each voxel
VEG_INSERT_LOOP: DO VXI=1,NVOX
CALL INSERT_VOLUMETRIC_PARTICLES
ENDDO VEG_INSERT_LOOP
CLOSE(LU_VEG_IN)
IN%ALREADY_INSERTED(NM) = .TRUE.
ELSE
WRITE(MESSAGE,'(A,I0,A,A,A)') 'ERROR: INIT ',INIT_INDEX,', could not read binary bulk density file ', &
TRIM(IN%BULK_DENSITY_FILE),'. Check file exists.'
CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.); RETURN
ENDIF
ENDIF
ENDIF
ENDDO
! Insert volumetric particles that are invoked by a particular SURF type.
IF (INIT_INVOKED_BY_SURF) THEN
DO IW=1,N_INTERNAL_WALL_CELLS+N_EXTERNAL_WALL_CELLS
WC => MESHES(NM)%WALL(IW)
BC => MESHES(NM)%BOUNDARY_COORD(WC%BC_INDEX)
SF => SURFACE(WC%SURF_INDEX)
IF (BC%IOR==3 .AND. SF%INIT_INDICES(1)>0) THEN
CALL RANDOM_NUMBER(RN)
X_OFFSET = BC%X + (RN-0.5_FB)*DX(BC%IIG)
CALL RANDOM_NUMBER(RN)
Y_OFFSET = BC%Y + (RN-0.5_FB)*DY(BC%JJG)
Z_OFFSET = BC%Z
B1 => MESHES(NM)%BOUNDARY_PROP1(WC%B1_INDEX)
INIT_CALLS = REAL(SF%INIT_PER_AREA*B1%AREA)
ITERATIONS = INT(INIT_CALLS)
CALL RANDOM_NUMBER(RN)
IF (RN<INIT_CALLS-ITERATIONS) ITERATIONS = ITERATIONS + 1
DO ITS=1,ITERATIONS
DO III=1,10
IF (SF%INIT_IDS(III)=='null') CYCLE
INIT_INDEX = SF%INIT_INDICES(III)
CALL INSERT_VOLUMETRIC_PARTICLES
ENDDO
ENDDO
ENDIF
ENDDO
DO INIT_INDEX=1,N_INIT
IF (INITIALIZATION(INIT_INDEX)%INVOKED_BY_SURF) INITIALIZATION(INIT_INDEX)%ALREADY_INSERTED = .TRUE.
ENDDO
ENDIF
! Insert particles in ducts
IF (DUCT_HT) CALL INSERT_DUCT_PARTICLES
! Reset particle/PARTICLE insertion clocks
DO N=1,N_INIT
IN => INITIALIZATION(N)
IF (IN%SINGLE_INSERTION) CYCLE
IF (T >= IN%PARTICLE_INSERT_CLOCK(NM)) IN%PARTICLE_INSERT_CLOCK(NM) = IN%PARTICLE_INSERT_CLOCK(NM) + IN%DT_INSERT
IF (T >= IN%PARTICLE_INSERT_CLOCK(NM)) INSERT_ANOTHER_BATCH = .TRUE.
ENDDO
DO N=1,N_SURF
SF => SURFACE(N)
ILPC = SF%PART_INDEX
IF (ILPC < 1 .AND. SF%N_LPC == 0) CYCLE
IF (T >= SF%PARTICLE_INSERT_CLOCK(NM)) SF%PARTICLE_INSERT_CLOCK(NM) = SF%PARTICLE_INSERT_CLOCK(NM) + SF%DT_INSERT
IF (T >= SF%PARTICLE_INSERT_CLOCK(NM)) INSERT_ANOTHER_BATCH = .TRUE.
ENDDO
IF (.NOT.INSERT_ANOTHER_BATCH) EXIT OVERALL_INSERT_LOOP
ENDDO OVERALL_INSERT_LOOP
! Compute initial particle CFL
IF (PARTICLE_CFL) THEN
DO IP=1,NLP
LP => MESHES(NM)%LAGRANGIAN_PARTICLE(IP)
LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
IF (LPC%MASSLESS_TRACER .OR. LPC%MASSLESS_TARGET) CYCLE
BC => MESHES(NM)%BOUNDARY_COORD(LP%BC_INDEX)
CALL GET_IJK(BC%X,BC%Y,BC%Z,NM,XI,YJ,ZK,II,JJ,KK)
PART_UVWMAX = MAX(PART_UVWMAX,MAX(ABS(LP%U)*RDX(II),ABS(LP%V)*RDY(JJ),ABS(LP%W)*RDZ(KK)))
ENDDO
ENDIF
! If any of the newly inserted particles finds itself in a neighboring mesh, set a flag to indicate that an MPI exchange must
! be done to transfer that particle to the MPI process that controls the neighboring mesh.
DO N=1,N_NEIGHBORING_MESHES
IF (OMESH(NEIGHBORING_MESH(N))%PARTICLE_SEND_BUFFER%N_ITEMS>0) EXCHANGE_INSERTED_PARTICLES = .TRUE.
ENDDO
T_USED(8)=T_USED(8)+CURRENT_TIME()-TNOW
CONTAINS
!> \brief Insert droplets for sprinklers and nozzles
SUBROUTINE INSERT_SPRAY_PARTICLES
INTEGER :: I,OI
! Loop over all devices, but look for sprinklers or nozzles. Count actuated sprinklers for output purposes.
N_ACTUATED_SPRINKLERS = 0
SPRINKLER_INSERT_LOOP: DO KS=1,N_DEVC
DV => DEVICE(KS)
PY => PROPERTY(DV%PROP_INDEX)
IF (PY%PART_ID == 'null') CYCLE SPRINKLER_INSERT_LOOP
IF (.NOT.DV%CURRENT_STATE) CYCLE SPRINKLER_INSERT_LOOP
IF (PY%QUANTITY=='SPRINKLER LINK TEMPERATURE') N_ACTUATED_SPRINKLERS = N_ACTUATED_SPRINKLERS + 1
IF (DV%MESH/=NM) CYCLE SPRINKLER_INSERT_LOOP
LPC=>LAGRANGIAN_PARTICLE_CLASS(PY%PART_INDEX)
IF (ABS(DV%T_CHANGE-T)<=TWO_EPSILON_EB) THEN
DV%T = T
CYCLE SPRINKLER_INSERT_LOOP
ENDIF
N_INSERT = NINT(REAL(PY%PARTICLES_PER_SECOND,EB)*(T-DV%T))
IF (N_INSERT<1) CYCLE SPRINKLER_INSERT_LOOP
! Determine sprinkler/nozzle flow rate
IF (ABS(T_BEGIN-DV%T_CHANGE)<=TWO_EPSILON_EB .AND. PY%FLOW_RAMP_INDEX>=1) THEN
TSI = T
ELSE
TSI = T - DV%T_CHANGE
ENDIF
IF (PY%MASS_FLOW_RATE > 0._EB) THEN
FLOW_RATE = PY%MASS_FLOW_RATE ! kg/s
PIPE_PRESSURE = 1.0_EB
D_PRES_FACTOR = 1.0_EB
ELSE
IF (PY%PRESSURE_RAMP_INDEX>0) THEN
PIPE_PRESSURE = EVALUATE_RAMP(REAL(DEVC_PIPE_OPERATING(DV%PIPE_INDEX),EB),PY%PRESSURE_RAMP_INDEX)
D_PRES_FACTOR = (PY%OPERATING_PRESSURE/PIPE_PRESSURE)**ONTH
FLOW_RATE = PY%K_FACTOR*SQRT(PIPE_PRESSURE)
ELSE
PIPE_PRESSURE = PY%OPERATING_PRESSURE
D_PRES_FACTOR = 1.0_EB
FLOW_RATE = PY%FLOW_RATE
ENDIF
FLOW_RATE = FLOW_RATE*LPC%DENSITY/1000._EB/60._EB ! convert from L/min to kg/s
ENDIF
FLOW_RATE = EVALUATE_RAMP(TSI,PY%FLOW_RAMP_INDEX,TAU=PY%FLOW_TAU)*FLOW_RATE ! kg/s
IF (FLOW_RATE <= TWO_EPSILON_EB) THEN
DV%T = T
CYCLE SPRINKLER_INSERT_LOOP
ENDIF
! Direction initialization stuff
OI = DV%ORIENTATION_INDEX
TRIGT1 = ACOS(-ORIENTATION_VECTOR(3,OI))
IF (ABS(ORIENTATION_VECTOR(2,OI))<=TWO_EPSILON_EB) THEN
TRIGT2 = ACOS(1._EB)
ELSE
TRIGT2 = ACOS(ABS(ORIENTATION_VECTOR(1,OI))/SQRT(ORIENTATION_VECTOR(1,OI)**2+ORIENTATION_VECTOR(2,OI)**2))
ENDIF
! PARTICLE insertion loop
MASS_SUM = 0._EB
DROP_SUM = 0
ALLOCATE(LP_INDEX_LOOKUP(N_INSERT))
LP_INDEX_LOOKUP = 0
PARTICLE_INSERT_LOOP: DO I=1,N_INSERT
! Insert a single particle
IF (NLP+1>MAXIMUM_PARTICLES) THEN
CALL REMOVE_OLDEST_PARTICLE(NM,PY%PART_INDEX,NLP,NEW_LP_INDEX)
IF (I>1) LP_INDEX_LOOKUP(I-1) = NEW_LP_INDEX
ELSE
NLP = NLP+1
ENDIF
LP_INDEX_LOOKUP(I) = NLP
CALL ALLOCATE_STORAGE(NM,LP_INDEX=NLP,LPC_INDEX=PY%PART_INDEX,SURF_INDEX=LPC%SURF_INDEX)
LAGRANGIAN_PARTICLE => MESHES(NM)%LAGRANGIAN_PARTICLE
LP=>LAGRANGIAN_PARTICLE(NLP)
LP%CLASS_INDEX = PY%PART_INDEX
BC=>MESHES(NM)%BOUNDARY_COORD(LP%BC_INDEX)
! Set PARTICLE properties
PARTICLE_TAG = PARTICLE_TAG + NMESHES
LP%TAG = PARTICLE_TAG
LP%T_INSERT = T
! Randomly choose particle direction angles, theta and phi
CHOOSE_COORDS: DO
PICK_PATTERN: IF(PY%SPRAY_PATTERN_INDEX>0) THEN ! Use spray pattern table
TA => TABLES(PY%SPRAY_PATTERN_INDEX)
CALL RANDOM_NUMBER(RN)
FIND_ROW: DO II=1,TA%NUMBER_ROWS
IF (REAL(RN,EB)>PY%TABLE_ROW(II)) CYCLE FIND_ROW
EXIT FIND_ROW
ENDDO FIND_ROW
CALL RANDOM_NUMBER(RN)
ETA_MAX=0.5_EB*(COS(TA%TABLE_DATA(II,1))+1._EB)
ETA_MIN=0.5_EB*(COS(TA%TABLE_DATA(II,2))+1._EB)
ETA=ETA_MIN+(ETA_MAX-ETA_MIN)*REAL(RN,EB)
THETA_RN=ACOS(2._EB*ETA-1._EB)
CALL RANDOM_NUMBER(RN)
PHI_RN = TA%TABLE_DATA(II,3) + REAL(RN,EB)*(TA%TABLE_DATA(II,4)-TA%TABLE_DATA(II,3))
IF (PY%PRESSURE_RAMP_INDEX>0) THEN
PARTICLE_SPEED = PY%V_FACTOR(II)*SQRT(PIPE_PRESSURE)
ELSE
PARTICLE_SPEED = TA%TABLE_DATA(II,5)
ENDIF
ELSE PICK_PATTERN ! Use conical spray
CALL RANDOM_CHOICE(PY%SPRAY_LON_CDF(:),PY%SPRAY_LON,NDC2,PHI_RN)
ILAT=MINLOC(ABS(PY%SPRAY_LON-PHI_RN),1)-1
CALL RANDOM_CHOICE(PY%SPRAY_LAT_CDF(:,ILAT),PY%SPRAY_LAT,NDC2,THETA_RN)
IF (PY%PRESSURE_RAMP_INDEX>0) THEN
PARTICLE_SPEED = PY%V_FACTOR(1)*SQRT(PIPE_PRESSURE)
ELSE
PARTICLE_SPEED = PY%PARTICLE_VELOCITY
ENDIF
ENDIF PICK_PATTERN
PHI_RN = PHI_RN + DV%ROTATION ! Adjust for rotation of head by rotating about z-axis
! Adjust for tilt of sprinkler pipe
SPHI = SIN(PHI_RN)
CPHI = COS(PHI_RN)
STHETA = SIN(THETA_RN)
CTHETA = COS(THETA_RN)
XTMP = STHETA*CPHI
YTMP = STHETA*SPHI
ZTMP = -CTHETA
! First rotate about y-axis away from x-axis
VLEN = SQRT(XTMP**2+ZTMP**2)
SHIFT1 = ACOS(ABS(XTMP)/VLEN)
SELECT CASE(INT(SIGN(1._EB,ZTMP)))
CASE (-1)
IF (XTMP<0) SHIFT1 = PI-SHIFT1
CASE ( 1)
SELECT CASE(INT(SIGN(1._EB,XTMP)))
CASE (-1)
SHIFT1 = SHIFT1+PI
CASE ( 1)
SHIFT1 = TWOPI - SHIFT1
END SELECT
END SELECT
SHIFT1 = SHIFT1 + TRIGT1
XTMP = VLEN * COS(SHIFT1)
ZTMP = -VLEN * SIN(SHIFT1)
! Second rotate about z-axis away from x-axis
VLEN = SQRT(XTMP**2+YTMP**2)
SHIFT1 = ACOS(ABS(XTMP)/VLEN)
SELECT CASE(INT(SIGN(1._EB,YTMP)))
CASE ( 1)
IF (XTMP<0) SHIFT1 = PI-SHIFT1
CASE (-1)
SELECT CASE(INT(SIGN(1._EB,XTMP)))
CASE (-1)
SHIFT1 = SHIFT1+PI
CASE ( 1)
SHIFT1 = TWOPI - SHIFT1
END SELECT
END SELECT
SHIFT2 = TRIGT2
OI = DV%ORIENTATION_INDEX
SELECT CASE(INT(SIGN(1._EB,ORIENTATION_VECTOR(1,OI))))
CASE (-1)
IF (ORIENTATION_VECTOR(2,OI)>0) SHIFT2 = TWOPI - SHIFT2
CASE ( 1)
SELECT CASE(INT(SIGN(1._EB,ORIENTATION_VECTOR(2,OI))))
CASE (-1)
SHIFT2 = PI-SHIFT2
CASE ( 1)
SHIFT2 = SHIFT2+ PI
END SELECT
END SELECT
SHIFT1=SHIFT1+SHIFT2
XTMP = VLEN * COS(SHIFT1)
YTMP = VLEN * SIN(SHIFT1)
! Compute initial position and velocity of PARTICLEs
LP%U = PARTICLE_SPEED*XTMP
LP%V = PARTICLE_SPEED*YTMP
LP%W = PARTICLE_SPEED*ZTMP
BC%X = DV%X + PY%OFFSET*XTMP
BC%Y = DV%Y + PY%OFFSET*YTMP
BC%Z = DV%Z + PY%OFFSET*ZTMP
IF (TWO_D) THEN
LP%V = 0._EB
BC%Y = DV%Y
ENDIF
! If the particle position is outside the current mesh, exit the loop and the particle will be sent to another mesh
! or eliminated by the call to REMOVE_PARTICLES at the end of the subroutine.
IF (BC%X<=XS .OR. BC%X>=XF .OR. BC%Y<=YS .OR. BC%Y>=YF .OR. BC%Z<=ZS .OR. BC%Z>=ZF) THEN
EXIT CHOOSE_COORDS
ELSE
CALL GET_IJK(BC%X,BC%Y,BC%Z,NM,XI,YJ,ZK,II,JJ,KK)
IC = CELL_INDEX(II,JJ,KK)
BC%IIG = II
BC%JJG = JJ
BC%KKG = KK
IF (.NOT.CELL(IC)%SOLID) EXIT CHOOSE_COORDS
ENDIF
ENDDO CHOOSE_COORDS
! Randomly choose PARTICLE size according to Cumulative Distribution Function (CDF)
CALL INITIALIZE_SINGLE_PARTICLE
! Adjust particle size to account for pressure dependence of nozzle
IF (LPC%LIQUID_DROPLET) THEN
LP%RADIUS = LP%RADIUS*D_PRES_FACTOR
LP%MASS = LP%MASS*D_PRES_FACTOR**3
ENDIF
! Sum up mass of liquid being introduced
MASS_SUM = MASS_SUM + LP%PWT*LP%MASS
DROP_SUM = DROP_SUM + 1
ENDDO PARTICLE_INSERT_LOOP
! Compute weighting factor for the PARTICLEs just inserted
IF (DROP_SUM > 0) THEN
IF (LPC%LIQUID_DROPLET) THEN
PWT0 = LPC%N_STRATA*FLOW_RATE/(LPC%DENSITY*LPC%MEAN_DROPLET_VOLUME*REAL(PY%PARTICLES_PER_SECOND,EB))/D_PRES_FACTOR**3
ELSE
PWT0 = FLOW_RATE*REAL(DROP_SUM,EB)/(MASS_SUM*REAL(PY%PARTICLES_PER_SECOND,EB))
ENDIF
DO I=1,N_INSERT
N = LP_INDEX_LOOKUP(I)
LAGRANGIAN_PARTICLE(N)%PWT = LAGRANGIAN_PARTICLE(N)%PWT * PWT0
ENDDO
ENDIF
DEALLOCATE(LP_INDEX_LOOKUP)
! Indicate that PARTICLEs from this device have been inserted at this time T
DV%T = T
ENDDO SPRINKLER_INSERT_LOOP
! There might be particles that were inserted outside the current mesh. Either remove them or transfer them to another mesh.
CALL REMOVE_PARTICLES(T,NM)
END SUBROUTINE INSERT_SPRAY_PARTICLES
!> \brief Loop through all boundary cells and insert particles if appropriate
SUBROUTINE INSERT_VENT_PARTICLES
DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
CALL PARTICLE_FACE_INSERT(WALL_INDEX=IW)
ENDDO
DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
CALL PARTICLE_FACE_INSERT(CFACE_INDEX=ICF)
ENDDO
END SUBROUTINE INSERT_VENT_PARTICLES
SUBROUTINE PARTICLE_FACE_INSERT(WALL_INDEX,CFACE_INDEX)
USE COMPLEX_GEOMETRY, ONLY : RANDOM_CFACE_XYZ
USE PHYSICAL_FUNCTIONS, ONLY: GET_PARTICLE_ENTHALPY
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX
INTEGER :: I,N_LPC,N_LPC_MAX,ITER,LPC_INDEX,INSERT_TYPE
REAL(EB):: CFA_X, CFA_Y, CFA_Z, RN, VEL_PART, C_S, H_1, H_2, TMP_PART, TMP_GUESS
LOGICAL :: RETURN_FLAG
TYPE (CFACE_TYPE), POINTER :: CFA
WALL_OR_CFACE_IF_1: IF (PRESENT(WALL_INDEX)) THEN
WC => MESHES(NM)%WALL(WALL_INDEX)
IF (WC%OD_INDEX>0) ONE_D => MESHES(NM)%BOUNDARY_ONE_D(WC%OD_INDEX)
B1 => MESHES(NM)%BOUNDARY_PROP1(WC%B1_INDEX)
BC => MESHES(NM)%BOUNDARY_COORD(WC%BC_INDEX)
SF => SURFACE(WC%SURF_INDEX)
IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) RETURN
II = BC%II
JJ = BC%JJ
KK = BC%KK
IC = CELL_INDEX(II,JJ,KK)
IF (.NOT.CELL(IC)%SOLID) RETURN
IOR = BC%IOR
ELSEIF (PRESENT(CFACE_INDEX)) THEN
CFA => MESHES(NM)%CFACE(CFACE_INDEX)
IF (CFA%OD_INDEX>0) ONE_D => MESHES(NM)%BOUNDARY_ONE_D(CFA%OD_INDEX)
B1 => MESHES(NM)%BOUNDARY_PROP1(CFA%B1_INDEX)
BC => MESHES(NM)%BOUNDARY_COORD(CFA%BC_INDEX)
SF => SURFACE(CFA%SURF_INDEX)
IF (CFA%BOUNDARY_TYPE/=SOLID_BOUNDARY) RETURN
ENDIF WALL_OR_CFACE_IF_1
ILPC = SF%PART_INDEX; IF (ILPC < 1 .AND. SF%N_LPC ==0) RETURN
IIG = BC%IIG
JJG = BC%JJG
KKG = BC%KKG
RETURN_FLAG = .FALSE.
IF (T < SF%PARTICLE_INSERT_CLOCK(NM)) RETURN
INSERT_TYPE_LOOP: DO INSERT_TYPE = 1,2
SELECT CASE (INSERT_TYPE)
CASE (1) ! PART_ID on SURF
IF (ILPC <=0) CYCLE INSERT_TYPE_LOOP
IF (T < B1%T_IGN) CYCLE INSERT_TYPE_LOOP
IF (SF%PARTICLE_SURFACE_DENSITY>0._EB .AND. T>T_BEGIN) CYCLE INSERT_TYPE_LOOP
IF (ANY(SF%EMBER_GENERATION_HEIGHT>=0._EB)) THEN
! specify generation only for regions of burning
IF (.NOT. B1%M_DOT_G_PP_ADJUST(REACTION(1)%FUEL_SMIX_INDEX)>0._EB) CYCLE INSERT_TYPE_LOOP
ENDIF
! If there is a flow ramp and the value is zero, return.
IF (SF%PARTICLE_MASS_FLUX > 0._EB) THEN
IF (ABS(B1%T_IGN-T_BEGIN)<=TWO_EPSILON_EB .AND. SF%RAMP(TIME_PART)%INDEX>=1) THEN
TSI = T
ELSE
TSI = T - B1%T_IGN
ENDIF
FLOW_RATE = EVALUATE_RAMP(TSI,SF%RAMP(TIME_PART)%INDEX,TAU=SF%RAMP(TIME_PART)%TAU)*SF%PARTICLE_MASS_FLUX
IF (FLOW_RATE < TWO_EPSILON_EB) RETURN
ENDIF
LPC => LAGRANGIAN_PARTICLE_CLASS(ILPC)
! Evalutate if we need to skip ILPC and return if no N_LPC
IF (LPC%DEVC_INDEX>0) THEN
IF (.NOT.DEVICE(LPC%DEVC_INDEX)%CURRENT_STATE) RETURN_FLAG = .TRUE.
ENDIF
IF (LPC%CTRL_INDEX>0) THEN
IF (.NOT.CONTROL(LPC%CTRL_INDEX)%CURRENT_STATE) RETURN_FLAG = .TRUE.
ENDIF
IF (NM > 1) THEN
IF (INTERPOLATED_MESH(IIG,JJG,KKG) > 0) RETURN_FLAG = .TRUE.
ENDIF
IF (RETURN_FLAG) CYCLE INSERT_TYPE_LOOP
N_LPC_MAX = 1
CASE (2) ! PART_ID on MATL
IF (SF%N_LPC==0) RETURN
N_LPC_MAX = SF%N_LPC
END SELECT
ALLOCATE(LP_INDEX_LOOKUP(SF%NPPC))
LPC_LOOP: DO N_LPC=1,N_LPC_MAX
LP_INDEX_LOOKUP = 0
IF (INSERT_TYPE==1) THEN
LPC_INDEX = ILPC
ELSE
LPC_INDEX = SF%MATL_PART_INDEX(N_LPC)
ENDIF
! Get particle temperature for MATL particle
IF (INSERT_TYPE==2) THEN
IF (ONE_D%PART_MASS(N_LPC) < TWO_EPSILON_EB) CYCLE
LPC => LAGRANGIAN_PARTICLE_CLASS(SF%MATL_PART_INDEX(N_LPC))
TMP_GUESS = ONE_D%TMP(1)
H_1 = ONE_D%PART_ENTHALPY(N_LPC)/ONE_D%PART_MASS(N_LPC)
C_S = H_1/TMP_GUESS
ITER = 0
T_SEARCH: DO
ITER = ITER + 1
H_2 = GET_PARTICLE_ENTHALPY(SF%MATL_PART_INDEX(N_LPC),TMP_GUESS)
C_S = H_2/TMP_GUESS
TMP_PART = TMP_GUESS - (H_2 - H_1)/C_S
IF (ABS(TMP_PART - TMP_GUESS) < 0.0001_EB) EXIT T_SEARCH
IF (ITER > 20) THEN
TMP_PART = 0.5_EB*(TMP_PART+TMP_GUESS)
EXIT T_SEARCH
ENDIF
TMP_GUESS = TMP_PART
ENDDO T_SEARCH
ENDIF
! Loop over all particles for the WALL_INDEX-th cell
MASS_SUM = 0._EB
PARTICLE_INSERT_LOOP: DO I=1,SF%NPPC
! Insert a single droplet at wall cell WALL_INDEX or cut-cell face CFACE_INDEXF
IF (NLP+1>MAXIMUM_PARTICLES) THEN
CALL REMOVE_OLDEST_PARTICLE(NM,LPC_INDEX,NLP,NEW_LP_INDEX)
IF (I>1) LP_INDEX_LOOKUP(I-1)=NEW_LP_INDEX
ELSE
NLP = NLP+1
ENDIF
LP_INDEX_LOOKUP(I) = NLP
CALL ALLOCATE_STORAGE(NM,LP_INDEX=NLP,LPC_INDEX=LPC_INDEX,SURF_INDEX=LPC%SURF_INDEX)
LAGRANGIAN_PARTICLE => MESHES(NM)%LAGRANGIAN_PARTICLE
LP => LAGRANGIAN_PARTICLE(NLP)
LP%CLASS_INDEX = LPC_INDEX
PARTICLE_TAG = PARTICLE_TAG + NMESHES
LP%TAG = PARTICLE_TAG
! Reassign pointers after calling ALLOCATE
BC => MESHES(NM)%BOUNDARY_COORD(LP%BC_INDEX)
IF (PRESENT(WALL_INDEX)) THEN
WC => MESHES(NM)%WALL(WALL_INDEX)
B1 => MESHES(NM)%BOUNDARY_PROP1(WC%B1_INDEX)
IF (WC%OD_INDEX>0) ONE_D => MESHES(NM)%BOUNDARY_ONE_D(WC%OD_INDEX)
ELSEIF (PRESENT(CFACE_INDEX)) THEN
CFA => MESHES(NM)%CFACE(CFACE_INDEX)
B1 => MESHES(NM)%BOUNDARY_PROP1(CFA%B1_INDEX)
IF (CFA%OD_INDEX>0) ONE_D => MESHES(NM)%BOUNDARY_ONE_D(CFA%OD_INDEX)
ENDIF
! Ember flag to be used for outputs
IF (INSERT_TYPE==1 .AND. ANY(SF%EMBER_GENERATION_HEIGHT>=0._EB)) LP%EMBER=.TRUE.
! Assign particle position on the cell face
CALL RANDOM_NUMBER(RN)
CALL RANDOM_NUMBER(RN2)
WALL_OR_CFACE_IF_2: IF (PRESENT(WALL_INDEX)) THEN
SELECT CASE (ABS(IOR))
CASE(1)
IF (IOR== 1) BC%X = X(II) + VENT_OFFSET*DX(II+1)
IF (IOR==-1) BC%X = X(II-1) - VENT_OFFSET*DX(II-1)
BC%Y = Y(JJ-1) + DY(JJ)*REAL(RN,EB)
BC%Z = Z(KK-1) + DZ(KK)*REAL(RN2,EB)
CASE(2)
IF (IOR== 2) BC%Y = Y(JJ) + VENT_OFFSET*DY(JJ+1)
IF (IOR==-2) BC%Y = Y(JJ-1) - VENT_OFFSET*DY(JJ-1)
BC%X = X(II-1) + DX(II)*REAL(RN,EB)
BC%Z = Z(KK-1) + DZ(KK)*REAL(RN2,EB)
CASE(3)
IF (IOR== 3) THEN
BC%Z = Z(KK) + VENT_OFFSET*DZ(KK+1)
IF (ANY(SF%EMBER_GENERATION_HEIGHT>=0._EB)) THEN
CALL RANDOM_NUMBER(RN3)
BC%Z = Z(KK) + SF%EMBER_GENERATION_HEIGHT(1) + &
(SF%EMBER_GENERATION_HEIGHT(2)-SF%EMBER_GENERATION_HEIGHT(1))*REAL(RN3,EB)
ENDIF
ENDIF
IF (IOR==-3) BC%Z = Z(KK-1) - VENT_OFFSET*DZ(KK-1)
BC%X = X(II-1) + DX(II)*REAL(RN,EB)
BC%Y = Y(JJ-1) + DY(JJ)*REAL(RN2,EB)
END SELECT
! Give particles an initial velocity
IF (.NOT.LPC%STATIC) THEN
IF (SF%VEL_PART >-999999._EB) THEN
VEL_PART = SF%VEL_PART
ELSE
VEL_PART = B1%U_NORMAL
ENDIF
SELECT CASE(IOR)
CASE( 1)
LP%U = -VEL_PART
LP%V = SF%VEL_T(1)
LP%W = SF%VEL_T(2)
CASE(-1)
LP%U = VEL_PART
LP%V = SF%VEL_T(1)
LP%W = SF%VEL_T(2)
CASE( 2)
LP%U = SF%VEL_T(1)
LP%V = -VEL_PART
LP%W = SF%VEL_T(2)
CASE(-2)
LP%U = SF%VEL_T(1)
LP%V = VEL_PART
LP%W = SF%VEL_T(2)
CASE( 3)
LP%U = SF%VEL_T(1)
LP%V = SF%VEL_T(2)
LP%W = -VEL_PART
CASE(-3)
LP%U = SF%VEL_T(1)
LP%V = SF%VEL_T(2)
LP%W = VEL_PART
END SELECT
ENDIF
ELSEIF (PRESENT(CFACE_INDEX)) THEN
CALL RANDOM_CFACE_XYZ(NM,CFA,CFA_X,CFA_Y,CFA_Z)
BC%X = CFA_X + CFA%NVEC(1)*VENT_OFFSET*DX(IIG)
BC%Y = CFA_Y + CFA%NVEC(2)*VENT_OFFSET*DY(JJG)
BC%Z = CFA_Z + CFA%NVEC(3)*VENT_OFFSET*DZ(KKG)
IF (INSERT_TYPE==1 .AND. LP%EMBER) THEN
CALL RANDOM_NUMBER(RN3)
BC%Z = CFA_Z + SF%EMBER_GENERATION_HEIGHT(1) + &
(SF%EMBER_GENERATION_HEIGHT(2)-SF%EMBER_GENERATION_HEIGHT(1))*REAL(RN3,EB)
ENDIF
LP%U = DOT_PRODUCT(CFA%NVEC,(/-B1%U_NORMAL,SF%VEL_T(1),SF%VEL_T(2)/))
LP%V = DOT_PRODUCT(CFA%NVEC,(/SF%VEL_T(1),-B1%U_NORMAL,SF%VEL_T(2)/))
LP%W = DOT_PRODUCT(CFA%NVEC,(/SF%VEL_T(1),SF%VEL_T(2),-B1%U_NORMAL/))
ENDIF WALL_OR_CFACE_IF_2
! Embers may not be generated in wall-adjacent cell
IF (INSERT_TYPE==1 .AND. LP%EMBER) THEN
CALL GET_IJK(BC%X,BC%Y,BC%Z,NM,XI,YJ,ZK,IIG,JJG,KKG)
BC%IIG = IIG
BC%JJG = JJG
BC%KKG = KKG
ELSE
BC%IIG = IIG
BC%JJG = JJG
BC%KKG = KKG
ENDIF
! Save the insertion time (TP) and scalar property (SP) for the particle
LP%T_INSERT = T
CALL INITIALIZE_SINGLE_PARTICLE
IF (INSERT_TYPE==2) MESHES(NM)%BOUNDARY_ONE_D(LP%OD_INDEX)%TMP = TMP_PART
IF (.NOT.LPC%MASSLESS_TRACER .AND. .NOT.LPC%MASSLESS_TARGET) MASS_SUM = MASS_SUM + LP%PWT*LP%MASS
ENDDO PARTICLE_INSERT_LOOP
! Adjust the particle weighting factors to get the right mass flux
MASS_CHECK: IF (MASS_SUM > 0._EB) THEN
SELECT CASE (INSERT_TYPE)
CASE (1) ! PART_ID on SURF
IF (MASS_SUM > 0._EB) THEN
IF (SF%PARTICLE_MASS_FLUX > 0._EB) THEN
DO I=1,SF%NPPC
LP => LAGRANGIAN_PARTICLE(LP_INDEX_LOOKUP(I))
LP%PWT = LP%PWT * FLOW_RATE*B1%AREA_ADJUST*B1%AREA*SF%DT_INSERT/MASS_SUM
ENDDO
ELSEIF (SF%PARTICLE_SURFACE_DENSITY > 0._EB) THEN
DO I=1,SF%NPPC
LP => LAGRANGIAN_PARTICLE(LP_INDEX_LOOKUP(I))
LP%PWT = LP%PWT * SF%PARTICLE_SURFACE_DENSITY*B1%AREA_ADJUST*B1%AREA/MASS_SUM
ENDDO
ENDIF
ENDIF
CASE (2) ! PART_ID on MATL
IF (MASS_SUM > 0._EB) THEN
DO I=1,SF%NPPC
LP => LAGRANGIAN_PARTICLE(LP_INDEX_LOOKUP(I))
LP%PWT = LP%PWT * ONE_D%PART_MASS(N_LPC)*B1%AREA_ADJUST*B1%AREA* &
MIN(1._EB,SF%DT_INSERT/B1%T_MATL_PART)/MASS_SUM
ENDDO
ENDIF
END SELECT
ENDIF MASS_CHECK
ENDDO LPC_LOOP
IF(ALLOCATED(LP_INDEX_LOOKUP)) DEALLOCATE(LP_INDEX_LOOKUP)
IF (INSERT_TYPE==2) THEN
! Decrement mass, enthalpy, and aggregation time for current insertion interval
ONE_D%PART_MASS = MAX(0._EB,1-SF%DT_INSERT/(B1%T_MATL_PART+TINY_EB))*ONE_D%PART_MASS
ONE_D%PART_ENTHALPY = MAX(0._EB,1-SF%DT_INSERT/(B1%T_MATL_PART+TINY_EB))*ONE_D%PART_ENTHALPY
B1%T_MATL_PART = MAX(0._EB,B1%T_MATL_PART-SF%DT_INSERT)
ENDIF
ENDDO INSERT_TYPE_LOOP
END SUBROUTINE PARTICLE_FACE_INSERT
!> \brief Loop over all INIT lines and look for particles inserted within a specified volume
SUBROUTINE INSERT_VOLUMETRIC_PARTICLES
USE COMPLEX_GEOMETRY, ONLY: CC_IDCC,CC_CGSC,CC_SOLID
INTEGER :: IIP,N_INSERT,I1,J1,K1,I2,J2,K2,N,N_PARTICLES_INSERT,ND,ICC
REAL(EB) :: XC1,XC2,YC1,YC2,ZC1,ZC2,X0,Y0,Z0,RR,RRI,HH,INSERT_VOLUME,INPUT_VOLUME,VOLUME_SPLIT_FACTOR,LP_X,LP_Y,LP_Z,RAMP_FACTOR,&
IN_X1,IN_X2,IN_Y1,IN_Y2,IN_Z1,IN_Z2,IN_X0,IN_Y0,IN_Z0,VCX,VCY,VCZ,MOIST_FRAC
LOGICAL :: CC_VALID
IN => INITIALIZATION(INIT_INDEX)
! Determine if the INITIALIZATION type involves particles. If not, return.
ILPC = IN%PART_INDEX
IF (ILPC<1) RETURN
IF (IN%SINGLE_INSERTION .AND. IN%ALREADY_INSERTED(NM)) RETURN
! If there is a RAMP for MASS_PER_TIME or MASS_PER_VOLUME, evaluate it now and if zero, return.
IF (IN%RAMP_PART_INDEX>0) THEN
RAMP_FACTOR = EVALUATE_RAMP(T,IN%RAMP_PART_INDEX)
IF (RAMP_FACTOR<TWO_EPSILON_EB) RETURN
ELSE
RAMP_FACTOR = 1._EB
ENDIF
! Determine if the particles/PARTICLEs are controlled by devices
LPC => LAGRANGIAN_PARTICLE_CLASS(ILPC)
IF (IN%DEVC_INDEX>0) THEN
IF (.NOT.DEVICE(IN%DEVC_INDEX)%CURRENT_STATE) THEN
IN%PARTICLE_INSERT_CLOCK(NM) = T
RETURN
ENDIF
ENDIF
IF (IN%CTRL_INDEX>0) THEN
IF (.NOT.CONTROL(IN%CTRL_INDEX)%CURRENT_STATE) THEN
IN%PARTICLE_INSERT_CLOCK(NM) = T
RETURN
ENDIF
ENDIF
! If it is not time to insert particles for this INITIALIZATION block, cycle.
IF (T < IN%PARTICLE_INSERT_CLOCK(NM)) RETURN
! Start processing the INITIALIZATION info
IF (IN%N_PARTICLES==0 .AND. IN%N_PARTICLES_PER_CELL==0) RETURN
! Adjust INIT values for voxelized bulk density data
IF (IN%BULK_DENSITY_FILE/='null') THEN
! Center of voxel
READ(LU_VEG_IN) VCX,VCY,VCZ
IN%X1=VCX-VDX/2._EB
IN%X2=VCX+VDX/2._EB
IN%Y1=VCY-VDY/2._EB
IN%Y2=VCY+VDY/2._EB
IN%Z1=VCZ-VDZ/2._EB
IN%Z2=VCZ+VDZ/2._EB
! Vegetation mass in voxel
READ(LU_VEG_IN) IN%MASS_PER_VOLUME
! Adjust with multiplier
IN%MASS_PER_VOLUME = MAX(0._EB,IN%BULK_DENSITY_FACTOR)*IN%MASS_PER_VOLUME
MOIST_FRAC=SURFACE(LAGRANGIAN_PARTICLE_CLASS(IN%PART_INDEX)%SURF_INDEX)%MOISTURE_FRACTION(1)
IF (MOIST_FRAC>=0._EB) IN%MASS_PER_VOLUME = IN%MASS_PER_VOLUME*(1._EB+MOIST_FRAC)
ENDIF
! Apply coordinate offset if needed
IN_X1 = X_OFFSET + IN%X1
IN_X2 = X_OFFSET + IN%X2
IN_Y1 = Y_OFFSET + IN%Y1
IN_Y2 = Y_OFFSET + IN%Y2
IN_Z1 = Z_OFFSET + IN%Z1
IN_Z2 = Z_OFFSET + IN%Z2
IN_X0 = X_OFFSET + IN%X0
IN_Y0 = X_OFFSET + IN%Y0
IN_Z0 = X_OFFSET + IN%Z0
! Cut off parts of the INIT region that are outside the current mesh
IF (IN_X1>XF .OR. IN_X2<XS .OR. IN_Y1>YF .OR. IN_Y2<YS .OR. IN_Z1>ZF .OR. IN_Z2<ZS) RETURN
! Skip mesh that is contained completely within a ring
IF (IN%SHAPE=='RING' .AND. IN_X1<XS .AND. IN_X2>XF .AND. IN_Y1<YS .AND. IN_Y2>YF .AND. IN_Z1<ZS .AND. IN_Z2>ZF) RETURN
X1 = MAX(IN_X1,XS)