-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpart.f90
2151 lines (1749 loc) · 73.9 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
MODULE PART
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE TRAN
IMPLICIT NONE
PRIVATE
PUBLIC INSERT_PARTICLES,UPDATE_PARTICLES,REMOVE_PARTICLES,INITIALIZE_PARTICLES,GET_REV_part
CHARACTER(255), PARAMETER :: partid='$Id: part.f90 9953 2012-02-01 12:57:01Z Topi.Sikanen $'
CHARACTER(255), PARAMETER :: partrev='$Revision: 9953 $'
CHARACTER(255), PARAMETER :: partdate='$Date: 2012-02-01 04:57:01 -0800 (Wed, 01 Feb 2012) $'
REAL(EB) :: INSERT_RATE = 0._EB
INTEGER :: INSERT_COUNT = 0
CONTAINS
SUBROUTINE INITIALIZE_PARTICLES(NM)
! Insert PARTICLEs into the domain at the start of calculation
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE COMP_FUNCTIONS, ONLY : SECOND
USE PHYSICAL_FUNCTIONS, ONLY : PARTICLE_SIZE_DISTRIBUTION
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_PARTICLES
REAL(EB) :: LL,UL,BIN_SIZE,TNOW,DD,DI
INTEGER :: I,J,IL,IU,IPC
INTEGER, INTENT(IN) :: NM
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC=>NULL()
TYPE (RAMPS_TYPE), POINTER :: RM=>NULL()
IF (N_PART==0) RETURN ! Don't waste time if no particles
IF (EVACUATION_ONLY(NM)) RETURN ! Don't waste time if an evac mesh
TNOW=SECOND()
CALL POINT_TO_MESH(NM)
PART_CLASS_LOOP: DO IPC=1,N_PART
PC=>PARTICLE_CLASS(IPC)
! If particles/PARTICLEs have a size distribution, initialize here
IF_SIZE_DISTRIBUTION: IF (.NOT.PC%MONODISPERSE .AND. (PC%DIAMETER > 0._EB)) THEN
IF(PC%CNF_RAMP_INDEX<0) THEN
CALL PARTICLE_SIZE_DISTRIBUTION(PC%DIAMETER,PC%R_CDF(:),PC%CDF(:),NDC,PC%GAMMA,PC%SIGMA,PC%DISTRIBUTION)
ELSE
RM=>RAMPS(PC%CNF_RAMP_INDEX)
DD=RM%SPAN/NDC
DO I=1,NDC
DI=RM%T_MIN+(I-0.5_EB)*DD
PC%R_CDF(I) = 0.5_EB*DI
PC%CDF(I) = EVALUATE_RAMP(DI,0._EB,PC%CNF_RAMP_INDEX)
ENDDO
PC%CDF=PC%CDF/PC%CDF(NDC)
ENDIF
BIN_SIZE = PC%R_CDF(NDC)/REAL(NSTRATA,EB)
STRATIFY: DO I=1,NSTRATA
LL = (I-1)*BIN_SIZE
UL = I *BIN_SIZE
IL = 1
LL_LOOP: DO J=1,NDC
IF (PC%R_CDF(J)>LL) THEN
IL = J-1
PC%IL_CDF(I) = J-1
EXIT LL_LOOP
ENDIF
ENDDO LL_LOOP
IU = NDC
UL_LOOP: DO J=NDC,1,-1
IF (PC%R_CDF(J)<=UL) THEN
IU = J
PC%IU_CDF(I) = J
EXIT UL_LOOP
ENDIF
ENDDO UL_LOOP
PC%W_CDF(I) = PC%CDF(IU) - PC%CDF(IL)
ENDDO STRATIFY
ENDIF IF_SIZE_DISTRIBUTION
! If pacticles/PARTICLEs can break up, compute normalized (median = 1) size distribution for child PARTICLEs
IF (PC%BREAKUP .AND. .NOT.PC%MONODISPERSE) THEN
IF(PC%CHILD_CNF_RAMP_INDEX<0) THEN
CALL PARTICLE_SIZE_DISTRIBUTION(1._EB,PC%CHILD_R_CDF(:),PC%CHILD_CDF(:),NDC, &
PC%BREAKUP_CHILD_GAMMA,PC%BREAKUP_CHILD_SIGMA,PC%CHILD_DISTRIBUTION)
ELSE
RM=>RAMPS(PC%CHILD_CNF_RAMP_INDEX)
DD=RM%SPAN/NDC
DO I=1,NDC
DI=RM%T_MIN+(I-0.5_EB)*DD
PC%CHILD_R_CDF(I) = 0.5_EB*DI
PC%CHILD_CDF(I) = EVALUATE_RAMP(DI,0._EB,PC%CHILD_CNF_RAMP_INDEX)
ENDDO
PC%CHILD_CDF=PC%CHILD_CDF/PC%CHILD_CDF(NDC)
ENDIF
ENDIF
ENDDO PART_CLASS_LOOP
TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE INITIALIZE_PARTICLES
SUBROUTINE INSERT_PARTICLES(T,NM)
! Insert sprinkler PARTICLEs and lagrangian particles into the domain every time step
USE COMP_FUNCTIONS, ONLY : SECOND
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_PARTICLES
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,RANDOM_CHOICE
USE GEOMETRY_FUNCTIONS, ONLY: RANDOM_RECTANGLE,RANDOM_CONE
USE TRAN, ONLY: GET_IJK
USE DEVICE_VARIABLES
USE CONTROL_VARIABLES
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
REAL :: RN
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,MASS_PER_TIME,MASS_PER_VOLUME,X1,X2,Y1,Y2,Z1,Z2,BLOCK_VOLUME, &
ETA,ETA_MAX,ETA_MIN,XI,YJ,ZK
REAL(EB), PARAMETER :: VENT_OFFSET=0.1
INTEGER :: I,KS,II,JJ,KK,IC,IL,IU,IPC,DROP_SUM,IIG,JJG,KKG,IW,IOR,STRATUM,IB
INTEGER :: N_INITIAL,N,N_INSERT,ILAT
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 (PARTICLE_TYPE), POINTER :: LP=>NULL()
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC=>NULL()
TYPE (INITIALIZATION_TYPE), POINTER :: IN=>NULL()
IF (EVACUATION_ONLY(NM)) RETURN ! Don't waste time if an evac mesh
IF (N_PART==0) RETURN ! Don't waste time if no particles
TNOW=SECOND()
CALL POINT_TO_MESH(NM)
OVERALL_INSERT_LOOP: DO
INSERT_ANOTHER_BATCH = .FALSE.
CALL DEVICE_PARTICLE_INSERT
CALL WALL_PARTICLE_INSERT
CALL VOLUME_PARTICLE_INSERT
! 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)
IPC = SF%PART_INDEX
IF (IPC < 1) 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
TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
CONTAINS
SUBROUTINE DEVICE_PARTICLE_INSERT
! Count active sprinklers and nozzles
N_OPEN_NOZZLES = 0
N_ACTUATED_SPRINKLERS = 0
COUNT_OPEN_NOZZLES_LOOP: DO KS=1,N_DEVC ! Loop over all devices, but look for sprinklers or nozzles
DV => DEVICE(KS)
PY => PROPERTY(DV%PROP_INDEX)
IF (.NOT. DV%CURRENT_STATE) CYCLE COUNT_OPEN_NOZZLES_LOOP
IF (PY%PART_ID == 'null') CYCLE COUNT_OPEN_NOZZLES_LOOP
N_OPEN_NOZZLES = N_OPEN_NOZZLES + 1
IF (PY%QUANTITY=='SPRINKLER LINK TEMPERATURE') N_ACTUATED_SPRINKLERS = N_ACTUATED_SPRINKLERS + 1
ENDDO COUNT_OPEN_NOZZLES_LOOP
! Count the number of PARTICLEs/particles inserted
INSERT_RATE = 0._EB
INSERT_COUNT = 0
! Add more than one batch of particles/PARTICLEs if FDS time step is large enough
! Loop over all devices, but look for sprinklers or nozzles
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 (DV%MESH/=NM) CYCLE SPRINKLER_INSERT_LOOP
IF (.NOT. DV%CURRENT_STATE) CYCLE SPRINKLER_INSERT_LOOP
IPC = PY%PART_INDEX
PC=>PARTICLE_CLASS(IPC)
IF (ABS(DV%T_CHANGE-T)<=ZERO_P) THEN
DV%T = T
CYCLE SPRINKLER_INSERT_LOOP
ENDIF
! IF (T < PY%PARTICLE_INSERT_CLOCK(NM)) CYCLE SPRINKLER_INSERT_LOOP
INSERT_RATE = INSERT_RATE + PY%PARTICLES_PER_SECOND
N_INSERT = NINT(REAL(PY%PARTICLES_PER_SECOND,EB)*DT)
INSERT_COUNT = INSERT_COUNT + N_INSERT
! Determine sprinkler/nozzle flow rate
IF (ABS(T_BEGIN-DV%T_CHANGE)<=ZERO_P .AND. PY%FLOW_RAMP_INDEX>=1) THEN
TSI = T
ELSE
TSI = T - DV%T_CHANGE
ENDIF
IF (PY%PRESSURE_RAMP_INDEX>0) THEN
PIPE_PRESSURE = EVALUATE_RAMP(REAL(N_OPEN_NOZZLES,EB),0._EB,PY%PRESSURE_RAMP_INDEX)
D_PRES_FACTOR = (PY%OPERATING_PRESSURE/PIPE_PRESSURE)**(1._EB/3._EB)
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 = EVALUATE_RAMP(TSI,PY%FLOW_TAU,PY%FLOW_RAMP_INDEX)*FLOW_RATE*(PC%DENSITY/1000._EB)/60._EB ! kg/s
IF (FLOW_RATE <= 0._EB) THEN
DV%T = T
CYCLE SPRINKLER_INSERT_LOOP
ENDIF
! Direction initialization stuff
TRIGT1 = ACOS(-DV%ORIENTATION(3))
IF (ABS(DV%ORIENTATION(2))<=ZERO_P) THEN
TRIGT2 = ACOS(1._EB)
ELSE
TRIGT2 = ACOS(ABS(DV%ORIENTATION(1))/SQRT(DV%ORIENTATION(1)**2+DV%ORIENTATION(2)**2))
ENDIF
! PARTICLE insertion loop
MASS_SUM = 0._EB
DROP_SUM = 0
PARTICLE_INSERT_LOOP: DO I=1,N_INSERT
IF (NLP+1>MAXIMUM_PARTICLES) EXIT PARTICLE_INSERT_LOOP
NLP = NLP+1
PARTICLE_TAG = PARTICLE_TAG + NMESHES
IF (NLP>NLPDIM) THEN
CALL RE_ALLOCATE_PARTICLES(1,NM,0,1000)
PARTICLE=>MESHES(NM)%PARTICLE
ENDIF
LP=>PARTICLE(NLP)
! Set PARTICLE properties
LP%TMP = PC%TMP_INITIAL ! Initial temperature
LP%T = T ! Time of insertion
LP%IOR = 0 ! Orientation index (0 means it is not attached to a solid)
LP%ACCEL_X = 0._EB ! x component of drag on the gas in units of m/s2
LP%ACCEL_Y = 0._EB ! y component of drag on the gas in units of m/s2
LP%ACCEL_Z = 0._EB ! z component of drag on the gas in units of m/s2
LP%CLASS = IPC ! The class of particles
LP%TAG = PARTICLE_TAG ! A unique identifying integer
IF (MOD(NLP,PC%SAMPLING)==0) THEN ! Decide whether to show the PARTICLE in Smokeview
LP%SHOW = .TRUE.
ELSE
LP%SHOW = .FALSE.
ENDIF
LP%SPLAT = .FALSE.
LP%WALL_INDEX = 0
! Randomly choose 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
END DO FIND_ROW
CALL RANDOM_NUMBER(RN)
!THETA_RN = TA%TABLE_DATA(II,1) + RN*(TA%TABLE_DATA(II,2)-TA%TABLE_DATA(II,1))
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_NUMBER(RN)
!THETA_RN = PY%SPRAY_ANGLE(1) + REAL(RN,EB)*(PY%SPRAY_ANGLE(2)-PY%SPRAY_ANGLE(1))
CALL RANDOM_CHOICE(PY%SPRAY_LON_CDF(:),PY%SPRAY_LON,NDC2,PHI_RN)
ILAT=MINLOC(ABS(PY%SPRAY_LON-PHI_RN),1)
CALL RANDOM_CHOICE(PY%SPRAY_LAT_CDF(:,ILAT),PY%SPRAY_LAT,NDC2,THETA_RN)
!CALL RANDOM_NUMBER(RN)
!PHI_RN = RN*TWOPI
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
SELECT CASE(INT(SIGN(1._EB,DV%ORIENTATION(1))))
CASE (-1)
IF (DV%ORIENTATION(2)>0) SHIFT2 = TWOPI - SHIFT2
CASE ( 1)
SELECT CASE(INT(SIGN(1._EB,DV%ORIENTATION(2))))
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
LP%X = DV%X + PY%OFFSET*XTMP
LP%Y = DV%Y + PY%OFFSET*YTMP
LP%Z = DV%Z + PY%OFFSET*ZTMP
IF (TWO_D) THEN
LP%V = 0._EB
LP%Y = DV%Y
ENDIF
IF (LP%X<=XS .OR. LP%X>=XF) CYCLE CHOOSE_COORDS
IF (LP%Y<=YS .OR. LP%Y>=YF) CYCLE CHOOSE_COORDS
IF (LP%Z<=ZS .OR. LP%Z>=ZF) CYCLE CHOOSE_COORDS
CALL GET_IJK(LP%X,LP%Y,LP%Z,NM,XI,YJ,ZK,II,JJ,KK)
IC = CELL_INDEX(II,JJ,KK)
IF (.NOT.SOLID(IC)) EXIT CHOOSE_COORDS
ENDDO CHOOSE_COORDS
! Randomly choose PARTICLE size according to Cumulative Distribution Function (CDF)
IF (PC%MONODISPERSE) THEN
LP%R = 0.5_EB*PC%DIAMETER*D_PRES_FACTOR
LP%PWT = 1._EB
ELSE
CALL RANDOM_NUMBER(RN)
STRATUM = NINT(REAL(NSTRATA,EB)*REAL(RN,EB)+0.5_EB)
IL = PC%IL_CDF(STRATUM)
IU = PC%IU_CDF(STRATUM)
CALL RANDOM_CHOICE(PC%CDF(IL:IU), PC%R_CDF(IL:IU), IU-IL,LP%R)
LP%PWT = PC%W_CDF(STRATUM)
LP%R = D_PRES_FACTOR*LP%R
IF (2._EB*LP%R > PC%MAXIMUM_DIAMETER) THEN
LP%PWT = LP%PWT*LP%R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
LP%R = 0.5_EB*PC%MAXIMUM_DIAMETER
ENDIF
ENDIF
! Sum up mass of liquid being introduced
MASS_SUM = MASS_SUM + LP%PWT*PC%FTPR*LP%R**3
DROP_SUM = DROP_SUM + 1
ENDDO PARTICLE_INSERT_LOOP
! Compute weighting factor for the PARTICLEs just inserted
IF (DROP_SUM > 0) THEN
PWT0 = FLOW_RATE*DT/MASS_SUM
PARTICLE(NLP-DROP_SUM+1:NLP)%PWT = PARTICLE(NLP-DROP_SUM+1:NLP)%PWT*PWT0
ENDIF
! Indicate that PARTICLEs from this device have been inserted at this time T
DV%T = T
ENDDO SPRINKLER_INSERT_LOOP
END SUBROUTINE DEVICE_PARTICLE_INSERT
SUBROUTINE WALL_PARTICLE_INSERT
TYPE(WALL_TYPE), POINTER :: WC=>NULL()
! Loop through all boundary cells and insert particles if appropriate
WALL_INSERT_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
WC => WALL(IW)
IF (T < WC%TW) CYCLE WALL_INSERT_LOOP
IF (WALL(IW)%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE WALL_INSERT_LOOP
SF => SURFACE(WC%SURF_INDEX)
IPC = SF%PART_INDEX
IF (IPC < 1) CYCLE WALL_INSERT_LOOP
PC => PARTICLE_CLASS(IPC)
IF (PC%DEVC_INDEX>0) THEN
IF (.NOT.DEVICE(PC%DEVC_INDEX)%CURRENT_STATE) CYCLE WALL_INSERT_LOOP
ENDIF
IF (PC%CTRL_INDEX>0) THEN
IF (.NOT.CONTROL(PC%CTRL_INDEX)%CURRENT_STATE) CYCLE WALL_INSERT_LOOP
ENDIF
IF (T < SF%PARTICLE_INSERT_CLOCK(NM)) CYCLE WALL_INSERT_LOOP
INSERT_RATE = INSERT_RATE + WALL(IW)%NPPCW/SF%DT_INSERT
INSERT_COUNT = INSERT_COUNT + WALL(IW)%NPPCW
IF (WC%UW >= -0.0001_EB) CYCLE WALL_INSERT_LOOP
IF (WC%QCONF >= 0.0_EB) CYCLE WALL_INSERT_LOOP !for Level Set, particles only at fire front
II = WC%II
JJ = WC%JJ
KK = WC%KK
IC = CELL_INDEX(II,JJ,KK)
IF (.NOT.SOLID(IC)) CYCLE WALL_INSERT_LOOP
IF (NM > 1) THEN
IIG = WC%IIG
JJG = WC%JJG
KKG = WC%KKG
IF (INTERPOLATED_MESH(IIG,JJG,KKG) > 0) CYCLE WALL_INSERT_LOOP
ENDIF
! Loop over all particles for the IW-th cell
IOR = WC%IOR
MASS_SUM = 0._EB
PARTICLE_INSERT_LOOP2: DO I=1,WALL(IW)%NPPCW
IF (NLP+1 > MAXIMUM_PARTICLES) EXIT PARTICLE_INSERT_LOOP2
NLP = NLP+1
PARTICLE_TAG = PARTICLE_TAG + NMESHES
IF (NLP > NLPDIM) THEN
CALL RE_ALLOCATE_PARTICLES(1,NM,0,1000)
PARTICLE => MESHES(NM)%PARTICLE
ENDIF
LP=>PARTICLE(NLP)
IF (ABS(IOR)==1) THEN
IF (IOR== 1) LP%X = X(II) + VENT_OFFSET*DX(II+1)
IF (IOR==-1) LP%X = X(II-1) - VENT_OFFSET*DX(II-1)
CALL RANDOM_NUMBER(RN)
LP%Y = Y(JJ-1) + DY(JJ)*REAL(RN,EB)
CALL RANDOM_NUMBER(RN)
LP%Z = Z(KK-1) + DZ(KK)*REAL(RN,EB)
ENDIF
IF (ABS(IOR)==2) THEN
IF (IOR== 2) LP%Y = Y(JJ) + VENT_OFFSET*DY(JJ+1)
IF (IOR==-2) LP%Y = Y(JJ-1) - VENT_OFFSET*DY(JJ-1)
CALL RANDOM_NUMBER(RN)
LP%X = X(II-1) + DX(II)*REAL(RN,EB)
CALL RANDOM_NUMBER(RN)
LP%Z = Z(KK-1) + DZ(KK)*REAL(RN,EB)
ENDIF
IF (ABS(IOR)==3) THEN
IF (IOR== 3) LP%Z = Z(KK) + VENT_OFFSET*DZ(KK+1)
IF (IOR==-3) LP%Z = Z(KK-1) - VENT_OFFSET*DZ(KK-1)
CALL RANDOM_NUMBER(RN)
LP%X = X(II-1) + DX(II)*REAL(RN,EB)
CALL RANDOM_NUMBER(RN)
LP%Y = Y(JJ-1) + DY(JJ)*REAL(RN,EB)
ENDIF
SELECT CASE(IOR) ! Give particles an initial velocity (if desired)
CASE( 1)
LP%U = -WALL(IW)%UW
LP%V = SF%VEL_T(1)
LP%W = SF%VEL_T(2)
CASE(-1)
LP%U = WALL(IW)%UW
LP%V = SF%VEL_T(1)
LP%W = SF%VEL_T(2)
CASE( 2)
LP%U = SF%VEL_T(2)
LP%V = -WALL(IW)%UW
LP%W = SF%VEL_T(1)
CASE(-2)
LP%U = SF%VEL_T(2)
LP%V = WALL(IW)%UW
LP%W = SF%VEL_T(1)
CASE( 3)
LP%U = SF%VEL_T(1)
LP%V = SF%VEL_T(2)
LP%W = -WALL(IW)%UW
CASE(-3)
LP%U = SF%VEL_T(1)
LP%V = SF%VEL_T(2)
LP%W = WALL(IW)%UW
END SELECT
! Save the insertion time (TP) and scalar property (SP) for the particle
LP%ACCEL_X = 0._EB
LP%ACCEL_Y = 0._EB
LP%ACCEL_Z = 0._EB
LP%CLASS = IPC
LP%TAG = PARTICLE_TAG
IF (MOD(NLP,PC%SAMPLING)==0) THEN
LP%SHOW = .TRUE.
ELSE
LP%SHOW = .FALSE.
ENDIF
LP%SPLAT = .FALSE.
LP%WALL_INDEX = 0
LP%R = 0.0_EB
LP%TMP = PC%TMP_INITIAL
LP%T = T
LP%PWT = 1._EB
LP%IOR = 0
IF (PC%DIAMETER > 0._EB) THEN
CALL PARTICLE_SIZE_WEIGHT(LP%R,LP%PWT)
MASS_SUM = MASS_SUM + LP%PWT*PC%FTPR*LP%R**3
ENDIF
ENDDO PARTICLE_INSERT_LOOP2
IF (MASS_SUM > 0._EB) THEN
IF (SF%PARTICLE_MASS_FLUX > 0._EB) THEN
IF (ABS(WC%TW-T_BEGIN)<=ZERO_P .AND. SF%RAMP_INDEX(TIME_PART)>=1) THEN
TSI = T
ELSE
TSI = T - WC%TW
ENDIF
FLOW_RATE = EVALUATE_RAMP(TSI,SF%TAU(TIME_PART),SF%RAMP_INDEX(TIME_PART))*SF%PARTICLE_MASS_FLUX
PARTICLE(NLP-WALL(IW)%NPPCW+1:NLP)%PWT = &
PARTICLE(NLP-WALL(IW)%NPPCW+1:NLP)%PWT*FLOW_RATE*WALL(IW)%AREA_ADJUST*wALL(IW)%AW*SF%DT_INSERT/MASS_SUM
ENDIF
ENDIF
ENDDO WALL_INSERT_LOOP
END SUBROUTINE WALL_PARTICLE_INSERT
SUBROUTINE VOLUME_PARTICLE_INSERT
! Loop over all INIT lines and look for particles inserted within a specified volume
VOLUME_INSERT_LOOP: DO IB=1,N_INIT
IN => INITIALIZATION(IB)
! Determine if the INITIALIZATION type involves particles or PARTICLEs. If not, cycle.
IPC = IN%PART_INDEX
IF (IPC<1) CYCLE VOLUME_INSERT_LOOP
IF (IN%SINGLE_INSERTION .AND. IN%ALREADY_INSERTED) CYCLE VOLUME_INSERT_LOOP
! Determine if the particles/PARTICLEs are controlled by devices
PC => PARTICLE_CLASS(IPC)
IF (IN%DEVC_INDEX>0) THEN
IF (.NOT.DEVICE(IN%DEVC_INDEX)%CURRENT_STATE) THEN
IN%PARTICLE_INSERT_CLOCK(NM) = T
CYCLE VOLUME_INSERT_LOOP
ENDIF
ENDIF
IF (IN%CTRL_INDEX>0) THEN
IF (.NOT.CONTROL(IN%CTRL_INDEX)%CURRENT_STATE) THEN
IN%PARTICLE_INSERT_CLOCK(NM) = T
CYCLE VOLUME_INSERT_LOOP
ENDIF
ENDIF
! If it is not time to insert particles/PARTICLEs for this INITIALIZATION block, cycle.
IF (T < IN%PARTICLE_INSERT_CLOCK(NM)) CYCLE VOLUME_INSERT_LOOP
! Start processing the INITIALIZATION info
N_INITIAL = IN%NUMBER_INITIAL_PARTICLES
IF (N_INITIAL==0) CYCLE VOLUME_INSERT_LOOP
MASS_PER_VOLUME = IN%MASS_PER_VOLUME
MASS_PER_TIME = IN%MASS_PER_TIME
SELECT CASE(IN%SHAPE)
CASE('BLOCK')
IF (IN%X1>XF .OR. IN%X2<XS .OR. IN%Y1>YF .OR. IN%Y2<YS .OR. IN%Z1>ZF .OR. IN%Z2<ZS) CYCLE VOLUME_INSERT_LOOP
X1 = MAX(IN%X1,XS)
X2 = MIN(IN%X2,XF)
Y1 = MAX(IN%Y1,YS)
Y2 = MIN(IN%Y2,YF)
Z1 = MAX(IN%Z1,ZS)
Z2 = MIN(IN%Z2,ZF)
BLOCK_VOLUME = (X2-X1)*(Y2-Y1)*(Z2-Z1)
IF (BLOCK_VOLUME<=0._EB .AND. (MASS_PER_VOLUME>0._EB .OR. MASS_PER_TIME>0._EB)) CYCLE VOLUME_INSERT_LOOP
CASE('CONE')
CASE('POINT')
IF (IN%X0>XF .OR. IN%X0<XS .OR. IN%Y0>YF .OR. IN%Y0<YS .OR. IN%Z0>ZF .OR. IN%Z0<ZS) CYCLE VOLUME_INSERT_LOOP
END SELECT
! Assign properties to the initial PARTICLEs/particles
IN%ALREADY_INSERTED = .TRUE.
MASS_SUM = 0._EB
INSERT_PARTICLE_LOOP: DO I=1,N_INITIAL
NLP = NLP + 1
PARTICLE_TAG = PARTICLE_TAG + NMESHES
IF (NLP>NLPDIM) THEN
CALL RE_ALLOCATE_PARTICLES(1,NM,0,1000)
PARTICLE=>MESHES(NM)%PARTICLE
ENDIF
LP=>PARTICLE(NLP)
POINTWISE_IF: IF (IN%POINTWISE_PARTICLE_INIT) THEN
LP%X = IN%X0
LP%Y = IN%Y0
LP%Z = IN%Z0
LP%U = IN%U0
LP%V = IN%V0
LP%W = IN%W0
LP%R = IN%RADIUS
LP%TMP = IN%TEMPERATURE
IF (LP%R>0._EB) THEN
LP%PWT = IN%DT_INSERT*IN%MASS_PER_TIME/(PC%FTPR*LP%R**3)
ELSE
LP%PWT = 0._EB
ENDIF
ELSE POINTWISE_IF
BLOCK_OUT_LOOP: DO
!! CALL RANDOM_CONE(LP%X,LP%Y,LP%Z,0.5_EB*(X2+X1),0.5_EB*(Y2+Y1),Z1,0.25_EB,0.5_EB)
CALL RANDOM_RECTANGLE(LP%X,LP%Y,LP%Z,X1,X2,Y1,Y2,Z1,Z2)
CALL GET_IJK(LP%X,LP%Y,LP%Z,NM,XI,YJ,ZK,II,JJ,KK)
IF (.NOT.SOLID(CELL_INDEX(II,JJ,KK))) EXIT BLOCK_OUT_LOOP
ENDDO BLOCK_OUT_LOOP
LP%U = IN%U0 ! Initial velocity (default zero)
LP%V = IN%V0
LP%W = IN%W0
LP%R = 0._EB ! Radius is zero unless DIAMETER has been specified
LP%TMP = PC%TMP_INITIAL ! Initial temperature
LP%PWT = 1._EB ! Weighting factor is one unless changed below
ENDIF POINTWISE_IF
LP%T = T ! Insertion time is current time
LP%IOR = 0 ! Orientation of solid surface (0 means the PARTICLE/particle is not attached)
LP%CLASS = IPC ! Class identifier
LP%TAG = PARTICLE_TAG ! Unique integer tag
IF (MOD(NLP,PC%SAMPLING)==0) THEN ! Decide whether to show or not show in Smokeview
LP%SHOW = .TRUE.
ELSE
LP%SHOW = .FALSE.
ENDIF
LP%SPLAT = .FALSE.
LP%WALL_INDEX = 0
IF (PC%DIAMETER>0._EB .AND. .NOT.IN%POINTWISE_PARTICLE_INIT) CALL PARTICLE_SIZE_WEIGHT(LP%R,LP%PWT)
MASS_SUM = MASS_SUM + LP%PWT*PC%FTPR*LP%R**3 ! if r=0 the sum will stay 0
! Process special particles that are associated with a particular SURFace type
IF (PC%SURF_INDEX>0) THEN
SF => SURFACE(PC%SURF_INDEX)
LP%WALL_INDEX = IN%WALL_INDEX_START(NM) + I - 1
IF (SF%THICKNESS>0._EB) THEN
LP%R = SF%THICKNESS
ELSEIF (SF%RADIUS>0._EB) THEN
LP%R = SF%RADIUS
ENDIF
LP%PWT = 1._EB
WALL(LP%WALL_INDEX)%XW = LP%X
WALL(LP%WALL_INDEX)%YW = LP%Y
WALL(LP%WALL_INDEX)%ZW = LP%Z
WALL(LP%WALL_INDEX)%II = II
WALL(LP%WALL_INDEX)%JJ = JJ
WALL(LP%WALL_INDEX)%KK = KK
WALL(LP%WALL_INDEX)%IOR = 0
WALL(LP%WALL_INDEX)%SURF_INDEX = PC%SURF_INDEX
WALL(LP%WALL_INDEX)%II = II
WALL(LP%WALL_INDEX)%JJ = JJ
WALL(LP%WALL_INDEX)%KK = KK
ENDIF
IF (PC%N_SPLIT>1) THEN
LP%SPLIT_IOR = MOD(I-1,PC%N_SPLIT)+1
LP%PWT = LP%PWT/REAL(PC%N_SPLIT,EB)
ENDIF
ENDDO INSERT_PARTICLE_LOOP
POINTWISE_IF2: IF (.NOT.IN%POINTWISE_PARTICLE_INIT) THEN
! Adjust particle weighting factor PWT so that desired MASS_PER_VOLUME is achieved
IF (MASS_PER_TIME>0._EB) MASS_PER_VOLUME = MASS_PER_TIME*IN%DT_INSERT/BLOCK_VOLUME
IF (MASS_PER_VOLUME>0._EB) PARTICLE(NLP-N_INITIAL+1:NLP)%PWT = &
PARTICLE(NLP-N_INITIAL+1:NLP)%PWT*MASS_PER_VOLUME*BLOCK_VOLUME/MASS_SUM
ENDIF POINTWISE_IF2
ENDDO VOLUME_INSERT_LOOP
END SUBROUTINE VOLUME_PARTICLE_INSERT
SUBROUTINE PARTICLE_SIZE_WEIGHT(R,PWT)
REAL(EB), INTENT(OUT):: R,PWT
IF (PC%MONODISPERSE) THEN
R = 0.5_EB*PC%DIAMETER
PWT = 1._EB
ELSE
CALL RANDOM_NUMBER(RN)
STRATUM = NINT(REAL(NSTRATA,EB)*REAL(RN,EB)+0.5_EB)
IL = PC%IL_CDF(STRATUM)
IU = PC%IU_CDF(STRATUM)
CALL RANDOM_CHOICE(PC%CDF(IL:IU),PC%R_CDF(IL:IU),IU-IL,R)
PWT = PC%W_CDF(STRATUM)
IF (2._EB*R > PC%MAXIMUM_DIAMETER) THEN
PWT = PWT*R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
R = 0.5_EB*PC%MAXIMUM_DIAMETER
ENDIF
ENDIF
END SUBROUTINE PARTICLE_SIZE_WEIGHT
END SUBROUTINE INSERT_PARTICLES
SUBROUTINE UPDATE_PARTICLES(T,NM)
USE COMP_FUNCTIONS, ONLY : SECOND
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
REAL(EB) :: TNOW
! Return if this is an evacuation mesh
IF (EVACUATION_ONLY(NM)) RETURN
! Zero out the number of the PARTICLEs in the "orphanage"; that is, the place to hold PARTICLEs transferring from mesh to mesh
MESHES(NM)%OMESH(1:NMESHES)%N_DROP_ORPHANS = 0
! Return if there are no particles in this mesh
IF (MESHES(NM)%NLP==0) RETURN
! Set the CPU timer and point to the current mesh variables
TNOW=SECOND()
CALL POINT_TO_MESH(NM)
! Zero out the contribution by lagrangian particles to divergence
IF (N_EVAP_INDICES>0 .AND. .NOT.EVACUATION_ONLY(NM) .AND. CORRECTOR) THEN
D_LAGRANGIAN = 0._EB
ENDIF
! Move the PARTICLEs/particles, then compute mass and energy transfer, then add PARTICLE momentum to gas
IF (CORRECTOR) CALL MOVE_PARTICLES(T,NM)
IF (CORRECTOR) CALL PARTICLE_MASS_ENERGY_TRANSFER(T,NM)
CALL PARTICLE_MOMENTUM_TRANSFER(NM)
TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE UPDATE_PARTICLES
SUBROUTINE MOVE_PARTICLES(T,NM)
! Momentum transfer from all particles and PARTICLEs
USE COMP_FUNCTIONS, ONLY : SECOND
USE TRAN, ONLY: GET_IJK
REAL(EB), INTENT(IN) :: T
REAL :: RN
REAL(EB) :: SURFACE_PARTICLE_DIAMETER,P_UVWMAX,XI,YJ,ZK,RD,RDS,RDC,&
X_OLD,Y_OLD,Z_OLD,THETA_RN,STEP_FRACTION(-3:3)
LOGICAL :: HIT_SOLID
INTEGER :: ICN,I,IIN,JJN,KKN,II,JJ,KK,IW,IWP1,IWM1,IWP2,IWM2,IWP3,IWM3,IOR_OLD,IC,IOR_FIRST,IML
INTEGER, INTENT(IN) :: NM
TYPE (PARTICLE_TYPE), POINTER :: LP=>NULL()
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC=>NULL()
REAL(EB), POINTER, DIMENSION(:,:,:) :: NDPC=>NULL() ! number of PARTICLEs per cell
SURFACE_PARTICLE_DIAMETER = 0.001_EB ! All PARTICLEs adjusted to this size when on solid (m)
P_UVWMAX = 0._EB
! Sum up the number of PARTICLEs/particles in each grid cell (NDPC -- Number PARTICLEs Per Cell)
NDPC=>WORK1
NDPC=0._EB
DO I=1,NLP
LP => PARTICLE(I)
CALL GET_IJK(LP%X,LP%Y,LP%Z,NM,XI,YJ,ZK,II,JJ,KK)
IF (LP%PWT>0._EB .AND. LP%IOR==0) NDPC(II,JJ,KK)=NDPC(II,JJ,KK)+LP%PWT
ENDDO
! Loop through all Lagrangian particles and move them one time step
PARTICLE_LOOP: DO I=1,NLP
! Assign particle (LP%) and particle class (PC%) shortcuts
LP => PARTICLE(I)
PC => PARTICLE_CLASS(LP%CLASS)
! Determine particle radius
RD = LP%R
IF (.NOT. PC%MASSLESS .AND. RD<=0._EB) CYCLE PARTICLE_LOOP
RDS = RD*RD
RDC = RD*RDS
! Determine the current coordinates of the particle
CALL GET_IJK(LP%X,LP%Y,LP%Z,NM,XI,YJ,ZK,II,JJ,KK)
IC = CELL_INDEX(II,JJ,KK)
X_OLD = LP%X
Y_OLD = LP%Y
Z_OLD = LP%Z
! Throw out particles that are inside a solid obstruction
IF (SOLID(IC)) THEN
LP%X = 1.E6_EB
CYCLE PARTICLE_LOOP
ENDIF
SOLID_GAS_MOVE: IF (LP%IOR/=0) THEN
CALL MOVE_ON_SOLID(I)
ELSE SOLID_GAS_MOVE
CALL MOVE_IN_GAS(I,NM,P_UVWMAX,T)
! If the particle does not move, but does drag, go on to the next particle
IF (PC%MASSLESS .OR. LP%PWT<=ZERO_P .OR. PC%STATIC) CYCLE PARTICLE_LOOP
ENDIF SOLID_GAS_MOVE
! PARTICLE hits the floor
IF (POROUS_FLOOR .AND. LP%Z<ZS .AND. PC%EVAP_INDEX>0) THEN
IC = CELL_INDEX(II,JJ,1)
IW = WALL_INDEX(IC,-3)
IF (WALL(IW)%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. ACCUMULATE_WATER .AND. .NOT.LP%SPLAT) THEN
WALL(IW)%A_LP_MPUA(PC%EVAP_INDEX) = WALL(IW)%A_LP_MPUA(PC%EVAP_INDEX) + LP%PWT*PC%FTPR*RDC*WALL(IW)%RAW
LP%SPLAT = .TRUE.
ENDIF
CYCLE PARTICLE_LOOP
ENDIF
IF (.NOT.POROUS_FLOOR .AND. LP%Z<ZS) THEN
IC = CELL_INDEX(II,JJ,1)
IW = WALL_INDEX(IC,-3)
IF (WALL(IW)%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
IF (ACCUMULATE_WATER .AND. .NOT.LP%SPLAT) THEN
WALL(IW)%A_LP_MPUA(PC%EVAP_INDEX) = WALL(IW)%A_LP_MPUA(PC%EVAP_INDEX) + LP%PWT*PC%FTPR*RDC*WALL(IW)%RAW
LP%SPLAT = .TRUE.
ENDIF
LP%Z = ZS + 0.05_EB*DZ(1)
LP%IOR = 3
CALL RANDOM_NUMBER(RN)
THETA_RN = TWOPI*REAL(RN,EB)
LP%U = PC%HORIZONTAL_VELOCITY*COS(THETA_RN)
LP%V = PC%HORIZONTAL_VELOCITY*SIN(THETA_RN)
LP%W = 0._EB
ENDIF
ENDIF
! Where is the PARTICLE now? Limit the location by UBOUND and LBOUND due to the possible super fast PARTICLEs
IIN = MAX(LBOUND(CELLSI,1),MIN(UBOUND(CELLSI,1),FLOOR((LP%X-XS)*RDXINT)))
JJN = MAX(LBOUND(CELLSJ,1),MIN(UBOUND(CELLSJ,1),FLOOR((LP%Y-YS)*RDYINT)))
KKN = MAX(LBOUND(CELLSK,1),MIN(UBOUND(CELLSK,1),FLOOR((LP%Z-ZS)*RDZINT)))
XI = CELLSI(IIN)
YJ = CELLSJ(JJN)
ZK = CELLSK(KKN)
IIN = FLOOR(XI+1._EB)
JJN = FLOOR(YJ+1._EB)
KKN = FLOOR(ZK+1._EB)
IF (IIN<0 .OR. IIN>IBP1) CYCLE PARTICLE_LOOP
IF (JJN<0 .OR. JJN>JBP1) CYCLE PARTICLE_LOOP
IF (KKN<0 .OR. KKN>KBP1) CYCLE PARTICLE_LOOP
ICN = CELL_INDEX(IIN,JJN,KKN)
IF (IC==0 .OR. ICN==0) CYCLE PARTICLE_LOOP
IF (LP%X<XS .AND. WALL(WALL_INDEX(IC,-1))%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PARTICLE_LOOP
IF (LP%X>XF .AND. WALL(WALL_INDEX(IC, 1))%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PARTICLE_LOOP
IF (LP%Y<YS .AND. WALL(WALL_INDEX(IC,-2))%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PARTICLE_LOOP
IF (LP%Y>YF .AND. WALL(WALL_INDEX(IC, 2))%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PARTICLE_LOOP
IF (LP%Z<ZS .AND. WALL(WALL_INDEX(IC,-3))%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PARTICLE_LOOP
IF (LP%Z>ZF .AND. WALL(WALL_INDEX(IC, 3))%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PARTICLE_LOOP
! If PARTICLE hits an obstacle, change its properties
AIR_TO_SOLID: IF (II/=IIN .OR. JJ/=JJN .OR. KK/=KKN) THEN
IOR_OLD = LP%IOR
HIT_SOLID = .FALSE.
! Check if any solid boundaries of original grid cell have been crossed
IWP1 = WALL_INDEX(IC, 1)
IWM1 = WALL_INDEX(IC,-1)
IWP2 = WALL_INDEX(IC, 2)
IWM2 = WALL_INDEX(IC,-2)
IWP3 = WALL_INDEX(IC, 3)
IWM3 = WALL_INDEX(IC,-3)
STEP_FRACTION = 1._EB
IF (KKN>KK .AND. WALL(IWP3)%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
LP%IOR=-3
HIT_SOLID = .TRUE.
STEP_FRACTION(LP%IOR) = MAX(0._EB,(Z(KK)-Z_OLD-0.05_EB*DZ(KK))/(LP%Z-Z_OLD))
ENDIF
IF (KKN<KK .AND. WALL(IWM3)%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
LP%IOR= 3
HIT_SOLID = .TRUE.
STEP_FRACTION(LP%IOR) = MAX(0._EB,(Z(KK-1)-Z_OLD+0.05_EB*DZ(KK-1))/(LP%Z-Z_OLD))
ENDIF
IF (IIN>II .AND. WALL(IWP1)%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
LP%IOR=-1
HIT_SOLID = .TRUE.
STEP_FRACTION(LP%IOR) = MAX(0._EB,(X(II)-X_OLD-0.05_EB*DX(II))/(LP%X-X_OLD))
ENDIF
IF (IIN<II .AND. WALL(IWM1)%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
LP%IOR= 1
HIT_SOLID = .TRUE.
STEP_FRACTION(LP%IOR) = MAX(0._EB,(X(II-1)-X_OLD+0.05_EB*DX(II-1))/(LP%X-X_OLD))
ENDIF
IF (JJN>JJ .AND. WALL(IWP2)%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN