-
Notifications
You must be signed in to change notification settings - Fork 312
/
CNFireLi2016Mod.F90
658 lines (610 loc) · 38.5 KB
/
CNFireLi2016Mod.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
module CNFireLi2016Mod
#include "shr_assert.h"
!-----------------------------------------------------------------------
! !DESCRIPTION:
! module for fire dynamics
! created in Nov, 2012 and revised in Apr, 2013 by F. Li and S. Levis
! based on Li et al. (2012a,b; 2013)
! revised in Apr, 2014 according to Li et al.(2014)
! revised in May, 2015, according to Li et al. (2015, in prep.)
! Fire-related parameters were calibrated or tuned in May, 2015 based on the
! 20th Century transient simulations at f19_g16 with a CLM4.5 version
! (clm50fire), CRUNCEPv5, and climatological lightning data.
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_CL
use shr_const_mod , only : SHR_CONST_PI,SHR_CONST_TKFRZ
use shr_infnan_mod , only : shr_infnan_isnan
use clm_varctl , only : iulog
use clm_varpar , only : nlevdecomp, ndecomp_pools, nlevdecomp_full
use clm_varcon , only : dzsoi_decomp
use pftconMod , only : noveg, pftcon
use abortutils , only : endrun
use decompMod , only : bounds_type
use subgridAveMod , only : p2c
use atm2lndType , only : atm2lnd_type
use CNDVType , only : dgvs_type
use CNVegStateType , only : cnveg_state_type
use CNVegCarbonStateType , only : cnveg_carbonstate_type
use CNVegCarbonFluxType , only : cnveg_carbonflux_type
use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type
use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type
use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con
use EnergyFluxType , only : energyflux_type
use SaturatedExcessRunoffMod , only : saturated_excess_runoff_type
use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type
use Wateratm2lndBulkType , only : wateratm2lndbulk_type
use GridcellType , only : grc
use ColumnType , only : col
use PatchType , only : patch
use SoilBiogeochemStateType , only : get_spinup_latitude_term
use CNFireMethodMod , only : cnfire_method_type
use CNFireBaseMod , only : cnfire_base_type, cnfire_const, cnfire_params
!
implicit none
private
!
! !PUBLIC TYPES:
public :: cnfire_li2016_type
!
type, extends(cnfire_base_type) :: cnfire_li2016_type
private
contains
!
! !PUBLIC MEMBER FUNCTIONS:
procedure, public :: need_lightning_and_popdens
procedure, public :: CNFireArea ! Calculate fire area
end type cnfire_li2016_type
!
! !PRIVATE MEMBER DATA:
!-----------------------------------------------------------------------
character(len=*), parameter, private :: sourcefile = &
__FILE__
contains
!-----------------------------------------------------------------------
function need_lightning_and_popdens(this)
! !ARGUMENTS:
class(cnfire_li2016_type), intent(in) :: this
logical :: need_lightning_and_popdens ! function result
!
! !LOCAL VARIABLES:
character(len=*), parameter :: subname = 'need_lightning_and_popdens'
!-----------------------------------------------------------------------
need_lightning_and_popdens = .true.
end function need_lightning_and_popdens
!-----------------------------------------------------------------------
subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &
atm2lnd_inst, energyflux_inst, saturated_excess_runoff_inst, waterdiagnosticbulk_inst, &
wateratm2lndbulk_inst, cnveg_state_inst, cnveg_carbonstate_inst, totlitc_col, decomp_cpools_vr_col, t_soi17cm_col)
!
! !DESCRIPTION:
! Computes column-level burned area
!
! !USES:
use clm_time_manager , only: get_step_size_real, get_days_per_year, get_curr_date, get_nstep
use clm_varpar , only: max_patch_per_col
use clm_varcon , only: secspday, secsphr
use clm_varctl , only: spinup_state
use pftconMod , only: nc4_grass, nc3crop, ndllf_evr_tmp_tree
use pftconMod , only: nbrdlf_evr_trp_tree, nbrdlf_dcd_trp_tree, nbrdlf_evr_shrub
use dynSubgridControlMod , only : run_has_transient_landcover
!
! !ARGUMENTS:
class(cnfire_li2016_type) :: this
type(bounds_type) , intent(in) :: bounds
integer , intent(in) :: num_soilc ! number of soil columns in filter
integer , intent(in) :: filter_soilc(:) ! filter for soil columns
integer , intent(in) :: num_soilp ! number of soil patches in filter
integer , intent(in) :: filter_soilp(:) ! filter for soil patches
type(atm2lnd_type) , intent(in) :: atm2lnd_inst
type(energyflux_type) , intent(in) :: energyflux_inst
type(saturated_excess_runoff_type) , intent(in) :: saturated_excess_runoff_inst
type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst
type(wateratm2lndbulk_type) , intent(in) :: wateratm2lndbulk_inst
type(cnveg_state_type) , intent(inout) :: cnveg_state_inst
type(cnveg_carbonstate_type) , intent(inout) :: cnveg_carbonstate_inst
real(r8) , intent(in) :: totlitc_col(bounds%begc:)
real(r8) , intent(in) :: decomp_cpools_vr_col(bounds%begc:,1:,1:)
real(r8) , intent(in) :: t_soi17cm_col(bounds%begc:)
!
! !LOCAL VARIABLES:
!
integer :: g,l,c,p,pi,j,fc,fp,kyr, kmo, kda, mcsec ! index variables
real(r8) :: dt ! time step variable (s)
real(r8) :: dayspyr ! days per year
real(r8) :: cli ! effect of climate on deforestation fires (0-1)
real(r8) :: cri ! thresholds used for cli, (mm/d), see Eq.(7) in Li et al.(2013)
real(r8) :: fb ! availability of fuel for regs A and C
real(r8) :: fhd ! impact of hd on agricultural fire
real(r8) :: fgdp ! impact of gdp on agricultural fire
real(r8) :: fire_m ! combustability of fuel for fire occurrence
real(r8) :: spread_m ! combustability of fuel for fire spread
real(r8) :: Lb_lf ! length-to-breadth ratio added by Lifang
integer :: i_cwd ! cwd pool
real(r8) :: lh ! anthro. ignitions (count/km2/hr)
real(r8) :: fs ! hd-dependent fires suppression (0-1)
real(r8) :: ig ! total ignitions (count/km2/hr)
real(r8) :: hdmlf ! human density
real(r8) :: arh, arh30 !combustability of fuel related to RH and RH30
real(r8) :: afuel !weight for arh and arh30
real(r8) :: btran_col(bounds%begc:bounds%endc)
logical :: transient_landcover ! whether this run has any prescribed transient landcover
real(r8), target :: prec60_col_target(bounds%begc:bounds%endc)
real(r8), target :: prec10_col_target(bounds%begc:bounds%endc)
real(r8), target :: rh30_col_target(bounds%begc:bounds%endc)
real(r8), pointer :: prec60_col(:)
real(r8), pointer :: prec10_col(:)
real(r8), pointer :: rh30_col(:)
!-----------------------------------------------------------------------
SHR_ASSERT_ALL_FL((ubound(totlitc_col) == (/bounds%endc/)) , sourcefile, __LINE__)
SHR_ASSERT_ALL_FL((ubound(decomp_cpools_vr_col) == (/bounds%endc,nlevdecomp_full,ndecomp_pools/)), sourcefile, __LINE__)
SHR_ASSERT_ALL_FL((ubound(t_soi17cm_col) == (/bounds%endc/)) , sourcefile, __LINE__)
associate( &
totlitc => totlitc_col , & ! Input: [real(r8) (:) ] (gC/m2) total lit C (column-level mean)
decomp_cpools_vr => decomp_cpools_vr_col , & ! Input: [real(r8) (:,:,:) ] (gC/m3) VR decomp. (litter, cwd, soil)
tsoi17 => t_soi17cm_col , & ! Input: [real(r8) (:) ] (K) soil T for top 0.17 m
lfuel => cnfire_const%lfuel , & ! Input: [real(r8) ] (gC/m2) Lower threshold of fuel mass
ufuel => cnfire_const%ufuel , & ! Input: [real(r8) ] (gC/m2) Upper threshold of fuel mass
rh_hgh => cnfire_const%rh_hgh , & ! Input: [real(r8) ] (%) High relative humidity
rh_low => cnfire_const%rh_low , & ! Input: [real(r8) ] (%) Low relative humidity
bt_min => cnfire_const%bt_min , & ! Input: [real(r8) ] (0-1) Minimum btran
bt_max => cnfire_const%bt_max , & ! Input: [real(r8) ] (0-1) Maximum btran
cli_scale => cnfire_const%cli_scale , & ! Input: [real(r8) ] (/d) global constant for deforestation fires
cropfire_a1 => cnfire_const%cropfire_a1 , & ! Input: [real(r8) ] (/hr) a1 parameter for cropland fire
non_boreal_peatfire_c => cnfire_const%non_boreal_peatfire_c , & ! Input: [real(r8) ] (/hr) c parameter for non-boreal peatland fire
pot_hmn_ign_counts_alpha => cnfire_const%pot_hmn_ign_counts_alpha , & ! Input: [real(r8) ] (/person/month) Potential human ignition counts
boreal_peatfire_c => cnfire_const%boreal_peatfire_c , & ! Input: [real(r8) ] (/hr) c parameter for boreal peatland fire
fsr_pft => pftcon%fsr_pft , & ! Input:
fd_pft => pftcon%fd_pft , & ! Input:
btran2 => energyflux_inst%btran2_patch , & ! Input: [real(r8) (:) ] root zone soil wetness
fsat => saturated_excess_runoff_inst%fsat_col , & ! Input: [real(r8) (:) ] fractional area with water table at surface
wf2 => waterdiagnosticbulk_inst%wf2_col , & ! Input: [real(r8) (:) ] soil water as frac. of whc for top 0.17 m
is_cwd => decomp_cascade_con%is_cwd , & ! Input: [logical (:) ] TRUE => pool is a cwd pool
spinup_factor => decomp_cascade_con%spinup_factor , & ! Input: [real(r8) (:) ] factor for AD spinup associated with each pool
forc_rh => wateratm2lndbulk_inst%forc_rh_grc , & ! Input: [real(r8) (:) ] relative humidity
forc_wind => atm2lnd_inst%forc_wind_grc , & ! Input: [real(r8) (:) ] atmospheric wind speed (m/s)
forc_t => atm2lnd_inst%forc_t_downscaled_col , & ! Input: [real(r8) (:) ] downscaled atmospheric temperature (Kelvin)
forc_rain => wateratm2lndbulk_inst%forc_rain_downscaled_col , & ! Input: [real(r8) (:) ] downscaled rain
forc_snow => wateratm2lndbulk_inst%forc_snow_downscaled_col , & ! Input: [real(r8) (:) ] downscaled snow
prec60 => wateratm2lndbulk_inst%prec60_patch , & ! Input: [real(r8) (:) ] 60-day running mean of tot. precipitation
prec10 => wateratm2lndbulk_inst%prec10_patch , & ! Input: [real(r8) (:) ] 10-day running mean of tot. precipitation
rh30 => wateratm2lndbulk_inst%rh30_patch , & ! Input: [real(r8) (:) ] 10-day running mean of tot. precipitation
dwt_smoothed => cnveg_state_inst%dwt_smoothed_patch , & ! Input: [real(r8) (:) ] change in patch weight (-1 to 1) on the gridcell, smoothed over the year
cropf_col => cnveg_state_inst%cropf_col , & ! Input: [real(r8) (:) ] cropland fraction in veg column
gdp_lf => cnveg_state_inst%gdp_lf_col , & ! Input: [real(r8) (:) ] gdp data
peatf_lf => cnveg_state_inst%peatf_lf_col , & ! Input: [real(r8) (:) ] peatland fraction data
abm_lf => cnveg_state_inst%abm_lf_col , & ! Input: [integer (:) ] prescribed crop fire time
baf_crop => cnveg_state_inst%baf_crop_col , & ! Output: [real(r8) (:) ] burned area fraction for cropland (/sec)
baf_peatf => cnveg_state_inst%baf_peatf_col , & ! Output: [real(r8) (:) ] burned area fraction for peatland (/sec)
burndate => cnveg_state_inst%burndate_patch , & ! Output: [integer (:) ] burn date for crop
fbac => cnveg_state_inst%fbac_col , & ! Output: [real(r8) (:) ] total burned area out of conversion (/sec)
fbac1 => cnveg_state_inst%fbac1_col , & ! Output: [real(r8) (:) ] burned area out of conversion region due to land use fire
farea_burned => cnveg_state_inst%farea_burned_col , & ! Output: [real(r8) (:) ] total fractional area burned (/sec)
nfire => cnveg_state_inst%nfire_col , & ! Output: [real(r8) (:) ] fire counts (count/km2/sec), valid only in Reg. C
fsr_col => cnveg_state_inst%fsr_col , & ! Output: [real(r8) (:) ] fire spread rate at column level
fd_col => cnveg_state_inst%fd_col , & ! Output: [real(r8) (:) ] fire duration rate at column level
lgdp_col => cnveg_state_inst%lgdp_col , & ! Output: [real(r8) (:) ] gdp limitation factor for nfire
lgdp1_col => cnveg_state_inst%lgdp1_col , & ! Output: [real(r8) (:) ] gdp limitation factor for baf per fire
lpop_col => cnveg_state_inst%lpop_col , & ! Output: [real(r8) (:) ] pop limitation factor for baf per fire
lfwt => cnveg_state_inst%lfwt_col , & ! Output: [real(r8) (:) ] fractional coverage of non-crop and non-bare-soil Patches
trotr1_col => cnveg_state_inst%trotr1_col , & ! Output: [real(r8) (:) ] patch weight of BET on the column (0-1)
trotr2_col => cnveg_state_inst%trotr2_col , & ! Output: [real(r8) (:) ] patch weight of BDT on the column (0-1)
dtrotr_col => cnveg_state_inst%dtrotr_col , & ! Output: [real(r8) (:) ] decreased frac. coverage of BET+BDT on grid for dt
lfc => cnveg_state_inst%lfc_col , & ! Output: [real(r8) (:) ] conversion area frac. of BET+BDT that haven't burned before
wtlf => cnveg_state_inst%wtlf_col , & ! Output: [real(r8) (:) ] fractional coverage of non-crop Patches
totvegc => cnveg_carbonstate_inst%totvegc_col , & ! Input: [real(r8) (:) ] totvegc at column level
deadcrootc => cnveg_carbonstate_inst%deadcrootc_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead coarse root C
deadcrootc_storage => cnveg_carbonstate_inst%deadcrootc_storage_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead coarse root C storage
deadcrootc_xfer => cnveg_carbonstate_inst%deadcrootc_xfer_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead coarse root C transfer
deadstemc => cnveg_carbonstate_inst%deadstemc_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead stem root C
frootc => cnveg_carbonstate_inst%frootc_patch , & ! Input: [real(r8) (:) ] (gC/m2) fine root C
frootc_storage => cnveg_carbonstate_inst%frootc_storage_patch , & ! Input: [real(r8) (:) ] (gC/m2) fine root C storage
frootc_xfer => cnveg_carbonstate_inst%frootc_xfer_patch , & ! Input: [real(r8) (:) ] (gC/m2) fine root C transfer
livecrootc => cnveg_carbonstate_inst%livecrootc_patch , & ! Input: [real(r8) (:) ] (gC/m2) live coarse root C
livecrootc_storage => cnveg_carbonstate_inst%livecrootc_storage_patch , & ! Input: [real(r8) (:) ] (gC/m2) live coarse root C storage
livecrootc_xfer => cnveg_carbonstate_inst%livecrootc_xfer_patch , & ! Input: [real(r8) (:) ] (gC/m2) live coarse root C transfer
leafc => cnveg_carbonstate_inst%leafc_patch , & ! Input: [real(r8) (:) ] (gC/m2) leaf C
leafc_storage => cnveg_carbonstate_inst%leafc_storage_patch , & ! Input: [real(r8) (:) ] (gC/m2) leaf C storage
leafc_xfer => cnveg_carbonstate_inst%leafc_xfer_patch , & ! Input: [real(r8) (:) ] (gC/m2) leaf C transfer
rootc_col => cnveg_carbonstate_inst%rootc_col , & ! Output: [real(r8) (:) ] root carbon
leafc_col => cnveg_carbonstate_inst%leafc_col , & ! Output: [real(r8) (:) ] leaf carbon at column level
deadstemc_col => cnveg_carbonstate_inst%deadstemc_col , & ! Output: [real(r8) (:) ] deadstem carbon at column level
fuelc => cnveg_carbonstate_inst%fuelc_col , & ! Output: [real(r8) (:) ] fuel load coutside cropland
fuelc_crop => cnveg_carbonstate_inst%fuelc_crop_col & ! Output: [real(r8) (:) ] fuel load for cropland
)
transient_landcover = run_has_transient_landcover()
!pft to column average
prec10_col =>prec10_col_target
call p2c(bounds, num_soilc, filter_soilc, &
prec10(bounds%begp:bounds%endp), &
prec10_col(bounds%begc:bounds%endc))
prec60_col =>prec60_col_target
call p2c(bounds, num_soilc, filter_soilc, &
prec60(bounds%begp:bounds%endp), &
prec60_col(bounds%begc:bounds%endc))
rh30_col =>rh30_col_target
call p2c(bounds, num_soilc, filter_soilc, &
rh30(bounds%begp:bounds%endp), &
rh30_col(bounds%begc:bounds%endc))
call p2c(bounds, num_soilc, filter_soilc, &
leafc(bounds%begp:bounds%endp), &
leafc_col(bounds%begc:bounds%endc))
call p2c(bounds, num_soilc, filter_soilc, &
deadstemc(bounds%begp:bounds%endp), &
deadstemc_col(bounds%begc:bounds%endc))
call get_curr_date (kyr, kmo, kda, mcsec)
dayspyr = get_days_per_year()
! Get model step size
dt = get_step_size_real()
!
! On first time-step, just set area burned to zero and exit
!
if ( get_nstep() == 0 )then
do fc = 1,num_soilc
c = filter_soilc(fc)
farea_burned(c) = 0._r8
baf_crop(c) = 0._r8
baf_peatf(c) = 0._r8
fbac(c) = 0._r8
fbac1(c) = 0._r8
cropf_col(c) = 0._r8
end do
return
end if
!
! Calculate fraction of crop (cropf_col) and non-crop and non-bare-soil
! vegetation (lfwt) in vegetated column
!
do fc = 1,num_soilc
c = filter_soilc(fc)
cropf_col(c) = 0._r8
lfwt(c) = 0._r8
end do
do pi = 1,max_patch_per_col
do fc = 1,num_soilc
c = filter_soilc(fc)
if (pi <= col%npatches(c)) then
p = col%patchi(c) + pi - 1
! For crop veg types
if( patch%itype(p) > nc4_grass )then
cropf_col(c) = cropf_col(c) + patch%wtcol(p)
end if
! For natural vegetation (non-crop and non-bare-soil)
if( patch%itype(p) >= ndllf_evr_tmp_tree .and. patch%itype(p) <= nc4_grass )then
lfwt(c) = lfwt(c) + patch%wtcol(p)
end if
end if
end do
end do
!
! Calculate crop fuel
!
do fc = 1,num_soilc
c = filter_soilc(fc)
fuelc_crop(c)=0._r8
end do
do pi = 1,max_patch_per_col
do fc = 1,num_soilc
c = filter_soilc(fc)
if (pi <= col%npatches(c)) then
p = col%patchi(c) + pi - 1
! For crop PFTs, fuel load includes leaf and litter; only
! column-level litter carbon
! is available, so we use leaf carbon to estimate the
! litter carbon for crop PFTs
if( patch%itype(p) > nc4_grass .and. patch%wtcol(p) > 0._r8 .and. leafc_col(c) > 0._r8 )then
fuelc_crop(c)=fuelc_crop(c) + (leafc(p) + leafc_storage(p) + &
leafc_xfer(p))*patch%wtcol(p)/cropf_col(c) + &
totlitc(c)*leafc(p)/leafc_col(c)*patch%wtcol(p)/cropf_col(c)
end if
end if
end do
end do
!
! Calculate noncrop column variables
!
do fc = 1,num_soilc
c = filter_soilc(fc)
fsr_col(c) = 0._r8
fd_col(c) = 0._r8
rootc_col(c) = 0._r8
lgdp_col(c) = 0._r8
lgdp1_col(c) = 0._r8
lpop_col(c) = 0._r8
btran_col(c) = 0._r8
wtlf(c) = 0._r8
trotr1_col(c)= 0._r8
trotr2_col(c)= 0._r8
if (transient_landcover) then
dtrotr_col(c)=0._r8
end if
end do
do pi = 1,max_patch_per_col
do fc = 1,num_soilc
c = filter_soilc(fc)
g = col%gridcell(c)
if (pi <= col%npatches(c)) then
p = col%patchi(c) + pi - 1
! For non-crop -- natural vegetation and bare-soil
if( patch%itype(p) < nc3crop .and. cropf_col(c) < 1.0_r8 )then
if( .not. shr_infnan_isnan(btran2(p))) then
if (btran2(p) <= 1._r8 ) then
btran_col(c) = btran_col(c)+btran2(p)*patch%wtcol(p)
wtlf(c) = wtlf(c)+patch%wtcol(p)
end if
end if
! NOTE(wjs, 2016-12-15) These calculations of the fraction of evergreen
! and deciduous tropical trees (used to determine if a column is
! tropical closed forest) use the current fractions. However, I think
! they are used in code that applies to land cover change. Note that
! land cover change is currently generated on the first time step of the
! year (even though the fire code sees the annually-smoothed dwt). Thus,
! I think that, for this to be totally consistent, this code should
! consider the fractional coverage of each PFT prior to the relevant
! land cover change event. (These fractions could be computed in the
! code that handles land cover change, so that the fire code remains
! agnostic to exactly how and when land cover change happens.)
!
! For example, if a year started with fractional coverages of
! nbrdlf_evr_trp_tree = 0.35 and nbrdlf_dcd_trp_tree = 0.35, but then
! the start-of-year land cover change reduced both of these to 0.2: The
! current code would consider the column to NOT be tropical closed
! forest (because nbrdlf_evr_trp_tree+nbrdlf_dcd_trp_tree < 0.6),
! whereas in fact the land cover change occurred when the column *was*
! tropical closed forest.
if( patch%itype(p) == nbrdlf_evr_trp_tree .and. patch%wtcol(p) > 0._r8 )then
trotr1_col(c)=trotr1_col(c)+patch%wtcol(p)
end if
if( patch%itype(p) == nbrdlf_dcd_trp_tree .and. patch%wtcol(p) > 0._r8 )then
trotr2_col(c)=trotr2_col(c)+patch%wtcol(p)
end if
if (transient_landcover) then
if( patch%itype(p) == nbrdlf_evr_trp_tree .or. patch%itype(p) == nbrdlf_dcd_trp_tree )then
if(dwt_smoothed(p) < 0._r8)then
! Land cover change in CLM happens all at once on the first time
! step of the year. However, the fire code needs deforestation
! rates throughout the year, in order to combine these
! deforestation rates with the current season's climate. So we
! use a smoothed version of dwt.
!
! This isn't ideal, because the carbon stocks that the fire code
! is operating on will have decreased by the full annual amount
! before the fire code does anything. But the biggest effect of
! these deforestation fires is as a trigger for other fires, and
! the C fluxes are merely diagnostic so don't need to be
! conservative, so this isn't a big issue.
!
! (Actually, it would be even better if the fire code had a
! realistic breakdown of annual deforestation into the
! different seasons. But having deforestation spread evenly
! throughout the year is much better than having it all
! concentrated on January 1.)
dtrotr_col(c)=dtrotr_col(c)-dwt_smoothed(p)
end if
end if
end if
if (spinup_state == 2) then
rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + &
frootc_xfer(p) + deadcrootc(p) * 10._r8 + &
deadcrootc_storage(p) + deadcrootc_xfer(p) + &
livecrootc(p)+livecrootc_storage(p) + &
livecrootc_xfer(p))*patch%wtcol(p)
else
rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + &
frootc_xfer(p) + deadcrootc(p) + &
deadcrootc_storage(p) + deadcrootc_xfer(p) + &
livecrootc(p)+livecrootc_storage(p) + &
livecrootc_xfer(p))*patch%wtcol(p)
endif
fsr_col(c) = fsr_col(c) + fsr_pft(patch%itype(p))*patch%wtcol(p)/(1.0_r8-cropf_col(c))
hdmlf=this%forc_hdm(g)
! all these constants are in Li et al. BG (2012a,b;2013)
if( hdmlf > 0.1_r8 )then
! For NOT bare-soil
if( patch%itype(p) /= noveg )then
! For shrub and grass (crop already excluded above)
if( patch%itype(p) >= nbrdlf_evr_shrub )then !for shurb and grass
lgdp_col(c) = lgdp_col(c) + (0.1_r8 + 0.9_r8* &
exp(-1._r8*SHR_CONST_PI* &
(gdp_lf(c)/8._r8)**0.5_r8))*patch%wtcol(p) &
/(1.0_r8 - cropf_col(c))
lgdp1_col(c) = lgdp1_col(c) + (0.2_r8 + 0.8_r8* &
exp(-1._r8*SHR_CONST_PI* &
(gdp_lf(c)/7._r8)))*patch%wtcol(p)/(1._r8 - cropf_col(c))
lpop_col(c) = lpop_col(c) + (0.2_r8 + 0.8_r8* &
exp(-1._r8*SHR_CONST_PI* &
(hdmlf/450._r8)**0.5_r8))*patch%wtcol(p)/(1._r8 - cropf_col(c))
else ! for trees
if( gdp_lf(c) > 20._r8 )then
lgdp_col(c) =lgdp_col(c)+cnfire_const%occur_hi_gdp_tree*patch%wtcol(p)/(1._r8 - cropf_col(c))
lgdp1_col(c) =lgdp1_col(c)+0.62_r8*patch%wtcol(p)/(1._r8 - cropf_col(c))
else
if( gdp_lf(c) > 8._r8 )then
lgdp_col(c)=lgdp_col(c)+0.79_r8*patch%wtcol(p)/(1._r8 - cropf_col(c))
lgdp1_col(c)=lgdp1_col(c)+0.83_r8*patch%wtcol(p)/(1._r8 - cropf_col(c))
else
lgdp_col(c) = lgdp_col(c)+patch%wtcol(p)/(1._r8 - cropf_col(c))
lgdp1_col(c)=lgdp1_col(c)+patch%wtcol(p)/(1._r8 - cropf_col(c))
end if
end if
lpop_col(c) = lpop_col(c) + (0.4_r8 + 0.6_r8* &
exp(-1._r8*SHR_CONST_PI* &
(hdmlf/125._r8)))*patch%wtcol(p)/(1._r8 -cropf_col(c))
end if
end if
else
lgdp_col(c) = lgdp_col(c)+patch%wtcol(p)/(1.0_r8 - cropf_col(c))
lgdp1_col(c) = lgdp1_col(c)+patch%wtcol(p)/(1.0_r8 -cropf_col(c))
lpop_col(c) = lpop_col(c)+patch%wtcol(p)/(1.0_r8 -cropf_col(c))
end if
fd_col(c) = fd_col(c) + fd_pft(patch%itype(p)) * patch%wtcol(p) * secsphr / (1.0_r8-cropf_col(c))
end if
end if
end do
end do
! estimate annual decreased fractional coverage of BET+BDT
! land cover conversion in CLM4.5 is the same for each timestep except for the beginning
if (transient_landcover) then
do fc = 1,num_soilc
c = filter_soilc(fc)
if( dtrotr_col(c) > 0._r8 )then
if( kmo == 1 .and. kda == 1 .and. mcsec == 0)then
lfc(c) = 0._r8
end if
if( kmo == 1 .and. kda == 1 .and. mcsec == dt)then
lfc(c) = dtrotr_col(c)*dayspyr*secspday/dt
end if
else
lfc(c)=0._r8
end if
end do
end if
!
! calculate burned area fraction in cropland
!
do fc = 1,num_soilc
c = filter_soilc(fc)
baf_crop(c)=0._r8
end do
do fp = 1,num_soilp
p = filter_soilp(fp)
if( kmo == 1 .and. kda == 1 .and. mcsec == 0 )then
burndate(p) = 10000 ! init. value; actual range [0 365]
end if
end do
do pi = 1,max_patch_per_col
do fc = 1,num_soilc
c = filter_soilc(fc)
g= col%gridcell(c)
hdmlf=this%forc_hdm(g)
if (pi <= col%npatches(c)) then
p = col%patchi(c) + pi - 1
! For crop
if( forc_t(c) >= SHR_CONST_TKFRZ .and. patch%itype(p) > nc4_grass .and. &
kmo == abm_lf(c) .and. &
burndate(p) >= 999 .and. patch%wtcol(p) > 0._r8 )then ! catch crop burn time
! calculate human density impact on ag. fire
fhd = 0.04_r8+0.96_r8*exp(-1._r8*SHR_CONST_PI*(hdmlf/350._r8)**0.5_r8)
! calculate impact of GDP on ag. fire
fgdp = 0.01_r8+0.99_r8*exp(-1._r8*SHR_CONST_PI*(gdp_lf(c)/10._r8))
! calculate burned area
fb = max(0.0_r8,min(1.0_r8,(fuelc_crop(c)-lfuel)/(ufuel-lfuel)))
! crop fire only for generic crop types at this time
! managed crops are treated as grasses if crop model is turned on
baf_crop(c) = baf_crop(c) + cropfire_a1/secsphr*fhd*fgdp*patch%wtcol(p)
if( fb*fhd*fgdp*patch%wtcol(p) > 0._r8)then
burndate(p)=kda
end if
end if
end if
end do
end do
!
! calculate peatland fire
!
do fc = 1, num_soilc
c = filter_soilc(fc)
g= col%gridcell(c)
if(grc%latdeg(g) < cnfire_const%borealat )then
baf_peatf(c) = non_boreal_peatfire_c/secsphr*max(0._r8, &
min(1._r8,(4.0_r8-prec60_col(c)*secspday)/ &
4.0_r8))**2*peatf_lf(c)*(1._r8-fsat(c))
else
baf_peatf(c) = boreal_peatfire_c/secsphr*exp(-SHR_CONST_PI*(max(wf2(c),0._r8)/0.3_r8))* &
max(0._r8,min(1._r8,(tsoi17(c)-SHR_CONST_TKFRZ)/10._r8))*peatf_lf(c)* &
(1._r8-fsat(c))
end if
end do
!
! calculate other fires
!
! Set the number of timesteps for e-folding.
! When the simulation has run fewer than this number of steps,
! re-scale the e-folding time to get a stable early estimate.
! find which pool is the cwd pool
i_cwd = 0
do l = 1, ndecomp_pools
if ( is_cwd(l) ) then
i_cwd = l
endif
end do
!
! begin column loop to calculate fractional area affected by fire
!
do fc = 1, num_soilc
c = filter_soilc(fc)
g = col%gridcell(c)
hdmlf=this%forc_hdm(g)
nfire(c) = 0._r8
if( cropf_col(c) < 1._r8 )then
fuelc(c) = totlitc(c)+totvegc(c)-rootc_col(c)-fuelc_crop(c)*cropf_col(c)
if (spinup_state == 2) then
fuelc(c) = fuelc(c) + ((10._r8 - 1._r8)*deadstemc_col(c))
do j = 1, nlevdecomp
fuelc(c) = fuelc(c)+decomp_cpools_vr(c,j,i_cwd) * dzsoi_decomp(j) * spinup_factor(i_cwd) &
* get_spinup_latitude_term(grc%latdeg(col%gridcell(c)))
end do
else
do j = 1, nlevdecomp
fuelc(c) = fuelc(c)+decomp_cpools_vr(c,j,i_cwd) * dzsoi_decomp(j)
end do
end if
fuelc(c) = fuelc(c)/(1._r8-cropf_col(c))
fb = max(0.0_r8,min(1.0_r8,(fuelc(c)-lfuel)/(ufuel-lfuel)))
if (trotr1_col(c)+trotr2_col(c)<=0.6_r8) then
afuel =min(1._r8,max(0._r8,(fuelc(c)-2500._r8)/(5000._r8-2500._r8)))
arh=1._r8-max(0._r8, min(1._r8,(forc_rh(g)-rh_low)/(rh_hgh-rh_low)))
arh30=1._r8-max(cnfire_params%prh30, min(1._r8,rh30_col(c)/90._r8))
if (forc_rh(g) < rh_hgh.and. wtlf(c) > 0._r8 .and. tsoi17(c)> SHR_CONST_TKFRZ)then
fire_m = ((afuel*arh30+(1._r8-afuel)*arh)**1.5_r8)*((1._r8 -max(0._r8,&
min(1._r8,(btran_col(c)/wtlf(c)-bt_min)/(bt_max-bt_min))))**0.5_r8)
else
fire_m = 0._r8
end if
lh = pot_hmn_ign_counts_alpha*6.8_r8*hdmlf**(0.43_r8)/30._r8/24._r8
fs = 1._r8-(0.01_r8+0.98_r8*exp(-0.025_r8*hdmlf))
ig = (lh+this%forc_lnfm(g)/(5.16_r8+2.16_r8*cos(SHR_CONST_PI/180._r8*3*min(60._r8,abs(grc%latdeg(g)))))* &
cnfire_params%ignition_efficiency)*(1._r8-fs)*(1._r8-cropf_col(c))
nfire(c) = ig/secsphr*fb*fire_m*lgdp_col(c) !fire counts/km2/sec
Lb_lf = 1._r8+10._r8*(1._r8-EXP(-0.06_r8*forc_wind(g)))
spread_m = fire_m**0.5_r8
farea_burned(c) = min(1._r8,(cnfire_const%g0*spread_m*fsr_col(c)* &
fd_col(c)/1000._r8)**2*lgdp1_col(c)* &
lpop_col(c)*nfire(c)*SHR_CONST_PI*Lb_lf+ &
baf_crop(c)+baf_peatf(c)) ! fraction (0-1) per sec
else
farea_burned(c)=min(1._r8,baf_crop(c)+baf_peatf(c))
end if
!
! if landuse change data is used, calculate deforestation fires and
! add it in the total of burned area fraction
!
if (transient_landcover) then
if( trotr1_col(c)+trotr2_col(c) > 0.6_r8 )then
if(( kmo == 1 .and. kda == 1 .and. mcsec == 0) .or. &
dtrotr_col(c) <=0._r8 )then
fbac1(c) = 0._r8
farea_burned(c) = baf_crop(c)+baf_peatf(c)
else
cri = (4.0_r8*trotr1_col(c)+1.8_r8*trotr2_col(c))/(trotr1_col(c)+trotr2_col(c))
cli = (max(0._r8,min(1._r8,(cri-prec60_col(c)*secspday)/cri))**0.5)* &
(max(0._r8,min(1._r8,(cri-prec10_col(c)*secspday)/cri))**0.5)* &
max(0.0005_r8,min(1._r8,19._r8*dtrotr_col(c)*dayspyr*secspday/dt-0.001_r8))* &
max(0._r8,min(1._r8,(0.25_r8-(forc_rain(c)+forc_snow(c))*secsphr)/0.25_r8))
farea_burned(c) = cli*(cli_scale/secspday)+baf_crop(c)+baf_peatf(c)
! burned area out of conversion region due to land use fire
fbac1(c) = max(0._r8,fb*cli*(cli_scale/secspday) - 2.0_r8*lfc(c)/dt)
end if
! total burned area out of conversion
fbac(c) = fbac1(c)+baf_crop(c)+baf_peatf(c)
else
fbac(c) = farea_burned(c)
end if
end if
else
farea_burned(c) = min(1._r8,baf_crop(c)+baf_peatf(c))
end if
end do ! end of column loop
end associate
end subroutine CNFireArea
end module CNFireLi2016Mod