-
Notifications
You must be signed in to change notification settings - Fork 0
/
func.f90
2231 lines (1723 loc) · 56.3 KB
/
func.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 COMP_FUNCTIONS
! I/O + OS functions
USE PRECISION_PARAMETERS
IMPLICIT NONE
CHARACTER(255), PARAMETER :: funcid='$Id: func.f90 9953 2012-02-01 12:57:01Z Topi.Sikanen $'
CHARACTER(255), PARAMETER :: funcrev='$Revision: 9953 $'
CHARACTER(255), PARAMETER :: funcdate='$Date: 2012-02-01 04:57:01 -0800 (Wed, 01 Feb 2012) $'
CONTAINS
SUBROUTINE GET_REV_func(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE
WRITE(MODULE_DATE,'(A)') funcrev(INDEX(funcrev,':')+1:LEN_TRIM(funcrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') funcdate
END SUBROUTINE GET_REV_func
REAL(EB) FUNCTION SECOND() ! Returns the CPU time in seconds.
!$ USE OMP_LIB
REAL(FB) CPUTIME
CALL CPU_TIME(CPUTIME)
SECOND = CPUTIME
!$ SECOND = OMP_GET_WTIME()
END FUNCTION SECOND
REAL(EB) FUNCTION WALL_CLOCK_TIME() ! Returns the number of seconds since January 1, 2000, including leap years
! Thirty days hath September,
! April, June, and November.
! February has twenty-eight alone;
! All the rest have thirty-one,
! Excepting Leap-Year, that's the time
! When February's days are twenty-nine.
INTEGER :: DATE_TIME(8),WALL_CLOCK_SECONDS
CHARACTER(10) :: BIG_BEN(3)
! X_1 = common year, X_2 = leap year
INTEGER, PARAMETER :: S_PER_YEAR_1=31536000,S_PER_YEAR_2=31622400,S_PER_DAY=86400,S_PER_HOUR=3600,S_PER_MIN=60
INTEGER, PARAMETER, DIMENSION(12) :: ACCUMULATED_DAYS_1=(/0,31,59,90,120,151,181,212,243,273,304,334/), &
ACCUMULATED_DAYS_2=(/0,31,60,91,121,152,182,213,244,274,305,335/)
INTEGER :: YEAR_COUNT
CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
WALL_CLOCK_SECONDS = 0._EB
DO YEAR_COUNT=2001,DATE_TIME(1)
!Leap year if divisible by 4 but not 100 unless by 400 (1900 no, 1904 yes, 2000 yes)
IF (MOD(YEAR_COUNT,4)==0 .AND. (MOD(YEAR_COUNT,100)/=0 .OR. MOD(YEAR_COUNT,400)==0)) THEN
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_YEAR_2
ELSE
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_YEAR_1
ENDIF
ENDDO
IF (MOD(DATE_TIME(1),4)==0 .AND. (MOD(DATE_TIME(1),100)/=0 .OR. MOD(DATE_TIME(1),400)==0 )) THEN
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_DAY*(ACCUMULATED_DAYS_2(DATE_TIME(2))+DATE_TIME(3))
ELSE
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_DAY*(ACCUMULATED_DAYS_1(DATE_TIME(2))+DATE_TIME(3))
ENDIF
WALL_CLOCK_SECONDS = WALL_CLOCK_SECONDS + S_PER_HOUR*DATE_TIME(5) + S_PER_MIN*DATE_TIME(6) + DATE_TIME(7)
WALL_CLOCK_TIME = WALL_CLOCK_SECONDS + DATE_TIME(8)*0.001_EB
END FUNCTION WALL_CLOCK_TIME
SUBROUTINE GET_DATE(DATE)
INTEGER :: DATE_TIME(8)
CHARACTER(10) :: BIG_BEN(3),MONTH
CHARACTER(30), INTENT(OUT) :: DATE
CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)
SELECT CASE(DATE_TIME(2))
CASE(1)
MONTH='January'
CASE(2)
MONTH='February'
CASE(3)
MONTH='March'
CASE(4)
MONTH='April'
CASE(5)
MONTH='May'
CASE(6)
MONTH='June'
CASE(7)
MONTH='July'
CASE(8)
MONTH='August'
CASE(9)
MONTH='September'
CASE(10)
MONTH='October'
CASE(11)
MONTH='November'
CASE(12)
MONTH='December'
END SELECT
WRITE(DATE,'(A,I3,A,I4,2X,I2.2,A,I2.2,A,I2.2)') TRIM(MONTH),DATE_TIME(3),', ',DATE_TIME(1), &
DATE_TIME(5),':',DATE_TIME(6),':',DATE_TIME(7)
END SUBROUTINE GET_DATE
SUBROUTINE SHUTDOWN(MESSAGE)
! Stops the code gracefully after writing a message
USE GLOBAL_CONSTANTS, ONLY: LU_ERR,FN_OUTPUT,LU_OUTPUT,CHID
CHARACTER(*) MESSAGE
LOGICAL :: EX
WRITE(LU_ERR,'(/A)') TRIM(MESSAGE)
! For those running FDS on a PC, ensure that error message gets written to the .out file, too.
INQUIRE(FILE=FN_OUTPUT,EXIST=EX)
IF (.NOT.EX) OPEN(LU_OUTPUT,FILE=TRIM(CHID)//'.out',STATUS='REPLACE',FORM='FORMATTED')
WRITE(LU_OUTPUT,'(/A)') TRIM(MESSAGE)
STOP
END SUBROUTINE SHUTDOWN
SUBROUTINE GET_INPUT_FILE ! Read the argument after the command
USE GLOBAL_CONSTANTS, ONLY: FN_INPUT
IF (FN_INPUT=='null') CALL GETARG(1,FN_INPUT)
END SUBROUTINE GET_INPUT_FILE
INTEGER FUNCTION GET_FILE_NUMBER()
USE GLOBAL_CONSTANTS, ONLY: MYID,FILE_COUNTER
FILE_COUNTER(MYID) = FILE_COUNTER(MYID) + 1
GET_FILE_NUMBER = FILE_COUNTER(MYID)
END FUNCTION GET_FILE_NUMBER
SUBROUTINE CHECKREAD(NAME,LU,IOS)
INTEGER :: II
INTEGER, INTENT(OUT) :: IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(4), INTENT(IN) :: NAME
CHARACTER(80) TEXT
IOS = 1
READLOOP: DO
READ(LU,'(A)',END=10) TEXT
TLOOP: DO II=1,72
IF (TEXT(II:II)/='&' .AND. TEXT(II:II)/=' ') EXIT TLOOP
IF (TEXT(II:II)=='&') THEN
IF (TEXT(II+1:II+4)==NAME) THEN
BACKSPACE(LU)
IOS = 0
EXIT READLOOP
ELSE
CYCLE READLOOP
ENDIF
ENDIF
ENDDO TLOOP
ENDDO READLOOP
10 RETURN
END SUBROUTINE CHECKREAD
SUBROUTINE SEARCH_INPUT_FILE(NAME,LU,FOUND)
INTEGER :: II
LOGICAL, INTENT(OUT) :: FOUND
INTEGER, INTENT(IN) :: LU
CHARACTER(*), INTENT(IN) :: NAME
CHARACTER(255) TEXT
FOUND = .FALSE.
REWIND(LU)
READLOOP: DO
READ(LU,'(A)',END=10) TEXT
TLOOP: DO II=1,255-LEN_TRIM(NAME)
IF (TEXT(II:II+LEN_TRIM(NAME)-1)==TRIM(NAME)) THEN
FOUND = .TRUE.
EXIT READLOOP
ENDIF
ENDDO TLOOP
ENDDO READLOOP
10 REWIND(LU)
RETURN
END SUBROUTINE SEARCH_INPUT_FILE
SUBROUTINE CHECK_XB(XB)
REAL(EB) :: DUMMY,XB(6)
INTEGER :: I
DO I=1,5,2
IF (XB(I)>XB(I+1)) THEN
DUMMY = XB(I)
XB(I) = XB(I+1)
XB(I+1) = DUMMY
ENDIF
ENDDO
END SUBROUTINE CHECK_XB
SUBROUTINE CHANGE_UNITS(QUANTITY,UNITS,STATISTICS,MYID,LU_ERR)
CHARACTER(30), INTENT(IN) :: QUANTITY,STATISTICS
INTEGER, INTENT(IN) :: MYID,LU_ERR
CHARACTER(30), INTENT(INOUT) ::UNITS
CHARACTER(150) :: MESSAGE
INTEGER :: I,UNIT_L
UNIT_L = LEN(TRIM(UNITS))
SELECT CASE (STATISTICS)
CASE('VOLUME INTEGRAL')
I = INDEX(UNITS,'/m3')
IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
CASE('AREA INTEGRAL')
I = INDEX(UNITS,'/m2')
IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
CASE('SURFACE INTEGRAL')
I = INDEX(UNITS,'/m2')
IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
CASE('TIME INTEGRAL')
I = INDEX(UNITS,'/s')
IF (I/=0) THEN
WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+2:UNIT_L)
ELSE
I = INDEX(UNITS,'kW')
IF (I==1) WRITE(UNITS,'(A,A)') 'kJ',UNITS(3:UNIT_L)
IF (I>1) WRITE(UNITS,'(A,A,A)') UNITS(1:I-1),'kJ',UNITS(I+2:UNIT_L)
ENDIF
IF (QUANTITY=='FED') I=1
CASE DEFAULT
RETURN
END SELECT
IF (I==0) THEN
WRITE(MESSAGE,'(A,A,A,A)') 'WARNING: Problem with units compatibility using STATISTICS of ',TRIM(STATISTICS), &
' with the QUANTITY ',TRIM(QUANTITY)
IF (MYID==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
ENDIF
RETURN
END SUBROUTINE CHANGE_UNITS
END MODULE COMP_FUNCTIONS
MODULE MEMORY_FUNCTIONS
USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE COMP_FUNCTIONS, ONLY : SHUTDOWN
IMPLICIT NONE
CONTAINS
SUBROUTINE ChkMemErr(CodeSect,VarName,IZERO)
! Memory checking routine
CHARACTER(*), INTENT(IN) :: CodeSect, VarName
INTEGER IZERO
CHARACTER(100) MESSAGE
IF (IZERO==0) RETURN
WRITE(MESSAGE,'(4A)') 'ERROR: Memory allocation failed for ', TRIM(VarName),' in the routine ',TRIM(CodeSect)
CALL SHUTDOWN(MESSAGE)
END SUBROUTINE ChkMemErr
SUBROUTINE RE_ALLOCATE_PARTICLES(CODE,NM,NOM,NEW_DROPS)
TYPE (PARTICLE_TYPE), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: CODE,NM,NOM,NEW_DROPS
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE(OMESH_TYPE), POINTER :: M2=>NULL()
SELECT CASE(CODE)
CASE(1)
M=>MESHES(NM)
ALLOCATE(DUMMY(1:M%NLPDIM))
DUMMY = M%PARTICLE
DEALLOCATE(M%PARTICLE)
ALLOCATE(M%PARTICLE(M%NLPDIM+NEW_DROPS))
M%PARTICLE(1:M%NLPDIM) = DUMMY(1:M%NLPDIM)
M%NLPDIM = M%NLPDIM+NEW_DROPS
CASE(2)
M2=>MESHES(NM)%OMESH(NOM)
ALLOCATE(DUMMY(1:M2%N_DROP_ORPHANS_DIM))
DUMMY = M2%PARTICLE
DEALLOCATE(M2%PARTICLE)
ALLOCATE(M2%PARTICLE(M2%N_DROP_ORPHANS_DIM+NEW_DROPS))
M2%PARTICLE(1:M2%N_DROP_ORPHANS_DIM) = DUMMY(1:M2%N_DROP_ORPHANS_DIM)
M2%N_DROP_ORPHANS_DIM = M2%N_DROP_ORPHANS_DIM + NEW_DROPS
END SELECT
DEALLOCATE(DUMMY)
END SUBROUTINE RE_ALLOCATE_PARTICLES
FUNCTION REALLOCATE(P,N1,N2)
REAL(EB), POINTER, DIMENSION(:) :: P, REALLOCATE
INTEGER, INTENT(IN) :: N1,N2
INTEGER :: NOLD, IERR
CHARACTER(100) :: MESSAGE
ALLOCATE(REALLOCATE(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE'
CALL SHUTDOWN(MESSAGE)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE(N1:NOLD+N1-1) = P(N1:NOLD+N1-1) ! Restore the contents of the reallocated array
DEALLOCATE(P)
END FUNCTION REALLOCATE
SUBROUTINE RE_ALLOCATE_STRINGS(NM)
CHARACTER(80), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M=>NULL()
M=>MESHES(NM)
ALLOCATE(DUMMY(1:M%N_STRINGS))
DUMMY = M%STRING
DEALLOCATE(M%STRING)
ALLOCATE(M%STRING(M%N_STRINGS_MAX+100))
M%STRING(1:M%N_STRINGS) = DUMMY(1:M%N_STRINGS)
M%N_STRINGS_MAX = M%N_STRINGS_MAX+100
DEALLOCATE(DUMMY)
END SUBROUTINE RE_ALLOCATE_STRINGS
SUBROUTINE ALLOCATE_WALL_ARRAYS(NM,IW)
USE GLOBAL_CONSTANTS, ONLY: RADIATION,SIGMA,TMPA4,OPEN_BOUNDARY,VIRTUAL_BOUNDARY,EVACUATION_ONLY,&
NRA=>NUMBER_RADIATION_ANGLES,NSB=>NUMBER_SPECTRAL_BANDS, &
TMPA,ANGLE_INCREMENT
INTEGER, INTENT(IN) :: NM,IW
INTEGER :: IZERO,II,IL,NN,N
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE (SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (WALL_TYPE), POINTER :: WC=>NULL()
M=>MESHES(NM)
WC => M%WALL(IW)
SF => SURFACE(WC%SURF_INDEX)
IF (EVACUATION_ONLY(NM)) RETURN
IF (RADIATION) THEN
IF (WC%BOUNDARY_TYPE==OPEN_BOUNDARY) THEN
ALLOCATE(WC%ILW(ANGLE_INCREMENT,NSB),STAT=IZERO)
CALL ChkMemErr('FUNC','ILW',IZERO)
ELSE
ALLOCATE(WC%ILW(NRA,NSB),STAT=IZERO)
CALL ChkMemErr('FUNC','ILW',IZERO)
ENDIF
IF (WC%BOUNDARY_TYPE==VIRTUAL_BOUNDARY) THEN
WC%ILW = 0._EB
ELSE
WC%ILW = SIGMA*TMPA4*RPI
ENDIF
WC%QRADIN = WC%E_WALL*SIGMA*TMPA4
WC%QRADOUT = WC%E_WALL*SIGMA*TMPA4
ENDIF
IF (WC%BOUNDARY_TYPE==OPEN_BOUNDARY) RETURN
IF (SF%THERMALLY_THICK) THEN
ALLOCATE(WC%TMP_S(0:SF%N_CELLS+1),STAT=IZERO)
CALL ChkMemErr('FUNC','TMP_S',IZERO)
ALLOCATE(WC%RHO_S(0:SF%N_CELLS+1,SF%N_MATL),STAT=IZERO)
CALL ChkMemErr('FUNC','RHO_S',IZERO)
WC%RHO_S = 0._EB
DO II=0,SF%N_CELLS+1
IL = SF%LAYER_INDEX(II)
IF (SF%TMP_INNER(IL)>0._EB) THEN
WC%TMP_S(II) = SF%TMP_INNER(IL)
ELSE
WC%TMP_S(II) = M%TMP_0(WC%KK)
ENDIF
DO NN=1,SF%N_LAYER_MATL(IL)
DO N=1,SF%N_MATL
IF (SF%LAYER_MATL_INDEX(IL,NN)==SF%MATL_INDEX(N)) &
WC%RHO_S(II,N) = SF%LAYER_MATL_FRAC(IL,NN)*SF%LAYER_DENSITY(IL)
ENDDO
ENDDO
ENDDO
IF (SF%SHRINK) THEN
ALLOCATE(WC%LAYER_THICKNESS(SF%N_LAYERS),STAT=IZERO)
CALL ChkMemErr('FUNC','LAYER_THICKNESS',IZERO)
WC%LAYER_THICKNESS = SF%LAYER_THICKNESS(1:SF%N_LAYERS)
ALLOCATE(WC%N_LAYER_CELLS(SF%N_LAYERS),STAT=IZERO)
CALL ChkMemErr('FUNC','N_LAYER_CELLS',IZERO)
WC%N_LAYER_CELLS = SF%N_LAYER_CELLS
ALLOCATE(WC%X_S(0:SF%N_CELLS),STAT=IZERO)
CALL ChkMemErr('FUNC','X_S',IZERO)
WC%X_S = SF%X_S
ENDIF
WC%SHRINKING = .FALSE.
WC%BURNAWAY = .FALSE.
WC%TMP_F = WC%TMP_S(0)
WC%TMP_B = WC%TMP_S(SF%N_CELLS+1)
ENDIF
IF (SF%VEGETATION) THEN
ALLOCATE(WC%VEG_ASHMASS_L(SF%NVEG_L),STAT=IZERO)
CALL ChkMemErr('FUNC','VEG_ASHMASS_L',IZERO)
ALLOCATE(WC%VEG_CHARMASS_L(SF%NVEG_L),STAT=IZERO)
CALL ChkMemErr('FUNC','VEG_CHARMASS_L',IZERO)
ALLOCATE(WC%VEG_FUELMASS_L(SF%NVEG_L),STAT=IZERO)
CALL ChkMemErr('FUNC','VEG_FUELMASS_L',IZERO)
ALLOCATE(WC%VEG_MOISTMASS_L(SF%NVEG_L),STAT=IZERO)
CALL ChkMemErr('FUNC','VEG_MOISTMASS_L',IZERO)
ALLOCATE(WC%VEG_TMP_L(0:SF%NVEG_L),STAT=IZERO)
CALL ChkMemErr('FUNC','VEG_TMP_L',IZERO)
WC%VEG_ASHMASS_L(:) = 0.0_EB
WC%VEG_CHARMASS_L(:) = 0.0_EB
WC%VEG_FUELMASS_L(:) = SF%VEG_LOAD / SF%NVEG_L
WC%VEG_MOISTMASS_L(:) = SF%VEG_MOISTURE*WC%VEG_FUELMASS_L(:)
WC%VEG_TMP_L(:) = SF%VEG_INITIAL_TEMP + 273.15_EB
WC%TMP_F = SF%VEG_INITIAL_TEMP + 273.15_EB
IF(SF%VEG_NO_BURN) WC%TMP_F = TMPA
ENDIF
WC%ALREADY_ALLOCATED = .TRUE.
END SUBROUTINE ALLOCATE_WALL_ARRAYS
END MODULE MEMORY_FUNCTIONS
MODULE GEOMETRY_FUNCTIONS
! Functions for manipulating geometry
USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
IMPLICIT NONE
CONTAINS
SUBROUTINE SEARCH_OTHER_MESHES(XX,YY,ZZ,NOM,IIO,JJO,KKO)
! Given the point (XX,YY,ZZ), determine which other mesh it intersects and what its indices are.
REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(OUT) :: NOM,IIO,JJO,KKO
TYPE (MESH_TYPE), POINTER :: M2=>NULL()
OTHER_MESH_LOOP: DO NOM=1,NMESHES
IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP
M2=>MESHES(NOM)
IF (XX>=M2%XS .AND. XX<=M2%XF .AND. YY>=M2%YS .AND. YY<=M2%YF .AND. ZZ>=M2%ZS .AND. ZZ<=M2%ZF) THEN
XI = MAX( 1._EB , MIN( REAL(M2%IBAR,EB)+ALMOST_ONE , M2%CELLSI(FLOOR((XX-M2%XS)*M2%RDXINT)) + 1._EB ) )
YJ = MAX( 1._EB , MIN( REAL(M2%JBAR,EB)+ALMOST_ONE , M2%CELLSJ(FLOOR((YY-M2%YS)*M2%RDYINT)) + 1._EB ) )
ZK = MAX( 1._EB , MIN( REAL(M2%KBAR,EB)+ALMOST_ONE , M2%CELLSK(FLOOR((ZZ-M2%ZS)*M2%RDZINT)) + 1._EB ) )
IIO = FLOOR(XI)
JJO = FLOOR(YJ)
KKO = FLOOR(ZK)
RETURN
ENDIF
ENDDO OTHER_MESH_LOOP
NOM = 0
END SUBROUTINE SEARCH_OTHER_MESHES
SUBROUTINE BLOCK_CELL(NM,I1,I2,J1,J2,K1,K2,IVAL,OBST_INDEX)
! Indicate which cells are blocked off
INTEGER :: NM,I1,I2,J1,J2,K1,K2,IVAL,I,J,K,OBST_INDEX,IC
TYPE (MESH_TYPE), POINTER :: M=>NULL()
M => MESHES(NM)
DO K=K1,K2
DO J=J1,J2
DO I=I1,I2
IC = M%CELL_INDEX(I,J,K)
SELECT CASE(IVAL)
CASE(0)
M%SOLID(IC) = .FALSE.
M%OBST_INDEX_C(IC) = 0
CASE(1)
M%SOLID(IC) = .TRUE.
M%OBST_INDEX_C(IC) = OBST_INDEX
END SELECT
IF (OBST_INDEX==0) M%EXTERIOR(IC) = .TRUE.
ENDDO
ENDDO
ENDDO
END SUBROUTINE BLOCK_CELL
LOGICAL FUNCTION INTERIOR(XX,YY,ZZ)
INTEGER NM
REAL(EB), INTENT(IN) :: XX,YY,ZZ
INTERIOR = .FALSE.
DO NM=1,NMESHES
IF (XX>MESHES(NM)%XS .AND. XX<MESHES(NM)%XF .AND. &
YY>MESHES(NM)%YS .AND. YY<MESHES(NM)%YF .AND. &
ZZ>MESHES(NM)%ZS .AND. ZZ<MESHES(NM)%ZF) INTERIOR = .TRUE.
ENDDO
END FUNCTION INTERIOR
SUBROUTINE ASSIGN_PRESSURE_ZONE(NM,XX,YY,ZZ,I_ZONE)
! Given the point (XX,YY,ZZ) within Mesh NM, determine all mesh cells within the same pressure zone
REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(IN) :: NM,I_ZONE
INTEGER :: NN,IOR,IC,II,JJ,KK,III,JJJ,KKK,Q_N,IIO,JJO,KKO,NOM
INTEGER, ALLOCATABLE, DIMENSION(:) :: Q_I,Q_J,Q_K
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE (OBSTRUCTION_TYPE), POINTER :: OB=>NULL()
IF (EVACUATION_ONLY(NM)) RETURN
M=>MESHES(NM)
ALLOCATE(Q_I(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_J(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_K(M%IBAR*M%JBAR*M%KBAR))
! Find the cell indices corresponding to the given point
XI = MAX( 1._EB , MIN( REAL(M%IBAR,EB)+ALMOST_ONE , M%CELLSI(NINT((XX-M%XS)*M%RDXINT)) + 1._EB ) )
YJ = MAX( 1._EB , MIN( REAL(M%JBAR,EB)+ALMOST_ONE , M%CELLSJ(NINT((YY-M%YS)*M%RDYINT)) + 1._EB ) )
ZK = MAX( 1._EB , MIN( REAL(M%KBAR,EB)+ALMOST_ONE , M%CELLSK(NINT((ZZ-M%ZS)*M%RDZINT)) + 1._EB ) )
II = FLOOR(XI)
JJ = FLOOR(YJ)
KK = FLOOR(ZK)
! Add the first entry to "queue" of cells that need a pressure zone number
Q_I(1) = II
Q_J(1) = JJ
Q_K(1) = KK
Q_N = 1
M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
! Look to all cells adjacent to the starting cell and determine if they are in the ZONE as well.
! Repeat process until all cells are found.
SORT_QUEUE: DO
IF (Q_N<1) EXIT SORT_QUEUE
III = Q_I(Q_N)
JJJ = Q_J(Q_N)
KKK = Q_K(Q_N)
IC = M%CELL_INDEX(III,JJJ,KKK)
Q_N = Q_N - 1
SEARCH_LOOP: DO IOR=-3,3
IF (IOR==0) CYCLE SEARCH_LOOP
SELECT CASE(IOR)
CASE(-1)
IF (III==0) CYCLE SEARCH_LOOP
II = III-1
JJ = JJJ
KK = KKK
CASE( 1)
IF (III==M%IBP1) CYCLE SEARCH_LOOP
II = III+1
JJ = JJJ
KK = KKK
CASE(-2)
IF (JJJ==0) CYCLE SEARCH_LOOP
II = III
JJ = JJJ-1
KK = KKK
CASE( 2)
IF (JJJ==M%JBP1) CYCLE SEARCH_LOOP
II = III
JJ = JJJ+1
KK = KKK
CASE(-3)
IF (KKK==0) CYCLE SEARCH_LOOP
II = III
JJ = JJJ
KK = KKK-1
CASE( 3)
IF (KKK==M%KBP1) CYCLE SEARCH_LOOP
II = III
JJ = JJJ
KK = KKK+1
END SELECT
! If the cell is outside the computational domain, check if it is in another mesh
IF (II<1 .OR. II>M%IBAR .OR. JJ<1 .OR. JJ>M%JBAR .OR. KK<1 .OR. KK>M%KBAR) THEN
CALL SEARCH_OTHER_MESHES(M%XC(II),M%YC(JJ),M%ZC(KK),NOM,IIO,JJO,KKO)
IF (NOM>0) M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
CYCLE SEARCH_LOOP
ENDIF
! Look for thin obstructions bordering the current cell
DO NN=1,M%N_OBST
OB=>M%OBSTRUCTION(NN)
SELECT CASE(IOR)
CASE(-1)
IF (II== OB%I1 .AND. II== OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
CASE( 1)
IF (II-1==OB%I1 .AND. II-1==OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
CASE(-2)
IF (JJ== OB%J1 .AND. JJ== OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
CASE( 2)
IF (JJ-1==OB%J1 .AND. JJ-1==OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
CASE(-3)
IF (KK== OB%K1 .AND. KK== OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
CASE( 3)
IF (KK-1==OB%K1 .AND. KK-1==OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
END SELECT
ENDDO
! If an obstruction is found, assign its cells the current ZONE, just in case the obstruction is removed
IC = M%CELL_INDEX(II,JJ,KK)
IF (M%SOLID(IC) .AND. M%OBST_INDEX_C(IC)>0) THEN
OB => M%OBSTRUCTION(M%OBST_INDEX_C(IC))
M%PRESSURE_ZONE(OB%I1+1:OB%I2,OB%J1+1:OB%J2,OB%K1+1:OB%K2) = I_ZONE
CYCLE SEARCH_LOOP
ENDIF
! If the neighboring cell is not solid, assign the pressure zone
IF (.NOT.M%SOLID(IC) .AND. M%PRESSURE_ZONE(II,JJ,KK)<1) THEN
Q_N = Q_N+1
Q_I(Q_N) = II
Q_J(Q_N) = JJ
Q_K(Q_N) = KK
M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
ENDIF
ENDDO SEARCH_LOOP
END DO SORT_QUEUE
DEALLOCATE(Q_I)
DEALLOCATE(Q_J)
DEALLOCATE(Q_K)
END SUBROUTINE ASSIGN_PRESSURE_ZONE
SUBROUTINE GET_N_LAYER_CELLS(DIFFUSIVITY,THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR,N_CELLS,DX_MIN)
! Get number of wall cells in the layer
INTEGER, INTENT(OUT) :: N_CELLS
REAL(EB), INTENT(OUT) :: DX_MIN
REAL(EB), INTENT(IN) :: DIFFUSIVITY,THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR
REAL(EB) :: DDSUM
INTEGER :: N, I
DX_MIN = 0._EB
IF (THICKNESS<=ZERO_P) THEN
N_CELLS = 0
DX_MIN = 0._EB
RETURN
ENDIF
SHRINK_LOOP: DO N=1,999
DDSUM = 0._EB
SUM_LOOP: DO I=1,N
DDSUM = DDSUM + STRETCH_FACTOR**(MIN(I-1,N-I))
ENDDO SUM_LOOP
IF ((THICKNESS/DDSUM < CELL_SIZE_FACTOR*SQRT(DIFFUSIVITY)) .OR. (N==999)) THEN
N_CELLS = N
DX_MIN = THICKNESS/DDSUM
EXIT SHRINK_LOOP
ENDIF
ENDDO SHRINK_LOOP
END SUBROUTINE GET_N_LAYER_CELLS
SUBROUTINE GET_WALL_NODE_COORDINATES(N_CELLS,N_LAYERS,N_LAYER_CELLS,SMALLEST_CELL_SIZE,STRETCH_FACTOR,X_S)
! Get the wall internal coordinates
INTEGER, INTENT(IN) :: N_CELLS,N_LAYERS, N_LAYER_CELLS(N_LAYERS)
REAL(EB), INTENT(IN) :: SMALLEST_CELL_SIZE(N_LAYERS),STRETCH_FACTOR(N_LAYERS)
REAL(EB), INTENT(OUT) :: X_S(0:N_CELLS)
INTEGER I, II, NL
REAL(EB) DX_S
II = 0
X_S(0) = 0._EB
DO NL=1,N_LAYERS
DO I=1,N_LAYER_CELLS(NL)
II = II+1
DX_S = SMALLEST_CELL_SIZE(NL)*STRETCH_FACTOR(NL)**(MIN(I-1,N_LAYER_CELLS(NL)-I))
X_S(II) = X_S(II-1) + DX_S
ENDDO
ENDDO
END SUBROUTINE GET_WALL_NODE_COORDINATES
SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
LAYER_THICKNESS,GEOMETRY,X_S,LAYER_DIVIDE,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX,MF_FRAC)
! Get the wall internal coordinates
INTEGER, INTENT(IN) :: N_CELLS, N_LAYERS, N_LAYER_CELLS(N_LAYERS), GEOMETRY
REAL(EB), INTENT(IN) :: X_S(0:N_CELLS),LAYER_THICKNESS(1:N_LAYERS), LAYER_DIVIDE
INTEGER, INTENT(OUT) :: LAYER_INDEX(0:N_CELLS+1)
REAL(EB), INTENT(OUT) :: DX(1:N_CELLS),RDX(0:N_CELLS+1),RDXN(0:N_CELLS),DX_WGT(0:N_CELLS),DXF,DXB, &
MF_FRAC(1:N_CELLS)
INTEGER I, II, NL, I_GRAD
REAL(EB) R, THICKNESS, X_DIVIDE
THICKNESS = SUM(LAYER_THICKNESS)
SELECT CASE(GEOMETRY)
CASE(SURF_CARTESIAN)
I_GRAD = 0
CASE(SURF_CYLINDRICAL)
I_GRAD = 1
CASE(SURF_SPHERICAL)
I_GRAD = 2
END SELECT
II = 0
DO NL=1,N_LAYERS
DO I=1,N_LAYER_CELLS(NL)
II = II + 1
LAYER_INDEX(II) = NL
ENDDO
ENDDO
LAYER_INDEX(0) = 1
LAYER_INDEX(N_CELLS+1) = N_LAYERS
DXF = X_S(1) - X_S(0)
DXB = X_S(N_CELLS) - X_S(N_CELLS-1)
! Compute dx_weight for each node (dx_weight is the ratio of cell size to the combined size of the current and next cell)
DO I=1,N_CELLS-1
DX_WGT(I) = (X_S(I)-X_S(I-1))/(X_S(I+1)-X_S(I-1))
ENDDO
DX_WGT(0) = 0.5_EB
DX_WGT(N_CELLS) = 0.5_EB
! Compute dx and 1/dx for each node (dx is the distance from cell boundary to cell boundary)
DO I=1,N_CELLS
DX(I) = X_S(I)-X_S(I-1)
RDX(I) = 1._EB/DX(I)
ENDDO
! Adjust 1/dx_n to 1/rdr for cylindrical case and 1/r^2dr for spaherical
IF (GEOMETRY /=SURF_CARTESIAN) THEN
DO I=1,N_CELLS
R = THICKNESS-0.5_EB*(X_S(I)+X_S(I-1))
RDX(I) = RDX(I)/R**I_GRAD
ENDDO
ENDIF
RDX(0) = RDX(1)
RDX(N_CELLS+1) = RDX(N_CELLS)
! Compute 1/dx_n for each node (dx_n is the distance from cell center to cell center)
DO I=1,N_CELLS-1
RDXN(I) = 2._EB/(X_S(I+1)-X_S(I-1))
ENDDO
RDXN(0) = 1._EB/(X_S(1)-X_S(0))
RDXN(N_CELLS) = 1._EB/(X_S(N_CELLS)-X_S(N_CELLS-1))
! Adjust 1/dx_n to r/dr for cylindrical case and r^2/dr for spaherical
IF (GEOMETRY /= SURF_CARTESIAN) THEN
DO I=0,N_CELLS
R = THICKNESS-X_S(I)
RDXN(I) = RDXN(I)*R**I_GRAD
ENDDO
ENDIF
! Compute mass flux fraction array (array numbers indicate the fraction of mass flux that is added to the front
IF (LAYER_DIVIDE >= REAL(N_LAYERS,EB)) THEN
MF_FRAC = 1.0_EB
ELSE
MF_FRAC = 0._EB
X_DIVIDE = 0._EB
DO NL=1,N_LAYERS
IF (LAYER_DIVIDE>=REAL(NL,EB)) THEN
X_DIVIDE = X_DIVIDE + LAYER_THICKNESS(NL)
ELSE
X_DIVIDE = X_DIVIDE + MOD(LAYER_DIVIDE,1.0_EB)*LAYER_THICKNESS(NL)
EXIT
ENDIF
ENDDO
II = 0
DIVILOOP: DO NL=1,N_LAYERS
DO I=1,N_LAYER_CELLS(NL)
II = II + 1
IF (X_S(II) < X_DIVIDE) THEN
MF_FRAC(II) = 1._EB
ELSEIF (X_S(II-1) < X_DIVIDE) THEN
MF_FRAC(II) = (X_DIVIDE-X_S(II-1))/DX(II)
EXIT DIVILOOP
ENDIF
ENDDO
ENDDO DIVILOOP
ENDIF
END SUBROUTINE GET_WALL_NODE_WEIGHTS
SUBROUTINE GET_INTERPOLATION_WEIGHTS(N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS,N_LAYER_CELLS_NEW,X_S,X_S_NEW,INT_WGT)
INTEGER, INTENT(IN) :: N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS(N_LAYERS),N_LAYER_CELLS_NEW(N_LAYERS)
REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW)
REAL(EB), INTENT(OUT) :: INT_WGT(:,:)
REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,DX_NEW
INTEGER I, J, II, JJ, I_BASE, J_BASE, J_OLD,N
II = 0
JJ = 0
I_BASE = 0
J_BASE = 0
INT_WGT = 0._EB
DO N = 1,N_LAYERS
J_OLD = 1
DO I = 1,N_LAYER_CELLS_NEW(N)
II = I_BASE + I
XUP_NEW = X_S_NEW(II)
XLOW_NEW = X_S_NEW(II-1)
DX_NEW = XUP_NEW - XLOW_NEW
DO J = J_OLD,N_LAYER_CELLS(N)
JJ = J_BASE + J
XUP = X_S(JJ)
XLOW = X_S(JJ-1)
INT_WGT(II,JJ) = (MIN(XUP,XUP_NEW)-MAX(XLOW,XLOW_NEW))/DX_NEW
IF (XUP >= XUP_NEW) EXIT
ENDDO
J_OLD = J
ENDDO
I_BASE = I_BASE + N_LAYER_CELLS_NEW(N)
J_BASE = J_BASE + N_LAYER_CELLS(N)
ENDDO
END SUBROUTINE GET_INTERPOLATION_WEIGHTS
SUBROUTINE INTERPOLATE_WALL_ARRAY(NWP,NWP_NEW,INT_WGT,ARR)
INTEGER, INTENT(IN) :: NWP,NWP_NEW
REAL(EB), INTENT(IN) :: INT_WGT(:,:)
REAL(EB) ARR(NWP),TMP(NWP)
INTEGER I,J
TMP = ARR
ARR = 0._EB
DO I = 1,NWP_NEW
DO J = 1,NWP
ARR(I) = ARR(I) + INT_WGT(I,J)*TMP(J)
ENDDO
ENDDO
END SUBROUTINE INTERPOLATE_WALL_ARRAY
SUBROUTINE RANDOM_RECTANGLE(XX,YY,ZZ,X1,X2,Y1,Y2,Z1,Z2)
! Choose a random point (XX,YY,ZZ) from within a rectangular volume bounded by X1, X2, ...
REAL(EB), INTENT(IN) :: X1,X2,Y1,Y2,Z1,Z2
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL :: RN
CALL RANDOM_NUMBER(RN)
XX = X1 + REAL(RN,EB)*(X2-X1)
CALL RANDOM_NUMBER(RN)
YY = Y1 + REAL(RN,EB)*(Y2-Y1)
CALL RANDOM_NUMBER(RN)
ZZ = Z1 + REAL(RN,EB)*(Z2-Z1)
END SUBROUTINE RANDOM_RECTANGLE
SUBROUTINE RANDOM_CONE(XX,YY,ZZ,X0,Y0,Z0,DD,HH)
! Choose a random point (XX,YY,ZZ) from within a vertically oriented cone
REAL(EB), INTENT(IN) :: X0,Y0,Z0,DD,HH
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: THETA,RR
REAL :: RN
CALL RANDOM_NUMBER(RN)
ZZ = Z0 + HH*REAL(RN,EB)**2
CALL RANDOM_NUMBER(RN)
THETA = TWOPI*REAL(RN,EB)
CALL RANDOM_NUMBER(RN)
RR = 0.5_EB*DD*REAL(RN,EB)*(1._EB-(ZZ-Z0)/HH)
XX = X0 + RR*COS(THETA)
YY = Y0 + RR*SIN(THETA)
END SUBROUTINE RANDOM_CONE
END MODULE GEOMETRY_FUNCTIONS
MODULE PHYSICAL_FUNCTIONS
! Functions for physical quantities
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
IMPLICIT NONE
CONTAINS
SUBROUTINE GET_MASS_FRACTION(Z_IN,INDEX,Y_OUT)
! Y_OUT returns the mass fraction of species INDEX
INTEGER, INTENT(IN) :: INDEX
REAL(EB) :: Z_IN(0:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT
IF (N_TRACKED_SPECIES>0) THEN
Z_IN(0) = 1._EB - SUM(Z_IN(1:N_TRACKED_SPECIES))
ELSE
Z_IN(0) = 1._EB
ENDIF
Y_OUT = DOT_PRODUCT(Z2Y(INDEX,:),Z_IN)
Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))
END SUBROUTINE GET_MASS_FRACTION