@@ -12,6 +12,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
12
12
coef_bb_dc , fire_hist , hwp , hwp_prevd , &
13
13
swdown ,ebb_dcycle , ebu_in , ebu ,fire_type ,&
14
14
q_vap , add_fire_moist_flux , &
15
+ sc_factor , &
15
16
ids ,ide , jds ,jde , kds ,kde , &
16
17
ims ,ime , jms ,jme , kms ,kme , &
17
18
its ,ite , jts ,jte , kts ,kte ,mpiid )
@@ -34,7 +35,6 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
34
35
real (kind_phys), DIMENSION (ims:ime,jms:jme), INTENT (IN ) :: hwp, peak_hr, fire_end_hr, ebu_in ! RAR: Shall we make fire_end integer?
35
36
real (kind_phys), DIMENSION (ims:ime,jms:jme), INTENT (INOUT ) :: coef_bb_dc ! RAR:
36
37
real (kind_phys), DIMENSION (ims:ime,jms:jme), INTENT (IN ) :: hwp_prevd
37
-
38
38
real (kind_phys), DIMENSION (ims:ime,kms:kme,jms:jme), INTENT (IN ) :: dz8w,rho_phy ! ,rel_hum
39
39
real (kind_phys), INTENT (IN ) :: dtstep, gmt
40
40
real (kind_phys), INTENT (IN ) :: time_int, pi, ebb_min ! RAR: time in seconds since start of simulation
@@ -55,12 +55,17 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
55
55
real (kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm ! For BB emis. diurnal cycle calculation
56
56
57
57
! For Gaussian diurnal cycle
58
- real (kind_phys), PARAMETER :: sc_factor= 1 . ! to scale up the wildfire emissions, TBD later
58
+ real (kind_phys), INTENT ( IN ) :: sc_factor ! to scale up the wildfire emissions
59
59
real (kind_phys), PARAMETER :: rinti= 2.1813936e-8 , ax2= 3400 ., const2= 130 ., &
60
60
coef2= 10.6712963e-4 , cx2= 7200 ., timeq_max= 3600 .* 24 .
61
61
! >-- Fire parameters: Fores west, Forest east, Shrubland, Savannas, Grassland, Cropland
62
62
real (kind_phys), dimension (1 :5 ), parameter :: avg_fire_dur = (/ 8.9 , 4.2 , 3.3 , 3.0 , 1.4 / )
63
63
real (kind_phys), dimension (1 :5 ), parameter :: sigma_fire_dur = (/ 8.7 , 6.0 , 5.5 , 5.2 , 2.4 / )
64
+ ! For fire diurnal cycle calculation
65
+ ! real(kind_phys), parameter :: avgx1=-2.0, sigmx1=0.7, C1=0.083 ! Ag fires
66
+ ! real(kind_phys), parameter :: avgx2=-0.1, sigmx2=0.8, C2=0.55 ! Grass fires, slash burns
67
+ real (kind_phys), parameter :: avgx1= 0 ., sigmx1= 2.2 , C1= 0.2 ! Ag fires
68
+ real (kind_phys), parameter :: avgx2= 0.5 , sigmx2= 0.8 , C2= 1.1 ! Grass fires, slash burns
64
69
65
70
timeq= gmt* 3600._kind_phys + real (time_int,4 )
66
71
timeq= mod (timeq,timeq_max)
@@ -70,34 +75,31 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
70
75
71
76
do j= jts,jte
72
77
do i= its,ite
73
- fire_age= time_int/ 3600._kind_phys + (fire_end_hr(i,j)- 1._kind_phys ) ! One hour delay is due to the latency of the RAVE files
74
- fire_age= MAX (0.1_kind_phys ,fire_age) ! in hours
78
+ fire_age= MAX (0.01_kind_phys ,time_int/ 3600 . + (fire_end_hr(i,j)- 2.0 )) ! One hour delay is due to the latency of the RAVE files, hours; one more hour subtracted to have fire_end_hr in the range of 0-24 instead of 0-25
75
79
76
80
SELECT CASE ( fire_type(i,j) ) ! Ag, urban fires, bare land etc.
77
81
CASE (1 )
78
82
! these fires will have exponentially decreasing diurnal cycle,
79
- coef_bb_dc(i,j) = coef_con* 1._kind_phys / (sigma_fire_dur(5 ) * fire_age) * &
80
- exp (- ( log (fire_age) - avg_fire_dur(5 ))** 2 / (2._kind_phys * sigma_fire_dur(5 )** 2 ))
83
+ ! coef_bb_dc(i,j) = coef_con*1._kind_phys/(sigma_fire_dur(5) *fire_age) * &
84
+ ! exp(- ( log(fire_age) - avg_fire_dur(5))**2 /(2._kind_phys*sigma_fire_dur(5)**2 ))
85
+ coef_bb_dc(i,j)= C1/ (sigmx1* fire_age)* exp (- (log (fire_age) - avgx1)** 2 / (2 .* sigmx1** 2 ) )
81
86
82
87
IF ( dbg_opt .AND. time_int< 5000 .) then
83
88
WRITE (6 ,* ) ' i,j,peak_hr(i,j) ' ,i,j,peak_hr(i,j)
84
89
WRITE (6 ,* ) ' coef_bb_dc(i,j) ' ,coef_bb_dc(i,j)
85
90
END IF
86
91
87
- CASE (2 ) ! Savanna and grassland fires
88
- coef_bb_dc(i,j) = coef_con* 1._kind_phys / (sigma_fire_dur(4 ) * fire_age) * &
89
- exp (- ( log (fire_age) - avg_fire_dur(4 ))** 2 / (2._kind_phys * sigma_fire_dur(4 )** 2 ))
92
+ CASE (2 ) ! Savanna and grassland fires, or fires in the eastern US
93
+ ! coef_bb_dc(i,j) = coef_con*1._kind_phys/(sigma_fire_dur(4) *fire_age) * &
94
+ ! exp(- ( log(fire_age) - avg_fire_dur(4))**2 /(2._kind_phys*sigma_fire_dur(4)**2 ))
95
+ coef_bb_dc(i,j)= C2/ (sigmx2* fire_age)* exp (- (log (fire_age) - avgx2)** 2 / (2 .* sigmx2** 2 ) )
90
96
91
97
IF ( dbg_opt .AND. time_int< 5000 .) then
92
98
WRITE (6 ,* ) ' i,j,peak_hr(i,j) ' ,i,j,peak_hr(i,j)
93
99
WRITE (6 ,* ) ' coef_bb_dc(i,j) ' ,coef_bb_dc(i,j)
94
100
END IF
95
101
96
-
97
-
98
- CASE (3 )
99
- ! age_hr= fire_age/3600._kind_phys
100
-
102
+ CASE (3 ,4 ) ! wildfires
101
103
IF (swdown(i,j)<.1 .AND. fire_age> 12 . .AND. fire_hist(i,j)>0.75 ) THEN
102
104
fire_hist(i,j)= 0.75_kind_phys
103
105
ENDIF
@@ -113,15 +115,15 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
113
115
dc_hwp= MAX (0._kind_phys ,dc_hwp)
114
116
dc_hwp= MIN (20._kind_phys ,dc_hwp)
115
117
116
- ! RAR: Gaussian profile for wildfires
117
- dt1= abs (timeq - peak_hr(i,j))
118
- dt2= timeq_max - peak_hr(i,j) + timeq ! peak hour is always <86400.
119
- dtm= MIN (dt1,dt2)
120
- dc_gp = rinti* ( ax2 * exp (- dtm** 2 / (2._kind_phys * cx2** 2 ) ) + const2 - coef2* timeq )
121
- dc_gp = MAX (0._kind_phys ,dc_gp)
118
+ ! RAR: Gaussian profile for wildfires, to be used later
119
+ ! dt1= abs(timeq - peak_hr(i,j))
120
+ ! dt2= timeq_max - peak_hr(i,j) + timeq ! peak hour is always <86400.
121
+ ! dtm= MIN(dt1,dt2)
122
+ ! dc_gp = rinti*( ax2 * exp(- dtm**2/(2._kind_phys*cx2**2) ) + const2 - coef2*timeq )
123
+ ! dc_gp = MAX(0._kind_phys,dc_gp)
122
124
123
125
! dc_fn = MIN(dc_hwp/dc_gp,3._kind_phys)
124
- coef_bb_dc(i,j) = fire_hist(i,j)* dc_hwp
126
+ coef_bb_dc(i,j) = sc_factor * fire_hist(i,j)* dc_hwp ! RAR: scaling factor is applied to the forest fires only, except the eastern US
125
127
126
128
IF ( dbg_opt .AND. time_int< 5000 .) then
127
129
WRITE (6 ,* ) ' i,j,fire_hist(i,j),peak_hr(i,j) ' , i,j,fire_hist(i,j),peak_hr(i,j)
@@ -146,7 +148,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
146
148
if (ebb_dcycle== 1 ) then
147
149
conv= dtstep/ (rho_phy(i,k,j)* dz8w(i,k,j))
148
150
elseif (ebb_dcycle== 2 ) then
149
- conv= sc_factor * coef_bb_dc(i,j)* dtstep/ (rho_phy(i,k,j)* dz8w(i,k,j))
151
+ conv= coef_bb_dc(i,j)* dtstep/ (rho_phy(i,k,j)* dz8w(i,k,j))
150
152
endif
151
153
dm_smoke= conv* ebu(i,k,j)
152
154
chem(i,k,j,p_smoke) = chem(i,k,j,p_smoke) + dm_smoke
0 commit comments