PageRenderTime 56ms CodeModel.GetById 22ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_fr_sfire_phys.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1755 lines | 965 code | 237 blank | 553 comment | 75 complexity | bc832edc5d8b31d4f09e343acb31d622 MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

    

  1. !
    2. 

  2. 

  3. !*** Jan Mandel August-October 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com
    4. 

  4. !
    5. 

  5. ! This file contains parts copied and/or adapted from earlier codes by 
    6. 

  6. ! Terry Clark, Janice Coen, Don Latham, and Net Patton.
    7. 

  7. 

  8. module module_fr_sfire_phys
    9. 

  9. 

  10. use module_model_constants, only: cp,xlv
    11. 

  11. use module_fr_sfire_util
    12. 

  12. 

  13. implicit none
    14. 

  14. PRIVATE
    15. 

  15. 

  16. 

  17. type fire_params
    18. 

  18. real,pointer,dimension(:,:):: vx,vy                ! wind velocity (m/s)
    19. 

  19. real,pointer,dimension(:,:):: zsf                  ! terrain height (m)
    20. 

  20. real,pointer,dimension(:,:):: dzdxf,dzdyf          ! terrain grad (1)
    21. 

  21. real,pointer,dimension(:,:):: bbb,phisc,phiwc,r_0 ! spread formula coefficients
    22. 

  22. real,pointer,dimension(:,:):: fgip                 ! init mass of surface fuel (kg/m^2)
    23. 

  23. real,pointer,dimension(:,:):: ischap               ! 1 if chapparal
    24. 

  24. real,pointer,dimension(:,:):: fuel_time            ! time to burn to 1/e (s)
    25. 

  25. real,pointer,dimension(:,:):: fmc_g                ! fuel moisture contents, ground (1)
    26. 

  26. real,pointer,dimension(:,:):: nfuel_cat            ! fuel category (integer values)
    27. 

  27. end type fire_params
    28. 

  28. 

  29. ! subroutines and functions 
    30. 

  30. PUBLIC:: init_fuel_cats,fire_ros,heat_fluxes,set_nfuel_cat,set_fire_params, &
    31. 

  31. write_fuels_m,fire_risk,fire_intensity,fuel_moisture,advance_moisture,fuel_name,&
    32. 

  32. fire_rate_of_spread
    33. 

  33. 

  34. ! types
    35. 

  35. public:: fire_params
    36. 

  36. 

  37. ! variables
    38. 

  38. PUBLIC:: fire_wind_height,fcz0,fcwh,have_fuel_cats,nfuelcats,no_fuel_cat,no_fuel_cat2,windrf,moisture_classes
    39. 

  39. PUBLIC:: mfuelcats
    40. 

  40. ! NOTE: fcwh and fcz0 are called fwh and fz0 in read/write statements
    41. 

  41. 

  42. 

  43. ! moisture behavior, see Mandel et al EGU 2012
    44. 

  44. 

  45. !! To add moisture classes:
    46. 

  46. !  1. change parameter max_moisture_classes below
    47. 

  47. !  2. change the default of nfmc to the same value in Registry/registry.fire
    48. 

  48. !  3. add the appropriate lines real::fmc_gw<number>= <default values>
    49. 

  49. !  4. add dfault 
    50. 

  50. 

  51. !*** dimensions
    52. 

  52.    INTEGER, PARAMETER :: mfuelcats = 30     ! max number of fuel categories 
    53. 

  53.   integer, parameter::max_moisture_classes=5
    54. 

  54. !***
    55. 

  55. 

  56.   integer, parameter::zm=max_moisture_classes - 3
    57. 

  57.   integer:: moisture_classes=3
    58. 

  58.   real, dimension(max_moisture_classes):: drying_lag,wetting_lag,saturation_moisture,saturation_rain, &
    59. 

  59.          rain_threshold,rec_drying_lag_sec,rec_wetting_lag_sec
    60. 

  60.   integer, dimension(max_moisture_classes):: drying_model,wetting_model,fmc_gc_initialization
    61. 

  61.    ! relative weights of moisture class for each fuel category
    62. 

  62.    integer::itmp
    63. 

  63.    CHARACTER (len=80), DIMENSION(max_moisture_classes), save :: moisture_class_name
    64. 

  64.    real, dimension(mfuelcats):: &  ! should sum up to one
    65. 

  65.      fmc_gw01=(/ (1.0, itmp=1,mfuelcats) /), & 
    66. 

  66.      fmc_gw02=(/ (0.0, itmp=1,mfuelcats) /), & 
    67. 

  67.      fmc_gw03=(/ (0.0, itmp=1,mfuelcats) /), & 
    68. 

  68.      fmc_gw04=(/ (0.0, itmp=1,mfuelcats) /), & 
    69. 

  69.      fmc_gw05=(/ (0.0, itmp=1,mfuelcats) /) 
    70. 

  70.     
    71. 

  71.   data moisture_class_name /'1-hour fuel','10-hour fuel','100-hour fuel',zm*'NOT USED'/
    72. 

  72.   data drying_lag                 /1., 10., 100. , zm*0./  ! time lag (h) approaching equilibrium moisture 
    73. 

  73.   data wetting_lag                   /14, 140., 1400., zm*0./  ! time lag (h) for approaching saturation in rain
    74. 

  74.   data saturation_moisture        /2.5, 2.5, 2.5 , zm*0./  ! saturation moisture contents (1) in rain
    75. 

  75.   data saturation_rain            /8.0, 8.0, 8.0 , zm*0./  ! stronger rain matters only in duration (mm/h) 
    76. 

  76.   data rain_threshold             /0.05, 0.05, 0.05, zm*0 /! rain intensity this small is same as nothing
    77. 

  77.   data drying_model               /1,   1,   1,    zm*1 /
    78. 

  78.   data wetting_model              /1,   1,   1,    zm*1 /
    79. 

  79.   data fmc_gc_initialization      /2,   2,   2,    zm*2 /  ! initialization 0=input, 1=from fuelmc_g, 2=from equilibrium
    80. 

  80.   real, dimension(7)::eq_p
    81. 

  81.   data eq_p/  1.035e-09, &   !(3.893e-10, 1.681e-09) ! coefficients of the equilibrium fuel moisture polynomial
    82. 

  82.               -2.62e-07, &   !(-4.593e-07, -6.473e-08) ! fitted from the graph in Schroeder and Buck (1970)
    83. 

  83.               2.507e-05, &   !(2.194e-06, 4.795e-05)
    84. 

  84.               -0.001107, &   !(-0.002353, 0.000139)
    85. 

  85.                 0.02245, &   !(-0.009188, 0.05409)
    86. 

  86.                -0.05901, &   !(-0.3721, 0.254)
    87. 

  87.                   3.043/   !(2.17, 3.915)
    88. 

  88.    
    89. 

  89.    ! =========================================================================
    90. 

  90. 

  91. ! Following table copied from module_fr_cawfe_fuel by Ned Patton with minor changes. 
    92. 

  92. ! Based on:           Clark, T. L., J. L. Coen and D. Latham: 2004, 
    93. 

  93. !                      "Description of a coupled atmosphere-fire model",
    94. 

  94. !                      International Journal of Wildland Fire, 13, 49-63.
    95. 

  95. !
    96. 

  96. ! edited by Jan Mandel   jmandel@ucar.edu  September 2007
    97. 

  97. !
    98. 

  98. ! - moved all fuel related constants and the initialization subroutine here
    99. 

  99. ! - copied descriptions for fuel categories from fire_sfc.m4 in the original CAWFE code 
    100. 

  100. ! This file had to be copied under a new name because packages in wrf physics
    101. 

  101. ! layer are not allowed to call each other.
    102. 

  102. 

  103. !D in col 2 means quantity derived from the others
    104. 

  104. !
    105. 

  105. !  Scalar constants (data same for all fuel categories):
    106. 

  106. !       HFGL           SURFACE FIRE HEAT FLUX THRESHOLD TO IGNITE CANOPY (W/m^2)
    107. 

  107. !       CMBCNST        JOULES PER KG OF DRY FUEL
    108. 

  108. !       FUELHEAT       FUEL PARTICLE LOW HEAT CONTENT, BTU/LB
    109. 

  109. !       FUELMC_G       FUEL PARTICLE (SURFACE) MOISTURE CONTENT
    110. 

  110. !D      BMST           RATIO OF LATENT TO SENSIBLE HEAT FROM SFC BURN:
    111. 

  111. !                        % of total fuel mass that is water (not quite
    112. 

  112. !                        = % fuel moisture).    BMST= (H20)/(H20+DRY)
    113. 

  113. !                        so BMST = FUELMC_G / (1 + FUELMC_G)  where
    114. 

  114. !                        FUELMC_G = ground fuel moisture
    115. 

  115. !
    116. 

  116. !  Data arrays indexed by fuel category:
    117. 

  117. !       FGI            INITIAL TOTAL MASS OF SURFACE FUEL (KG/M**2)
    118. 

  118. !       FUELDEPTHM     FUEL DEPTH, IN M  (CONVERTED TO FT)              
    119. 

  119. !       SAVR           FUEL PARTICLE SURFACE-AREA-TO-VOLUME RATIO, 1/FT
    120. 

  120. !                         GRASS: 3500., 10 hr fuel: 109.,  100 hr fuel: 30.
    121. 

  121. !       FUELMCE        MOISTURE CONTENT OF EXTINCTION; 0.30 FOR MANY DEAD FUELS; 0.15 FOR GRASS
    122. 

  122. !       FUELDENS       OVENDRY PARTICLE DENSITY, LB/FT^3
    123. 

  123. !       ST             FUEL PARTICLE TOTAL MINERAL CONTENT
    124. 

  124. !       SE             FUEL PARTICLE EFFECTIVE MINERAL CONTENT
    125. 

  125. !       WEIGHT         WEIGHTING PARAMETER THAT DETERMINES THE SLOPE OF THE MASS LOSS CURVE
    126. 

  126. !                        RANGES FROM ~5 (FAST BURNUP) TO 1000 ( ~40% DECR OVER 10 MIN).
    127. 

  127. !       FCI_D          INITIAL DRY MASS OF CANOPY FUEL
    128. 

  128. !       FCT            BURN OUT TIME FOR CANOPY FUEL, AFTER DRY (S)
    129. 

  129. !       ichap          1 if chaparral, 0 if not
    130. 

  130. !D      FCI            INITIAL TOTAL MASS OF CANOPY FUEL
    131. 

  131. !D      FCBR           FUEL CANOPY BURN RATE (KG/M**2/S) 
    132. 

  132. 

  133. !
    134. 

  134. 

  135. ! scalar fuel coefficients
    136. 

  136.    REAL, SAVE:: cmbcnst,hfgl,fuelmc_g,fuelmc_c, fire_wind_height
    137. 

  137. ! computed values
    138. 

  138.    REAL, SAVE:: fuelheat
    139. 

  139. 

  140. ! defaults, may be changed in init_fuel_cats
    141. 

  141.    DATA cmbcnst  / 17.433e+06/             ! J/kg dry fuel
    142. 

  142.    DATA hfgl     / 17.e4 /                ! W/m^2
    143. 

  143.    DATA fuelmc_g / 0.08  /                ! set = 0 for dry ground fuel
    144. 

  144.    DATA fuelmc_c / 1.00  /                ! set = 0 for dry canopy
    145. 

  145.    DATA fire_wind_height/6.096/           ! m, 6.096m Behave, 0 to use fwh in each category
    146. 

  146. !  REAL, PARAMETER :: bmst     = fuelmc_g/(1+fuelmc_g)
    147. 

  147. !  REAL, PARAMETER :: fuelheat = cmbcnst * 4.30e-04     ! convert J/kg to BTU/lb
    148. 

  148. !  real, parameter :: xlv      = 2.5e6                  ! to make it selfcontained
    149. 

  149. !  real, parameter :: cp      =  7.*287./2              ! to make it selfcontained
    150. 

  150. 

  151. 

  152. ! fuel categorytables
    153. 

  153.    INTEGER, PARAMETER :: nf=14              ! fuel cats in data stmts, for fillers only`
    154. 

  154.    INTEGER, SAVE      :: nfuelcats = 13     ! number of fuel categories, can be reset from namelist.fire
    155. 

  155.    INTEGER, PARAMETER :: zf = mfuelcats-nf  ! number of zero fillers in data stmt 
    156. 

  156.    INTEGER, SAVE      :: no_fuel_cat = 14   ! special no fuel category outside of 1:nfuelcats
    157. 

  157.    INTEGER, SAVE      :: no_fuel_cat2 = 14  ! all categories between no_fuel_cat and no_fuel_cat2 are no fuel
    158. 

  158.    INTEGER, SAVE      :: ibeh=1             ! type of spread formula
    159. 

  159.    CHARACTER (len=80), DIMENSION(mfuelcats ), save :: fuel_name
    160. 

  160.    INTEGER, DIMENSION( mfuelcats ), save :: ichap
    161. 

  161.    REAL   , DIMENSION( mfuelcats ), save :: windrf,weight,fgi,fci,fci_d,fct,fcbr, &
    162. 

  162.                                             fueldepthm,fueldens,fuelmce,   &
    163. 

  163.                                             fcwh,fcz0, ffw, &
    164. 

  164.                                             savr,st,se,adjr0,adjrw,adjrs, &
    165. 

  165.                                             fmc_gl_stdev,fmc_gl_ndwi_0,fmc_gl_ndwi_rate,fmc_gl_ndwi_stdev
    166. 

  166.    REAL   , DIMENSION( mfuelcats , max_moisture_classes), save :: fmc_gw
    167. 

  167. ! =============================================================================
    168. 

  168. ! Standard 13 fire behavior fuel models (for surface fires), along with some
    169. 

  169. !          estimated canopy properties (for crown fire).
    170. 

  170. ! =============================================================================
    171. 

  171.    DATA fuel_name /  &
    172. 

  172.      'FUEL MODEL 1: Short grass (1 ft)', &
    173. 

  173.      'FUEL MODEL 2: Timber (grass and understory)', &
    174. 

  174.      'FUEL MODEL 3: Tall grass (2.5 ft)', &
    175. 

  175.      'FUEL MODEL 4: Chaparral (6 ft)', &
    176. 

  176.      'FUEL MODEL 5: Brush (2 ft) ', &
    177. 

  177.      'FUEL MODEL 6: Dormant brush, hardwood slash', &
    178. 

  178.      'FUEL MODEL 7: Southern rough', &
    179. 

  179.      'FUEL MODEL 8: Closed timber litter', &
    180. 

  180.      'FUEL MODEL 9: Hardwood litter', &
    181. 

  181.      'FUEL MODEL 10: Timber (litter + understory)', &
    182. 

  182.      'FUEL MODEL 11: Light logging slash', &
    183. 

  183.      'FUEL MODEL 12: Medium logging slash', &
    184. 

  184.      'FUEL MODEL 13: Heavy logging slash', &
    185. 

  185.      'FUEL MODEL 14: no fuel', &
    186. 

  186.      zf*' '/
    187. 

  187.    DATA windrf /0.36, 0.36, 0.44,  0.55,  0.42,  0.44,  0.44, &     
    188. 

  188.                 0.36, 0.36, 0.36,  0.36,  0.43,  0.46,  1e-7, zf*0 / ! added jmandel October 2010
    189. 

  189.    DATA fgi / 0.166, 0.897, 0.675, 2.468, 0.785, 1.345, 1.092, &
    190. 

  190.               1.121, 0.780, 2.694, 2.582, 7.749, 13.024, 1.e-7, zf*0.  /
    191. 

  191.    DATA fueldepthm /0.305,  0.305,  0.762, 1.829, 0.61,  0.762,0.762, &
    192. 

  192.                     0.0610, 0.0610, 0.305, 0.305, 0.701, 0.914, 0.305,zf*0. /
    193. 

  193.    DATA savr / 3500., 2784., 1500., 1739., 1683., 1564., 1562.,  &
    194. 

  194.                1889., 2484., 1764., 1182., 1145., 1159., 3500., zf*0. /
    195. 

  195.    DATA fuelmce / 0.12, 0.15, 0.25, 0.20, 0.20, 0.25, 0.40,  &
    196. 

  196.                   0.30, 0.25, 0.25, 0.15, 0.20, 0.25, 0.12 , zf*0. / 
    197. 

  197.    DATA fueldens / nf * 32., zf*0. /   ! 32 if solid, 19 if rotten.
    198. 

  198.    DATA st / nf* 0.0555 , zf*0./
    199. 

  199.    DATA se / nf* 0.010 , zf*0./
    200. 

  200. ! ----- Notes on weight: (4) - best fit of Latham data;
    201. 

  201. !                 (5)-(7) could be 60-120; (8)-(10) could be 300-1600;
    202. 

  202. !                 (11)-(13) could be 300-1600
    203. 

  203.    DATA weight / 7.,  7.,  7., 180., 100., 100., 100.,  &
    204. 

  204.               900., 900., 900., 900., 900., 900., 7. , zf*0./ 
    205. 

  205. ! ----- 1.12083 is 5 tons/acre.  5-50 tons/acre orig., 100-300 after blowdown
    206. 

  206.    DATA fci_d / 0., 0., 0., 1.123, 0., 0., 0.,  &
    207. 

  207.             1.121, 1.121, 1.121, 1.121, 1.121, 1.121, 0., zf*0./
    208. 

  208.    DATA fct / 60., 60., 60., 60., 60., 60., 60.,  &
    209. 

  209.             60., 120., 180., 180., 180., 180. , 60. , zf*0.   /
    210. 

  210.    DATA ichap / 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , zf*0/ 
    211. 

  211.    DATA fcwh  / 6.096, 6.096, 6.096, 6.096, 6.096, 6.096, 6.096, &
    212. 

  212.          6.096, 6.096, 6.096, 6.096, 6.096, 6.096, 6.096, zf*0. / ! consistent with BEHAVE
    213. 

  213. ! roughness length 0.13*fueldepthm except cat 3 fz0=0.1 for consistency with landuse
    214. 

  214.  ! fz0 =   0.0396,0.0396,0.0991,0.2378,0.0793,0.0991,0.0991,
    215. 

  215.    DATA fcz0 /  0.0396,0.0396,0.1000,0.2378,0.0793,0.0991,0.0991, &
    216. 

  216.          0.0079,0.0079,0.0396,0.0396,0.0911,0.1188,0.0396, zf * 0. /
    217. 

  217.    !DATA fcwh /0.6 , 0.6,  1.5, 36,  1.2,  1.5, 1.5, &
    218. 

  218.    !           0.12, 0.12, 0.6, 0.6, 1.38, 1.8, 1.8, zf*0 / ! fuel wind height, added jm 2/23/11
    219. 

  219.    !DATA fcz0 /0.3, 0.3, 0.75, 18.,  0.6,  0.75, 0.75, &
    220. 

  220.    !           0.06, 0.06, 0.3, 0.3, 0.69, 0.9, 0.9, zf*0 /  ! fuel roughness height, added jm 2/23/11
    221. 

  221.    DATA ffw /nf* 0.9, zf*0/
    222. 

  222.    DATA fmc_gl_ndwi_0 /nf*0.1, zf*0./
    223. 

  223.    DATA fmc_gl_ndwi_rate /nf*0.6, zf*0./
    224. 

  224.    DATA fmc_gl_ndwi_stdev /nf*0.2, zf*0./
    225. 

  225.    DATA fmc_gl_stdev /nf*0.2, zf*0./
    226. 

  226.    DATA adjr0 /mfuelcats*1./
    227. 

  227.    DATA adjrw /mfuelcats*1./
    228. 

  228.    DATA adjrs /mfuelcats*1./
    229. 

  229.    
    230. 

  230.    ! =========================================================================
    231. 

    232. 

  232.   logical, save :: have_fuel_cats=.false.
    233. 

    234. 

  234. contains
    235. 

    236. 

  236. subroutine fuel_moisture(                &
    237. 

  237.     id,                                  & ! for prints and maybe file names
    238. 

  238.     nfmc,                                &
    239. 

  239.     ids,ide, jds,jde,                    & ! atm grid dimensions
    240. 

  240.     ims,ime, jms,jme,                    &
    241. 

  241.     ips,ipe,jps,jpe,                     &
    242. 

  242.     its,ite,jts,jte,                     &
    243. 

  243.     ifds, ifde, jfds, jfde,              & ! fire grid dimensions
    244. 

  244.     ifms, ifme, jfms, jfme,              &
    245. 

  245.     ifts,ifte,jfts,jfte,                 &
    246. 

  246.     ir,jr,                               & ! atm/fire grid ratio
    247. 

  247.     nfuel_cat,                           & ! fuel data
    248. 

  248.     fndwi,                                & ! satellite sensing on fire grid
    249. 

  249.     fmc_gc,                              & ! moisture contents by class on atmospheric grid
    250. 

  250.     fmc_g                                & ! weighted fuel moisture contents on fire grid
    251. 

  251.     )
    252. 

    253. 

  253. implicit none
    254. 

    255. 

  255. !**** arguments
    256. 

  256. integer, intent(in)::                    &
    257. 

  257.     id,nfmc,                             &
    258. 

  258.     ids,ide, jds,jde,                    & ! atm grid dimensions
    259. 

  259.     ims,ime, jms,jme,                    &
    260. 

  260.     ips,ipe,jps,jpe,                     &
    261. 

  261.     its,ite,jts,jte,                     &
    262. 

  262.     ifds, ifde, jfds, jfde,              & ! fire grid dimensions
    263. 

  263.     ifms, ifme, jfms, jfme,              &
    264. 

  264.     ifts,ifte,jfts,jfte,                 &
    265. 

  265.     ir,jr                                  ! atm/fire grid ratio
    266. 

    267. 

    268. 

  268. real,intent(in),dimension(ifms:ifme,jfms:jfme):: nfuel_cat, & ! fuel data
    269. 

  269.                                        fndwi                  ! satellite sensing interpolated to fire grid
    270. 

  270. real,intent(inout),dimension(ims:ime,nfmc,jms:jme):: fmc_gc
    271. 

  271. real,intent(out),dimension(ifms:ifme,jfms:jfme):: fmc_g ! fuel data
    272. 

    273. 

  273. !**** local
    274. 

  274. real, dimension(its-1:ite+1,jts-1:jte+1):: fmc_k  ! copy of fmc_gc(:,k,:) 
    275. 

  275. real, dimension(ifts:ifte,jfts:jfte):: fmc_f,  &   ! interpolation of fmc_gc(:,k,:) to the fire grid 
    276. 

  276.                                        nwdi_f      ! inerpolation of nwdi to the fire grid
    277. 

  277. integer::i,j,k,n
    278. 

  278. integer::ibs,ibe,jbs,jbe
    279. 

  279. real::f1,w1,w2,f2,fa,fc
    280. 

    281. 

  281. character(len=128)::msg
    282. 

    283. 

  283. call check_mesh_2dim(ifts,ifte,jfts,jfte,ifds,ifde,jfds,jfde) ! check if fire tile fits into domain
    284. 

  284. call check_mesh_2dim(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme) ! check if fire tile fits into memory
    285. 

    286. 

  286. do j=jfts,jfte
    287. 

  287.     do i=ifts,ifte
    288. 

  288.         fmc_g(i,j)=0.               ! initialize sum over classes
    289. 

  289.     enddo
    290. 

  290. enddo
    291. 

    292. 

  292. ! one beyond the tile but not beyond the domain boundary
    293. 

  293. ibs=max(ids,its-1)
    294. 

  294. ibe=min(ide,ite+1)
    295. 

  295. jbs=max(jds,jts-1)
    296. 

  296. jbe=min(jde,jte+1)
    297. 

    298. 

  298. call check_mesh_2dim(ibs,ibe,jbs,jbe,ims,ime,jms,jme) ! check if tile with halo fits into memory
    299. 

    300. 

  300. do k=1,moisture_classes
    301. 

    302. 

  302.     ! copy halo beyond the tile but not beyond the domain boundary
    303. 

  303.     do j=jbs,jbe
    304. 

  304.         do i=ibs,ibe
    305. 

  305.             fmc_k(i,j)=fmc_gc(i,k,j)      ! copy slice to 2d array
    306. 

  306.         enddo
    307. 

  307.     enddo
    308. 

    309. 

  309.     call print_2d_stats(ibs,ibe,jbs,jbe,its-1,ite+1,jts-1,jte+1,fmc_k,'fuel_moisture: fmc_k')
    310. 

    311. 

  311.     ! interpolate moisture contents in the class k to the fire mesh
    312. 

  312.     call interpolate_z2fire(id,0,        & ! for debug output, <= 0 no output
    313. 

  313.     ids,ide,jds,jde,                     & ! atm grid dimensions
    314. 

  314.     its-1,ite+1,jts-1,jte+1,             & ! memory dimensions
    315. 

  315.     ips,ipe,jps,jpe,                     &
    316. 

  316.     its,ite,jts,jte,                     &
    317. 

  317.     ifds, ifde, jfds, jfde,              & ! fire grid dimensions
    318. 

  318.     ifts, ifte, jfts, jfte,              &
    319. 

  319.     ifts,ifte,  jfts,jfte,               &
    320. 

  320.     ir,jr,                               & ! atm/fire grid ratio
    321. 

  321.     fmc_k,                               & ! atm grid arrays in
    322. 

  322.     fmc_f)                                 ! fire grid arrays out
    323. 

    324. 

  324.     call print_2d_stats(ifts,ifte,jfts,jfte,ifts,ifte,jfts,jfte,fmc_f,'fuel_moisture: fmc_f')
    325. 

    326. 

  326.     if(k .eq. kfmc_ndwi)then  ! if live moisture, assimilate ndwi
    327. 

  327.         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fndwi,'fuel_moisture: fndwi')
    328. 

  328.         write(msg,'(a,i4)')'Assimilating NDWI in fuel moisture class ',k
    329. 

  329.         call message(msg)
    330. 

  330.     endif
    331. 

    332. 

  332.     ! add moisture contents for class k to the fuel moisture
    333. 

  333.     do j=jfts,jfte
    334. 

  334.         do i=ifts,ifte
    335. 

  335.             n = nfuel_cat(i,j)
    336. 

  336.             if(n > 0)then
    337. 

  337.                 if(k .ne. kfmc_ndwi)then 
    338. 

  338.                    fmc_g(i,j)=fmc_g(i,j)+fmc_gw(n,k)*fmc_f(i,j)      ! add to sum over classes
    339. 

  339.                 else                                                    ! if live moisture, assimilate ndwi
    340. 

  340.                    f1=fmc_f(i,j)
    341. 

  341.                    w1 = fmc_gl_stdev(n)                            
    342. 

  342.                    w1 = 1./(w1*w1)                                      ! weight of forecast
    343. 

  343.                    w2 = fmc_gl_ndwi_stdev(n)
    344. 

  344.                    w2 = 1./(w2*w2)                                      ! weight of update
    345. 

  345.                    f2 = fmc_gl_ndwi_0(n) + fmc_gl_ndwi_rate(n) * fndwi(i,j)   ! from regression
    346. 

  346.                    fa = (w1*f1 + w2*f2) / (w1 + w2)                     ! updated value
    347. 

  347.                    fc = fmc_gw(n,k)*fa                                  ! times proportion of live fuel
    348. 

  348.                    fmc_g(i,j)=fmc_g(i,j)+fc      ! add to sum over classes
    349. 

  349.                    ! write(*,*)'NDWI:',i,j,f1,f2,w1,w2,f1,fa,fmc_gw(n,k),fc,fmc_g(i,j)
    350. 

  350.                 endif
    351. 

  351.             endif
    352. 

  352.         enddo
    353. 

  353.     enddo
    354. 

    355. 

  355. enddo
    356. 

    357. 

  357. call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fmc_g,'fuel_moisture: fmc_g')
    358. 

    359. 

    360. 

  360. end subroutine fuel_moisture
    361. 

    362. 

  362. subroutine advance_moisture(    &
    363. 

  363.     initialize,                 & ! initialize timestepping. true on the first call at time 0, then false
    364. 

  364.     ims,ime,  jms,jme,          & ! memory dimensions
    365. 

  365.     its,ite,  jts,jte,          & ! tile dimensions
    366. 

  366.     nfmc,                       & ! dimension of moisture fields
    367. 

  367.     moisture_dt,                & ! timestep = time step time elapsed from the last call
    368. 

  368.     fmep_decay_tlag,            & ! moisture extended model assimilated diffs. decay time lag
    369. 

  369.     rainc, rainnc,              & ! accumulated rain 
    370. 

  370.     t2, q2, psfc,               & ! temperature (K), vapor contents (kg/kg), pressure (Pa) at the surface
    371. 

  371.     rain_old,                   & ! previous value of accumulated rain
    372. 

  372.     t2_old, q2_old, psfc_old,   & ! previous values of the atmospheric state at surface
    373. 

  373.     rh_fire,                    & ! relative humidity at surface, for diagnostic only
    374. 

  374.     fmc_gc,                     & ! fuel moisture by class, updated
    375. 

  375.     fmep,                       & ! fuel moisture extended model parameters
    376. 

  376.     fmc_equi,                   & ! fuel moisture equilibrium by class, for diagnostics only
    377. 

  377.     fmc_lag                     & ! fuel moisture tendency by classe, for diagnostics only
    378. 

  378.     )
    379. 

    380. 

  380. implicit none
    381. 

    382. 

  382. !*** arguments
    383. 

  383. logical, intent(in):: initialize
    384. 

  384. integer, intent(in)::           &
    385. 

  385.     ims,ime,  jms,jme,          & ! memory dimensions
    386. 

  386.     its,ite,  jts,jte,          & ! tile dimensions
    387. 

  387.     nfmc                          ! number of moisture fields
    388. 

  388. real, intent(in):: moisture_dt, fmep_decay_tlag
    389. 

  389. real, intent(in), dimension(ims:ime,jms:jme):: t2, q2, psfc, rainc, rainnc
    390. 

  390. real, intent(inout), dimension(ims:ime,jms:jme):: t2_old, q2_old, psfc_old, rain_old 
    391. 

  391. real, intent(inout), dimension(ims:ime,nfmc,jms:jme):: fmc_gc
    392. 

  392. real, intent(inout), dimension(ims:ime,2,jms:jme):: fmep
    393. 

  393. real, intent(out), dimension(ims:ime,nfmc,jms:jme):: fmc_equi, fmc_lag
    394. 

  394. real, intent(out), dimension(ims:ime,jms:jme)::rh_fire 
    395. 

    396. 

  396. !*** global
    397. 

  397. ! fuel properties moisture set by init_fuel_cats
    398. 

    399. 

  399. !*** local
    400. 

  400. integer:: i,j,k
    401. 

  401. real::rain_int, T, P, Q, QRS, ES, RH, tend, EMC_d, EMC_w, EMC, R, rain_diff, fmc, rlag, equi, &
    402. 

  402.     d, w, rhmax, rhmin, change, rainmax,rainmin, fmc_old, H, deltaS, deltaE
    403. 

  403. real, parameter::tol=1e-2 ! relative change larger than that will switch to exponential ode solver 
    404. 

  404. character(len=256)::msg
    405. 

  405. logical::bad_wrf
    406. 

  406. integer::msglevel=2
    407. 

  407. logical, parameter::check_rh=.false.
    408. 

  408. integer::check_data=2  ! 0=nothing, 1=replace quietly, 2=warning (also printed if msglevel>2), 3=crash
    409. 

  409. real::epsilon,Pws,Pw,t2_min,q2_min,psfc_min
    410. 

  410. real::t2_floor=200.    ! minimum allowed temperature (K)
    411. 

  411. real::q2_floor=1e-8    ! miniumu allowed moisture contents (kg/kg)
    412. 

  412. real::psfc_floor=1000. ! mimimum allowed surface pressiure (Pa)
    413. 

    414. 

    415. 

  415. !*** executable
    416. 

    417. 

  417. ! check arguments
    418. 

    419. 

  419. if(msglevel>1)then
    420. 

  420. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    421. 

  421.     write(msg,'(a,f10.2,a,i4,a,i4)')'advance_moisture dt=',moisture_dt,'s using ',moisture_classes,' classes from possible ',nfmc
    422. 

  422. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    423. 

  423.     call message(msg,level=2)
    424. 

  424. endif
    425. 

    426. 

  426. if(moisture_classes > nfmc .or. moisture_classes > max_moisture_classes)then
    427. 

  427. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    428. 

  428.     write(msg,*)'advance_moisture: moisture_classes=',moisture_classes, &
    429. 

  429.        ' > nfmc=',nfmc,' or >  max_moisture_classes=',max_moisture_classes
    430. 

  430. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    431. 

  431.     call crash(msg)
    432. 

  432. endif
    433. 

    434. 

  434. call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,t2,'T2')
    435. 

  435. call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,q2,'Q2')
    436. 

  436. call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,psfc,'PSFC')
    437. 

    438. 

  438. if(initialize) then 
    439. 

  439.     call message('advance_moisture: initializing, copying surface variables to old')
    440. 

  440.     call copy2old
    441. 

  441. else
    442. 

  442.     call print_3d_stats_by_slice(its,ite,1,moisture_classes,jts,jte,ims,ime,1,nfmc,jms,jme,fmc_gc,'before advance fmc_gc')
    443. 

  443. endif
    444. 

    445. 

  445. if(check_data.ge.2 .or. msglevel.ge.2)then
    446. 

  446.     t2_min = huge(t2_min)
    447. 

  447.     q2_min = huge(q2_min)
    448. 

  448.     psfc_min = huge(psfc_min)
    449. 

  449.     do j=jts,jte
    450. 

  450.         do i=its,ite
    451. 

  451.             t2_min=min(t2(i,j),t2_min)
    452. 

  452.             q2_min=min(q2(i,j),q2_min)
    453. 

  453.             psfc_min=min(psfc(i,j),psfc_min)
    454. 

  454.         enddo
    455. 

  455.     enddo
    456. 

  456.     bad_wrf = ( t2_min<t2_floor .or. psfc_min<psfc_floor .or. q2_min<q2_floor)
    457. 

  457.     if(bad_wrf .or. msglevel.ge.2)then
    458. 

  458. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    459. 

  459. 91      format(a,3(2x,a,e11.3))
    460. 

  460.         write(msg,91)'minimal  ','t2',t2_min,'q2',q2_min,'psfc',psfc_min
    461. 

  461.         call message(msg,level=0) 
    462. 

  462.         write(msg,91)'floor    ','t2',t2_floor,'q2',q2_floor,'psfc',psfc_floor
    463. 

  463.         call message(msg,level=0) 
    464. 

  464. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    465. 

  465.     endif
    466. 

  466.     if(bad_wrf)then
    467. 

  467.         if(check_data.ge.3)then
    468. 

  468.             call crash('advance_moisture: invalid data passed from WRF')
    469. 

  469.         else
    470. 

  470.             call message('WARNING: advance_moisture: nonphysical input values replaced by floor',level=0)
    471. 

  471.         endif
    472. 

  472.     endif
    473. 

  473. endif
    474. 

    475. 

  475. ! one time step
    476. 

    477. 

  477. rhmax=-huge(rhmax)
    478. 

  478. rhmin=huge(rhmin)
    479. 

  479. rainmax=-huge(rainmax)
    480. 

  480. rainmin= huge(rainmin)
    481. 

  481. do j=jts,jte
    482. 

  482.     do k=1,moisture_classes
    483. 

  483.         do i=its,ite
    484. 

  484.             ! old fuel moisture contents
    485. 

  485.             ! compute the rain intensity from the difference of accumulated rain
    486. 

  486.             rain_diff = ((rainc(i,j) + rainnc(i,j)) - rain_old(i,j))
    487. 

  487.             if(moisture_dt > 0.)then
    488. 

  488.                 rain_int  = 3600. * rain_diff / moisture_dt 
    489. 

  489.             else
    490. 

  490.                 rain_int  = 0.
    491. 

  491.             endif
    492. 

  492.             rainmax = max(rainmax,rain_int)
    493. 

  493.             rainmin = min(rainmin,rain_int)
    494. 

  494.             R = rain_int - rain_threshold(k)
    495. 

    496. 

  496.             ! average the inputs for second order accuracy
    497. 

  497.             T = 0.5 * (t2_old(i,j) + t2(i,j))
    498. 

  498.             P = 0.5 * (psfc_old(i,j) + psfc(i,j))
    499. 

  499.             Q = 0.5 * (q2_old(i,j) + q2(i,j))
    500. 

    501. 

  501.             ! replace nonphysical values by floor
    502. 

  502.             if(check_data .ge. 1)then
    503. 

  503.                 T = max(T,t2_floor)
    504. 

  504.                 P = max(P,psfc_floor)
    505. 

  505.                 Q = max(Q,q2_floor)
    506. 

  506. 	    endif
    507. 

    508. 

  508.             ! compute the relative humidity
    509. 

  509.             ! ES=610.78*exp(17.269*(T-273.161)/(T-35.861))
    510. 

  510.             ! QRS=0.622*ES/(P-0.378*ES)
    511. 

  511.             ! RH = Q/QRS
    512. 

  512.             ! function rh_from_q from Adam Kochanski following Murphy and Koop, Q.J.R. Meteorol. Soc (2005) 131 1539-1565 eq. (10)
    513. 

  513.             epsilon = 0.622 ! Molecular weight of water (18.02 g/mol) to molecular weight of dry air (28.97 g/mol)
    514. 

  514.             ! vapor pressure [Pa]
    515. 

  515.             Pw=q*P/(epsilon+(1-epsilon)*q); 
    516. 

  516.             ! saturation vapor pressure [Pa]
    517. 

  517.             Pws= exp( 54.842763 - 6763.22/T - 4.210 * log(T) + 0.000367*T + &
    518. 

  518.                 tanh(0.0415*(T - 218.8)) * (53.878 - 1331.22/T - 9.44523 * log(T) + 0.014025*T))
    519. 

  519.             !realtive humidity [1]
    520. 

  520.             RH = Pw/Pws
    521. 

  521.             rh_fire(i,j)=RH
    522. 

  522.             rhmax=max(RH,rhmax)
    523. 

  523.             rhmin=min(RH,rhmin)
    524. 

    525. 

  525.             deltaE = fmep(i,1,j)
    526. 

  526.             deltaS = fmep(i,2,j)
    527. 

    528. 

  528.             if(.not.check_rh)then
    529. 

  529.                 RH = min(RH,1.0)
    530. 

  530.             else
    531. 

  531.                 if(RH < 0.0 .or. RH > 1.0 .or. RH .ne. RH )then
    532. 

  532. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    533. 

  533.                     write(msg,'(a,2i6,5(a,f10.2))')'At i,j ',i,j,' RH=',RH, &
    534. 

  534.                         ' from T=',T,' P=',P,' Q=',Q
    535. 

  535.                     call message(msg) 
    536. 

  536.                     call crash('Relative humidity must be between 0 and 1, saturated water contents must be >0')
    537. 

  537. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    538. 

  538.                 endif
    539. 

  539.             endif 
    540. 

  540.             !print *,'ADV_MOIST i=',i,' j=',j,' T=',T,' P=',P,' Q=',Q,' ES=',ES,' QRS=',QRS,' RH=',RH
    541. 

    542. 

  542.             if (R > 0.) then
    543. 

  543.                 select case(wetting_model(k))
    544. 

  544.                 case(1) ! saturation_moisture=2.5 wetting_lag=14h saturation_rain=8 mm/h calibrated to VanWagner&Pickett 1985 per 24 hours
    545. 

  545.                     EMC_w=saturation_moisture(k) + deltaS
    546. 

  546.                     EMC_d=saturation_moisture(k) + deltaS
    547. 

  547.                     rlag=rec_wetting_lag_sec(k) * (1. - exp(-R/saturation_rain(k)))
    548. 

  548.                 end select
    549. 

  549.             else ! not raining
    550. 

  550.                 select case(drying_model(k))
    551. 

  551.                 case(1) ! Van Wagner and Pickett (1972, 1985) per Viney (1991) eq (7) (8)
    552. 

  552.                     H = RH * 100.
    553. 

  553.                     d=0.942*H**0.679 + 0.4994e-4*exp(0.1*H) + 0.18*(21.1+273.15-T)*(1-exp(-0.115*H)) ! equilibrium moisture for drying
    554. 

  554.                     w=0.618*H**0.753 + 0.4540e-4*exp(0.1*H) + 0.18*(21.1+273.15-T)*(1-exp(-0.115*H)) ! equilibrium moisture for adsorbtion
    555. 

  555.                     if(d.ne.d.or.w.ne.w)call crash('equilibrium moisture calculation failed, result is NaN')
    556. 

  556.                     d = d*0.01
    557. 

  557.                     w = w*0.01
    558. 

  558.                     EMC_d = max(max(d,w)+deltaE,0.0)
    559. 

  559.                     EMC_w = max(min(d,w)+deltaE,0.0)
    560. 

  560.                     rlag=rec_drying_lag_sec(k)
    561. 

  561.                 end select
    562. 

  562.             endif
    563. 

  563.             !*** MODELS THAT ARE NOT OF THE EXPONENTIAL TIME LAG KIND 
    564. 

  564.             ! ARE RESPONSIBLE FOR THEIR OWN LOGIC, THESE MODELS
    565. 

  565.             ! SHOULD COMPUTE fmc_gc(i,k,j) DIRECTLY AND SET TLAG < 0
    566. 

  566.             !
    567. 

  567.             if(rlag > 0.0)then
    568. 

    569. 

  569.                 if(.not.initialize .or. fmc_gc_initialization(k).eq.0)then ! take old from before, no initialization
    570. 

  570.                     fmc_old = fmc_gc(i,k,j)
    571. 

  571.                 elseif(fmc_gc_initialization(k).eq.1)then ! from scalar fuelmc_g
    572. 

  572.                     fmc_old = fuelmc_g
    573. 

  573.                 elseif(fmc_gc_initialization(k).eq.2)then ! from computed equilibrium
    574. 

  574.                     fmc_old=0.5*(EMC_d+EMC_w)
    575. 

  575.                 else
    576. 

  576.                     call crash('bad value of fmc_gc_initialization(k), must be between 0 and 2')
    577. 

  577.                 endif
    578. 

  578.                 equi = max(min(fmc_old, EMC_d),EMC_w) ! take lower or upper equilibrium value 
    579. 

    580. 

  580.                 change = moisture_dt * rlag
    581. 

    582. 

  582.                 if(change  < tol)then
    583. 

  583.                      if(fire_print_msg.ge.3)call message('midpoint method')
    584. 

  584.                      fmc = fmc_old + (equi - fmc_old)*change*(1.0 - 0.5*change)  ! 2nd order Taylor
    585. 

  585.                 else
    586. 

  586.                      if(fire_print_msg.ge.3)call message('exponential method')
    587. 

  587.                      fmc = fmc_old + (equi - fmc_old)*(1 - exp(-change))
    588. 

  588.                 endif
    589. 

  589.                 fmc_gc(i,k,j) = fmc
    590. 

    591. 

  591.                 ! diagnostics out
    592. 

  592.                 fmc_equi(i,k,j)=equi
    593. 

  593.                 fmc_lag(i,k,j)=1.0/(3600.0*rlag)
    594. 

    595. 

  595.                 ! diagnostic prints
    596. 

  596.                 if(fire_print_msg.ge.3)then
    597. 

  597. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    598. 

  598. 		    write(msg,*)'i=',i,' j=',j,'EMC_w=',EMC_w,' EMC_d=',EMC_d
    599. 

  599.                     call message(msg)
    600. 

  600.                     write(msg,*)'fmc_old=',fmc,' equi=',equi,' change=',change,' fmc=',fmc
    601. 

  601.                     call message(msg)
    602. 

  602. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    603. 

  603.                 endif
    604. 

    605. 

  605.             endif
    606. 

  606.       enddo
    607. 

  607.    enddo
    608. 

  608. enddo
    609. 

    610. 

    611. 

  611. ! assimilated differences decay
    612. 

  612. do j=jts,jte
    613. 

  613.   do k=1,2
    614. 

  614.     do i=its,ite
    615. 

  615.       change = moisture_dt / (fmep_decay_tlag * 3600.)
    616. 

  616.       if(change < tol) then
    617. 

  617.         fmep(i,k,j) = fmep(i,k,j)*(1.0 - change * (1.0 - 0.5 * change))
    618. 

  618.       else
    619. 

  619.         fmep(i,k,j) = fmep(i,k,j)*exp(-change)
    620. 

  620.       endif
    621. 

  621.     enddo
    622. 

  622.   enddo
    623. 

  623. enddo
    624. 

    625. 

    626. 

  626. if(fire_print_msg.ge.2)then
    627. 

  627. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    628. 

  628.     write(msg,2)'Rain intensity    min',rainmin,  ' max',rainmax,' mm/h'
    629. 

  629.     call message(msg) 
    630. 

  630.     if(rainmin <0.)then
    631. 

  631.         call message('WARNING rain accumulation must increase')
    632. 

  632.     endif
    633. 

  633.     write(msg,2)'Relative humidity min',100*rhmin,' max',100*rhmax,'%'
    634. 

  634.     call message(msg) 
    635. 

  635.     if(.not.(rhmax<=1.0 .and. rhmin>=0))then
    636. 

  636.         call message('WARNING Relative humidity must be between 0 and 100%')
    637. 

  637.     endif
    638. 

  638. 2   format(2(a,f10.2),a)
    639. 

  639. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    640. 

  640. endif
    641. 

    642. 

  642. call print_3d_stats_by_slice(its,ite,1,moisture_classes,jts,jte,ims,ime,1,nfmc,jms,jme,fmc_equi,'equilibrium fmc_equi')
    643. 

  643. call print_3d_stats_by_slice(its,ite,1,moisture_classes,jts,jte,ims,ime,1,nfmc,jms,jme,fmc_lag,'time lag')
    644. 

  644. call print_3d_stats_by_slice(its,ite,1,moisture_classes,jts,jte,ims,ime,1,nfmc,jms,jme,fmc_gc,'after advance fmc_gc')
    645. 

    646. 

  646. call copy2old
    647. 

    648. 

  648. return
    649. 

    650. 

  650. contains
    651. 

    652. 

  652. subroutine copy2old
    653. 

    654. 

  654. do j=jts,jte
    655. 

  655.     do i=its,ite
    656. 

  656.         rain_old(i,j) = rainc(i,j) + rainnc(i,j)
    657. 

  657.         t2_old(i,j) = t2(i,j)
    658. 

  658.         q2_old(i,j) = q2(i,j)
    659. 

  659.         psfc_old(i,j) = psfc(i,j)
    660. 

  660.     enddo
    661. 

  661. enddo
    662. 

    663. 

  663. end subroutine copy2old
    664. 

    665. 

  665. subroutine get_equi_moist
    666. 

  666. end subroutine get_equi_moist
    667. 

    668. 

  668. end subroutine advance_moisture
    669. 

    670. 

    671. 

  671. subroutine init_fuel_cats(init_fuel_moisture)
    672. 

  672. implicit none
    673. 

  673. !*** purpose: initialize fuel tables and variables by constants
    674. 

  674. !*** arguments: 
    675. 

  675. logical, intent(in)::init_fuel_moisture
    676. 

  676. logical, external:: wrf_dm_on_monitor
    677. 

  677. !$ integer, external:: OMP_GET_THREAD_NUM
    678. 

  678. !*** local
    679. 

  679. integer:: i,j,k,ii,iounit,ierr,kk
    680. 

  680. character(len=128):: msg
    681. 

  681. REAL   , DIMENSION( mfuelcats ) :: fwh, fz0
    682. 

  682. !*** executable
    683. 

    684. 

  684. ! read 
    685. 

  685. namelist /fuel_scalars/ cmbcnst,hfgl,fuelmc_g,fuelmc_c,nfuelcats,no_fuel_cat,no_fuel_cat2,fire_wind_height,ibeh
    686. 

  686. namelist /fuel_categories/ fuel_name,windrf,fgi,fueldepthm,savr, &
    687. 

  687.     fuelmce,fueldens,st,se,weight,fci_d,fct,ichap,fwh,fz0,ffw, &
    688. 

  688.     fmc_gl_ndwi_0,fmc_gl_ndwi_rate,fmc_gl_ndwi_stdev, fmc_gl_stdev, &
    689. 

  689.     adjr0,adjrw,adjrs,fmc_gw01,fmc_gw02,fmc_gw03,fmc_gw04,fmc_gw05
    690. 

    691. 

  691. namelist /moisture/ moisture_classes,drying_lag,wetting_lag,saturation_moisture,saturation_rain,rain_threshold, &
    692. 

  692.     drying_model,wetting_model, moisture_class_name,fmc_gc_initialization
    693. 

    694. 

  694. !$  if (OMP_GET_THREAD_NUM() .ne. 0)then
    695. 

  695. !$     call crash('init_fuel_cats: must be called from master thread')
    696. 

  696. !$  endif
    697. 

  697.  
    698. 

  698. IF ( wrf_dm_on_monitor() ) THEN
    699. 

  699.     ! we are the master task
    700. 

  700.     
    701. 

  701.     ! copy in defaults
    702. 

  702.     fwh=fcwh
    703. 

  703.     fz0=fcz0
    704. 

  704.     
    705. 

  705.     ! read the file
    706. 

  706.     iounit=open_text_file('namelist.fire','read')
    707. 

  707.     read(iounit,fuel_scalars,iostat=ierr)
    708. 

  708.     if(ierr.ne.0)call crash('init_fuel_cats: error reading namelist fuel_scalars in file namelist.fire')
    709. 

  709.     read(iounit,fuel_categories,iostat=ierr)
    710. 

  710.     if(ierr.ne.0)call crash('init_fuel_cats: error reading namelist fuel_categories in file namelist.fire')
    711. 

  711.     if(init_fuel_moisture)then
    712. 

  712.         read(iounit,moisture,iostat=ierr)
    713. 

  713.         if(ierr.ne.0)call crash('init_fuel_cats: error reading namelist moisture in file namelist.fire')
    714. 

  714.     endif
    715. 

  715.     fmc_gw(1:mfuelcats,1)=fmc_gw01
    716. 

  716.     fmc_gw(1:mfuelcats,2)=fmc_gw02
    717. 

  717.     fmc_gw(1:mfuelcats,3)=fmc_gw03
    718. 

  718.     fmc_gw(1:mfuelcats,4)=fmc_gw04
    719. 

  719.     fmc_gw(1:mfuelcats,5)=fmc_gw05
    720. 

  720.     CLOSE(iounit)
    721. 

  721.     
    722. 

  722.     ! copy out to permanent names
    723. 

  723.     fcwh=fwh
    724. 

  724.     fcz0=fz0
    725. 

  725.     
    726. 

  726.     if (nfuelcats>mfuelcats) then
    727. 

  727.         write(msg,*)'nfuelcats=',nfuelcats,' too large, increase mfuelcats'
    728. 

  728.         call crash(msg)
    729. 

  729.     endif
    730. 

  730.     if (no_fuel_cat >= 1 .and. no_fuel_cat <= nfuelcats)then
    731. 

  731.         write(msg,*)'no_fuel_cat=',no_fuel_cat,' may not be between 1 and nfuelcats=',nfuelcats
    732. 

  732.         call crash(msg)
    733. 

  733.     endif
    734. 

  734.     if (no_fuel_cat > no_fuel_cat2)then
    735. 

  735.         write(msg,*)'no_fuel_cat=',no_fuel_cat,' must not be larger than no_fuel_cat2=',no_fuel_cat2
    736. 

  736.         call crash(msg)
    737. 

  737.     endif
    738. 

  738. ENDIF
    739. 

    740. 

  740. ! broadcast fuel tables
    741. 

  741. call wrf_dm_bcast_real(cmbcnst,1)
    742. 

  742. call wrf_dm_bcast_real(hfgl,1)
    743. 

  743. call wrf_dm_bcast_real(fuelmc_g,1)
    744. 

  744. call wrf_dm_bcast_real(fuelmc_c,1)
    745. 

  745. call wrf_dm_bcast_real(fire_wind_height,1)
    746. 

  746. call wrf_dm_bcast_integer(nfuelcats,1)
    747. 

  747. call wrf_dm_bcast_integer(no_fuel_cat,1)
    748. 

  748. call wrf_dm_bcast_integer(no_fuel_cat2,1)
    749. 

  749. call wrf_dm_bcast_integer(ibeh,1)
    750. 

  750. call wrf_dm_bcast_real(windrf,    nfuelcats)
    751. 

  751. call wrf_dm_bcast_real(fgi,       nfuelcats)
    752. 

  752. call wrf_dm_bcast_real(fueldepthm,nfuelcats)
    753. 

  753. call wrf_dm_bcast_real(savr,      nfuelcats)
    754. 

  754. call wrf_dm_bcast_real(fuelmce,   nfuelcats)
    755. 

  755. call wrf_dm_bcast_real(fueldens,  nfuelcats)
    756. 

  756. call wrf_dm_bcast_real(st,        nfuelcats)
    757. 

  757. call wrf_dm_bcast_real(se,        nfuelcats)
    758. 

  758. call wrf_dm_bcast_real(weight,    nfuelcats)
    759. 

  759. call wrf_dm_bcast_real(fci_d,     nfuelcats)
    760. 

  760. call wrf_dm_bcast_real(fct,       nfuelcats)
    761. 

  761. call wrf_dm_bcast_integer(ichap,  nfuelcats)
    762. 

  762. call wrf_dm_bcast_real(fcwh,      nfuelcats)
    763. 

  763. call wrf_dm_bcast_real(fcz0,      nfuelcats)
    764. 

  764. call wrf_dm_bcast_real(ffw,       nfuelcats)
    765. 

  765. call wrf_dm_bcast_real(adjr0,     nfuelcats)
    766. 

  766. call wrf_dm_bcast_real(adjrw,     nfuelcats)
    767. 

  767. call wrf_dm_bcast_real(adjrs,     nfuelcats)
    768. 

  768. call wrf_dm_bcast_real(fmc_gl_ndwi_0,    nfuelcats)
    769. 

  769. call wrf_dm_bcast_real(fmc_gl_ndwi_rate, nfuelcats)
    770. 

  770. call wrf_dm_bcast_real(fmc_gl_ndwi_stdev,nfuelcats)
    771. 

  771. call wrf_dm_bcast_real(fmc_gl_stdev,     nfuelcats)
    772. 

  772. ! broadcast moisture tables
    773. 

  773. call wrf_dm_bcast_integer(moisture_classes,1)
    774. 

  774. call wrf_dm_bcast_real(drying_lag,     max_moisture_classes)
    775. 

  775. call wrf_dm_bcast_real(wetting_lag,     max_moisture_classes)
    776. 

  776. call wrf_dm_bcast_real(saturation_moisture,     max_moisture_classes)
    777. 

  777. call wrf_dm_bcast_real(saturation_rain,     max_moisture_classes)
    778. 

  778. call wrf_dm_bcast_real(rain_threshold,     max_moisture_classes)
    779. 

  779. call wrf_dm_bcast_integer(drying_model,     max_moisture_classes)
    780. 

  780. call wrf_dm_bcast_integer(wetting_model,     max_moisture_classes)
    781. 

  781. call wrf_dm_bcast_integer(fmc_gc_initialization,     max_moisture_classes)
    782. 

  782. call wrf_dm_bcast_real(fmc_gw,     mfuelcats*max_moisture_classes)
    783. 

    784. 

  784. ! moisture model derived scalars
    785. 

  785. do i=1,moisture_classes
    786. 

  786.     rec_drying_lag_sec(i)  = 1.0/(3600.0*drying_lag(i))
    787. 

  787.     rec_wetting_lag_sec(i) = 1.0/(3600.0*wetting_lag(i))
    788. 

  788. enddo
    789. 

    790. 

  790. !-------------------------------- fuel model
    791. 

  791. ! compute derived scalars
    792. 

    793. 

  793. fuelheat = cmbcnst * 4.30e-04     ! convert J/kg to BTU/lb
    794. 

    795. 

  795. ! compute derived fuel category coefficients 
    796. 

    797. 

  797. DO i = 1,nfuelcats
    798. 

  798.     fci(i) = (1.+fuelmc_c)*fci_d(i)
    799. 

  799.     if(fct(i) .ne.  0.)then
    800. 

  800.         fcbr(i) = fci_d(i)/fct(i) ! to avoid division by zero
    801. 

  801.     else
    802. 

  802.         fcbr(i) = 0
    803. 

  803.     endif
    804. 

  804. END DO
    805. 

    806. 

  806. ! prints
    807. 

    808. 

  808. call message('**********************************************************')
    809. 

  809. call message('FUEL COEFFICIENTS')
    810. 

  810. write(msg,8)'cmbcnst    ',cmbcnst
    811. 

  811. call message(msg)
    812. 

  812. write(msg,8)'hfgl       ',hfgl
    813. 

  813. call message(msg)
    814. 

  814. write(msg,8)'fuelmc_g   ',fuelmc_g
    815. 

  815. call message(msg)
    816. 

  816. write(msg,8)'fuelmc_c   ',fuelmc_c
    817. 

  817. call message(msg)
    818. 

  818. write(msg,8)'fuelheat   ',fuelheat
    819. 

  819. call message(msg)
    820. 

  820. write(msg,7)'nfuelcats  ',nfuelcats
    821. 

  821. call message(msg)
    822. 

  822. write(msg,7)'no_fuel_cat',no_fuel_cat
    823. 

  823. call message(msg)
    824. 

  824. write(msg,7)'no_fuel_cat2',no_fuel_cat2
    825. 

  825. call message(msg)
    826. 

  826. if(init_fuel_moisture)then
    827. 

  827.     write(msg,7)'moisture_classes',moisture_classes
    828. 

  828.     call message(msg)
    829. 

  829. endif
    830. 

    831. 

  831. j=1
    832. 

  832. 7 format(a,5(1x,i8,4x))
    833. 

  833. 8 format(a,5(1x,g12.5e2))
    834. 

  834. 9 format(a,5(1x,a))
    835. 

  835. 10 format(a,i2.2,2x,5(1x,g12.5e2))
    836. 

  836. do i=1,nfuelcats,j
    837. 

  837.     k=min(i+j-1,nfuelcats)
    838. 

  838.     call message(' ')
    839. 

  839.     write(msg,7)'CATEGORY  ',(ii,ii=i,k)
    840. 

  840.     call message(msg)
    841. 

  841.     write(msg,9)'fuel name ',(trim(fuel_name(ii)),ii=i,k)
    842. 

  842.     call message(msg)
    843. 

  843.     write(msg,8)'fwh       ',(fcwh(ii),ii=i,k)
    844. 

  844.     call message(msg)
    845. 

  845.     write(msg,8)'fz0       ',(fcz0(ii),ii=i,k)
    846. 

  846.     call message(msg)
    847. 

  847.     write(msg,8)'windrf    ',(windrf(ii),ii=i,k)
    848. 

  848.     call message(msg)
    849. 

  849.     write(msg,8)'fgi       ',(fgi(ii),ii=i,k)
    850. 

  850.     call message(msg)
    851. 

  851.     write(msg,8)'fueldepthm',(fueldepthm(ii),ii=i,k)
    852. 

  852.     call message(msg)
    853. 

  853.     write(msg,8)'savr      ',(savr(ii),ii=i,k)
    854. 

  854.     call message(msg)
    855. 

  855.     write(msg,8)'fuelmce   ',(fuelmce(ii),ii=i,k)
    856. 

  856.     call message(msg)
    857. 

  857.     write(msg,8)'fueldens  ',(fueldens(ii),ii=i,k)
    858. 

  858.     call message(msg)
    859. 

  859.     write(msg,8)'st        ',(st(ii),ii=i,k)
    860. 

  860.     call message(msg)
    861. 

  861.     write(msg,8)'se        ',(se(ii),ii=i,k)
    862. 

  862.     call message(msg)
    863. 

  863.     write(msg,8)'weight    ',(weight(ii),ii=i,k)
    864. 

  864.     call message(msg)
    865. 

  865.     write(msg,8)'fci_d     ',(fci_d(ii),ii=i,k)
    866. 

  866.     call message(msg)
    867. 

  867.     write(msg,8)'fct       ',(fct(ii),ii=i,k)
    868. 

  868.     call message(msg)
    869. 

  869.     write(msg,7)'ichap     ',(ichap(ii),ii=i,k)
    870. 

  870.     call message(msg)
    871. 

  871.     write(msg,8)'fci       ',(fci(ii),ii=i,k)
    872. 

  872.     call message(msg)
    873. 

  873.     write(msg,8)'fcbr      ',(fcbr(ii),ii=i,k)
    874. 

  874.     call message(msg)
    875. 

  875.     write(msg,8)'ffw       ',(ffw(ii),ii=i,k)
    876. 

  876.     call message(msg)
    877. 

  877.     write(msg,8)'adjr0     ',(adjr0(ii),ii=i,k)
    878. 

  878.     call message(msg)
    879. 

  879.     write(msg,8)'adjrw     ',(adjrw(ii),ii=i,k)
    880. 

  880.     call message(msg)
    881. 

  881.     write(msg,8)'adjrs     ',(adjrs(ii),ii=i,k)
    882. 

  882.     call message(msg)
    883. 

  883.     if(init_fuel_moisture)then
    884. 

  884.         do kk=1,moisture_classes 
    885. 

  885.             write(msg,10)'fmc_gw',kk,(fmc_gw(ii,kk),ii=i,k)
    886. 

  886.             call message(msg)
    887. 

  887.         enddo
    888. 

  888.     endif
    889. 

  889.     if(kfmc_ndwi>0)then
    890. 

  890.         write(msg,8)'fmc_gl_stdev     ',(fmc_gl_stdev(ii),ii=i,k)
    891. 

  891.         call message(msg)
    892. 

  892.         write(msg,8)'fmc_gl_ndwi_0    ',(fmc_gl_ndwi_0(ii),ii=i,k)
    893. 

  893.         call message(msg)
    894. 

  894.         write(msg,8)'fmc_gl_ndwi_rate ',(fmc_gl_ndwi_rate(ii),ii=i,k)
    895. 

  895.         call message(msg)
    896. 

  896.         write(msg,8)'fmc_gl_ndwi_stdev',(fmc_gl_ndwi_stdev(ii),ii=i,k)
    897. 

  897.         call message(msg)
    898. 

  898.     endif
    899. 

  899. enddo
    900. 

  900. call message(' ')
    901. 

  901. call message('**********************************************************')
    902. 

    903. 

  903. if(init_fuel_moisture)then
    904. 

  904.     j=1
    905. 

  905.     do i=1,moisture_classes,j
    906. 

  906.         k=min(i+j-1,nfuelcats)
    907. 

  907.         call message(' ')
    908. 

  908.         write(msg,7)'FUEL MOISTURE CLASS',(ii,ii=i,k)
    909. 

  909.         call message(msg)
    910. 

  910.         write(msg,9)'moisture class name    ',(trim(moisture_class_name(ii)),ii=i,k)
    911. 

  911.         call message(msg)
    912. 

  912.         write(msg,7)'drying_model           ',(drying_model(ii),ii=i,k)
    913. 

  913.         call message(msg)
    914. 

  914.         write(msg,8)'drying_lag (h)         ',(drying_lag(ii),ii=i,k)
    915. 

  915.         call message(msg)
    916. 

  916.         write(msg,7)'wetting_model          ',(wetting_model(ii),ii=i,k)
    917. 

  917.         call message(msg)
    918. 

  918.         write(msg,7)'fmc_gc_initialization  ',(fmc_gc_initialization(ii),ii=i,k)
    919. 

  919.         call message(msg)
    920. 

  920.         write(msg,8)'wetting_lag (h)        ',(wetting_lag(ii),ii=i,k)
    921. 

  921.         call message(msg)    
    922. 

  922.         write(msg,8)'saturation_moisture (1)',(saturation_moisture(ii),ii=i,k)
    923. 

  923.         call message(msg)    
    924. 

  924.         write(msg,8)'saturation_rain (mm/h) ',(saturation_rain(ii),ii=i,k)
    925. 

  925.         call message(msg)    
    926. 

  926.         write(msg,8)'rain_threshold (mm/h)  ',(rain_threshold(ii),ii=i,k)
    927. 

  927.         call message(msg)    
    928. 

  928.     enddo
    929. 

  929.     call message(' ')
    930. 

  930.     call message('**********************************************************')
    931. 

  931.     call message(' ')
    932. 

  932. endif
    933. 

  933. have_fuel_cats=.true.
    934. 

    935. 

  935. ! and print to file
    936. 

  936. IF ( wrf_dm_on_monitor() ) THEN
    937. 

  937.   call write_fuels_m(61,30.,1.)
    938. 

  938. ENDIF
    939. 

    940. 

    941. 

  941. end subroutine init_fuel_cats
    942. 

    943. 

    944. 

  944. subroutine write_fuels_m(nsteps,maxwind,maxslope)
    945. 

  945. implicit none
    946. 

  946. integer, intent(in):: nsteps   ! number of steps for speed computation
    947. 

  947. real, intent(in):: maxwind,maxslope ! computer from zero to these
    948. 

    949. 

  949. integer:: iounit,k,j,i,isave
    950. 

  950. type(fire_params)::fp
    951. 

  951. real, dimension(1:3,1:nsteps), target::vx,vy,zsf,dzdxf,dzdyf,bbb,phisc,phiwc,r_0,fgip,ischap,fmc_g,wind,nfuel_cat
    952. 

  952. real, dimension(1:3,1:nsteps)::fuel_time,ros,fwh,fz0
    953. 

  953. real::ros_back,ros_wind,ros_slope,propx,propy,r
    954. 

    955. 

  955. if(.not.have_fuel_cats)call crash('write_fuels_m: fuel categories not yet set')
    956. 

    957. 

  957. fp%vx=>vx
    958. 

  958. fp%vy=>vy
    959. 

  959. fp%dzdxf=>dzdxf
    960. 

  960. fp%dzdyf=>dzdyf
    961. 

  961. fp%bbb=>bbb
    962. 

  962. fp%phisc=>phisc
    963. 

  963. fp%phiwc=>phiwc
    964. 

  964. fp%r_0=>r_0
    965. 

  965. fp%fgip=>fgip
    966. 

  966. fp%ischap=>ischap
    967. 

  967. fp%fmc_g=>fmc_g
    968. 

  968. fp%nfuel_cat=>nfuel_cat
    969. 

    970. 

  970. iounit = open_text_file('fuels.m','write')
    971. 

    972. 

  972. 10 format('fuel(',i3,').',a,'=',"'",a,"'",';% ',a)
    973. 

  973. do k=1,nfuelcats
    974. 

  974.     write(iounit,10)k,'fuel_name',trim(fuel_name(k)),'FUEL MODEL NAME'
    975. 

  975.     call write_var(k,'windrf',windrf(k),'WIND REDUCTION FACTOR FROM FCWH TO MIDFLAME HEIGHT' )
    976. 

  976.     call write_var(k,'fwh',fcwh(k),'WIND HEIGHT TO INTERPOLATE VERTICALLY TO (M)' )
    977. 

  977.     call write_var(k,'fz0',fcz0(k),'ROUGHNESS LENGTH FOR VERTICAL WIND LOG INTERPOLATION (M)' )
    978. 

  978.     call write_var(k,'fgi',fgi(k),'INITIAL TOTAL MASS OF SURFACE FUEL (KG/M**2)' )
    979. 

  979.     call write_var(k,'fueldepthm',fueldepthm(k),'FUEL DEPTH (M)')
    980. 

  980.     call write_var(k,'savr',savr(k),'FUEL PARTICLE SURFACE-AREA-TO-VOLUME RATIO, 1/FT')
    981. 

  981.     call write_var(k,'fuelmce',fuelmce(k),'MOISTURE CONTENT OF EXTINCTION')
    982. 

  982.     call write_var(k,'fueldens',fueldens(k),'OVENDRY PARTICLE DENSITY, LB/FT^3')
    983. 

  983.     call write_var(k,'st',st(k),'FUEL PARTICLE TOTAL MINERAL CONTENT')
    984. 

  984.     call write_var(k,'se',se(k),'FUEL PARTICLE EFFECTIVE MINERAL CONTENT')
    985. 

  985.     call write_var(k,'weight',weight(k),'WEIGHTING PARAMETER THAT DETERMINES THE SLOPE OF THE MASS LOSS CURVE')
    986. 

  986.     call write_var(k,'fci_d',fci_d(k),'INITIAL DRY MASS OF CANOPY FUEL')
    987. 

  987.     call write_var(k,'fct',fct(k),'BURN OUT TIME FOR CANOPY FUEL, AFTER DRY (S)')
    988. 

  988.     call write_var(k,'ichap',float(ichap(k)),'1 if chaparral, 0 if not')
    989. 

  989.     call write_var(k,'fci',fci(k),'INITIAL TOTAL MASS OF CANOPY FUEL')
    990. 

  990.     call write_var(k,'fcbr',fcbr(k),'FUEL CANOPY BURN RATE (KG/M**2/S)')
    991. 

  991.     call write_var(k,'adjr0',adjr0(k),'MULTIPLICATIVE ADJUSTMENT OF BACKING SPREAD RATE')
    992. 

  992.     call write_var(k,'adjrw',adjrw(k),'MULTIPLICATIVE ADJUSTMENT OF WIND CONTRIBUTION TO SPREAD RATE')
    993. 

  993.     call write_var(k,'adjrs',adjrs(k),'MULTIPLICATIVE ADJUSTMENT OF SLOPE CONTRIBUTION TO SPREAD RATE')
    994. 

  994.     call write_var(k,'ffw',ffw(k),'FUEL FRACTION CONSUMED IN THE FLAMING ZONE')
    995. 

  995.     call write_var(k,'hfgl',hfgl,'SURFACE FIRE HEAT FLUX THRESHOLD TO IGNITE CANOPY (W/m^2)')
    996. 

  996.     call write_var(k,'cmbcnst',cmbcnst,'JOULES PER KG OF DRY FUEL')
    997. 

  997.     call write_var(k,'fuelheat',fuelheat,'FUEL PARTICLE LOW HEAT CONTENT, BTU/LB')
    998. 

  998.     call write_var(k,'fuelmc_g',fuelmc_g,'FUEL PARTICLE (SURFACE) MOISTURE CONTENT')
    999. 

  999.     call write_var(k,'fuelmc_c',fuelmc_c,'FUEL PARTICLE (CANOPY) MOISTURE CONTENT')
    1000. 

  1000.     ! set up fuel arrays 
    1001. 

  1001.     !subroutine set_fire_params( &
    1002. 

  1002.     !                       ifds,ifde,jfds,jfde, &
    1003. 

  1003.     !                       ifms,ifme,jfms,jfme, &
    1004. 

  1004.     !                       ifts,ifte,jfts,jfte, &
    1005. 

  1005.     !                       fdx,fdy,nfuel_cat0,  &
    1006. 

  1006.     !                       nfuel_cat,fuel_time, &
    1007. 

  1007.     !                       fp ) 
    1008. 

  1008.     nfuel_cat = k
    1009. 

  1009.     do j=1,nsteps ! set moisture - must be before set_fire_params
    1010. 

  1010.        fmc_g(1,j)=fuelmc_g
    1011. 

  1011.        fmc_g(2,j)=fuelmc_g
    1012. 

  1012.        fmc_g(3,j)=(fuelmce(k)*(j-1))/(nsteps-2)
    1013. 

  1013.     enddo
    1014. 

  1014.     isave=fire_fmc_read
    1015. 

  1015.     fire_fmc_read=0
    1016. 

  1016.     call set_fire_params( &
    1017. 

  1017.                            1,3,1,nsteps, &
    1018. 

  1018.                            1,3,1,nsteps, &
    1019. 

  1019.                            1,3,1,nsteps, &
    1020. 

  1020.                            0.,0.,k,  &
    1021. 

  1021.                            nfuel_cat,fuel_time, &
    1022. 

  1022.                            fp ) 
    1023. 

  1023.     fire_fmc_read=isave
    1024. 

  1024.     ! set up windspeed slope moisture table
    1025. 

  1025.     propx=1.
    1026. 

  1026.     propy=0.
    1027. 

  1027.     do j=1,nsteps
    1028. 

  1028.        r=float(j-1)/float(nsteps-1)
    1029. 

  1029.        ! line 1 varies windspeed (in x direction), zero slope
    1030. 

  1030.        wind(1,j)=maxwind*r
    1031. 

  1031.        vx(1,j)=wind(1,j)*windrf(k)
    1032. 

  1032.        vy(1,j)=0.
    1033. 

  1033.        dzdxf(1,j)=0.
    1034. 

  1034.        dzdyf(1,j)=0.
    1035. 

  1035.        ! line 2 varies slope (in x direction), zero slope
    1036. 

  1036.        vx(2,j)=0.
    1037. 

  1037.        vy(2,j)=0.
    1038. 

  1038.        dzdxf(2,j)=maxslope*r
    1039. 

  1039.        dzdyf(2,j)=0.
    1040. 

  1040.        ! line 3 varies moisture, zero slope and wind
    1041. 

  1041.        vx(3,j)=0.
    1042. 

  1042.        vy(3,j)=0.
    1043. 

  1043.        dzdxf(3,j)=0.
    1044. 

  1044.        dzdyf(3,j)=0.
    1045. 

  1045.     enddo
    1046. 

  1046.     do j=1,nsteps
    1047. 

  1047.        do i=1,3
    1048. 

  1048.           call fire_ros(ros_back,ros_wind,ros_slope, &
    1049. 

  1049.              propx,propy,i,j,fp)
    1050. 

  1050.           ros(i,j)=ros_back+ros_wind+ros_slope
    1051. 

  1051.        enddo
    1052. 

  1052.        write(iounit,13)k,'wind',j,wind(1,j),'wind speed at 6.1m'
    1053. 

  1053.        write(iounit,13)k,'ros_wind',j,ros(1,j),'rate of spread for the wind speed at 6.1m'
    1054. 

  1054.        write(iounit,13)k,'slope',j,dzdxf(2,j),'slope'
    1055. 

  1055.        write(iounit,13)k,'ros_slope',j,ros(2,j),'rate of spread for the slope'
    1056. 

  1056.        write(iounit,13)k,'fmc_g',j,fmc_g(3,j),'fuel moisture content'
    1057. 

  1057.        write(iounit,13)k,'ros_fmc_g',j,ros(3,j),'rate of spread for the fuel moisture content'
    1058. 

  1058.     enddo
    1059. 

  1059. enddo
    1060. 

  1060. 13 format('fuel(',i3,').',a,'(',i3,')=',g12.5e2,';% ',a)
    1061. 

  1061.  
    1062. 

  1062. close(iounit)
    1063. 

  1063. ! stop
    1064. 

    1065. 

  1065. contains
    1066. 

    1067. 

  1067. subroutine write_var(k,name,value,descr)
    1068. 

  1068. ! write entry for one variable
    1069. 

  1069. integer, intent(in)::k
    1070. 

  1070. character(len=*), intent(in)::name,descr
    1071. 

  1071. real, intent(in)::value
    1072. 

  1072. write(iounit,11)k,name,value
    1073. 

  1073. write(iounit,12)k,name,descr
    1074. 

  1074. 11 format('fuel(',i3,').',a,'=',g12.5e2,  ';')
    1075. 

  1075. 12 format('fuel(',i3,').',a,"_descr='",a,"';")
    1076. 

  1076. end subroutine write_var
    1077. 

    1078. 

  1078. end subroutine write_fuels_m
    1079. 

    1080. 

  1080. !
    1081. 

  1081. !*******************
    1082. 

  1082. !
    1083. 

    1084. 

  1084. subroutine set_fire_params( &
    1085. 

  1085.                            ifds,ifde,jfds,jfde, &
    1086. 

  1086.                            ifms,ifme,jfms,jfme, &
    1087. 

  1087.                            ifts,ifte,jfts,jfte, &
    1088. 

  1088.                            fdx,fdy,nfuel_cat0,  &
    1089. 

  1089.                            nfuel_cat,fuel_time, &
    1090. 

  1090.                            fp ) 
    1091. 

    1092. 

  1092. implicit none
    1093. 

    1094. 

  1094. !*** purpose: Set all fire model params arrays, constant values.
    1095. 

    1096. 

  1096. !*** arguments
    1097. 

  1097. integer, intent(in)::ifds,ifde,jfds,jfde                        ! fire domain bounds
    1098. 

  1098. integer, intent(in)::ifts,ifte,jfts,jfte                        ! fire tile bounds
    1099. 

  1099. integer, intent(in)::ifms,ifme,jfms,jfme                        ! memory array bounds
    1100. 

  1100. real, intent(in):: fdx,fdy                                      ! fire mesh spacing
    1101. 

  1101. integer,intent(in)::nfuel_cat0                                  ! default fuel category, if nfuel_cat=0
    1102. 

  1102. real, intent(inout),dimension(ifms:ifme, jfms:jfme)::nfuel_cat  ! fuel data
    1103. 

  1103. real, intent(out), dimension(ifms:ifme, jfms:jfme)::fuel_time   ! fire params arrays
    1104. 

  1104. type(fire_params),intent(inout)::fp
    1105. 

    1106. 

  1106. !*** local
    1107. 

    1108. 

  1108. real::  fuelload, fueldepth, rtemp1, rtemp2, &
    1109. 

  1109.         qig, epsilon, rhob, wn, betaop, e, c, &
    1110. 

  1110.         xifr, etas, etam, a, gammax, gamma, ratio, ir, &
    1111. 

  1111.         fuelloadm,fdxinv,fdyinv,betafl, bmst
    1112. 

  1112. integer:: i,j,k
    1113. 

  1113. integer::nerr
    1114. 

  1114. character(len=128)::msg
    1115. 

    1116. 

  1116. !*** executable
    1117. 

    1118. 

  1118. if(.not.have_fuel_cats)call crash('set_fire_params: fuel categories not yet set')
    1119. 

    1120. 

  1120. call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%fmc_g,'set_fire_params: fmc_g')
    1121. 

    1122. 

  1122. nerr=0
    1123. 

  1123. do j=jfts,jfte
    1124. 

  1124.    do i=ifts,ifte
    1125. 

  1125.      ! fuel category 
    1126. 

  1126.      k=int( nfuel_cat(i,j) )
    1127. 

  1127.      if(k.ge.no_fuel_cat.and.k.le.no_fuel_cat2)then   ! no fuel 
    1128. 

  1128.         fp%fgip(i,j)=0.            ! no mass 
    1129. 

  1129.         fp%ischap(i,j)=0.
    1130. 

  1130.         fp%phisc(i,j)=0.          ! 
    1131. 

  1131.         fp%bbb(i,j)=1.             !
    1132. 

  1132.         fuel_time(i,j)=7./0.85  ! does not matter, just what was there before
    1133. 

  1133.         fp%phiwc(i,j)=0.
    1134. 

  1134.         fp%r_0(i,j)=0.             ! no fuel, no spread.
    1135. 

  1135.      else
    1136. 

  1136.         ! if(k.eq.0.and.nfuel_cat0.ge.1.and.nfuel_cat0.le.nfuelcats)then
    1137. 

  1137.         !     ! replace k=0 by default
    1138. 

  1138.         !     k=nfuel_cat0
    1139. 

  1139.         !     nerr=nerr+1
    1140. 

  1140.         ! endif
    1141. 

  1141.    
    1142. 

  1142.         if(k.lt.1.or.k.gt.nfuelcats)then
    1143. 

  1143. !$OMP CRITICAL(SFIRE_PHYS_CRIT)
    1144. 

  1144.             write(msg,'(3(a,i5))')'nfuel_cat(', i ,',',j,')=',k
    1145. 

  1145. !$OMP END CRITICAL(SFIRE_PHYS_CRIT)
    1146. 

  1146.             call message(msg)
    1147. 

  1147.             if(k.eq.0)then
    1148. 

  1148.                 call message('Possibly nfuel_cat is uninitialized on input')
    1149. 

  1149.             endif
    1150. 

  1150.             call crash('set_fire_params: fuel category out of bounds')
    1151. 

  1151.         endif
    1152. 

    1153. 

  1153.         fuel_time(i,j)=weight(k)/0.85 ! cell based
    1154. 

  1154.         
    1155. 

  1155.         ! do not understand calculations of stime in binit.m4
    1156. 

  1156.         ! set fuel time constant: weight=1000=>40% decrease over 10 min
    1157. 

  1157.         ! fuel decreases as exp(-t/fuel_time) 
    1158. 

  1158.         ! exp(-600*0.85/1000) = approx 0.6 
    1159. 

    1160. 

  1160.         fp%ischap(i,j)=ichap(k)
    1161. 

  1161.         fp%fgip(i,j)=fgi(k)
    1162. 

  1162.         if(fire_fmc_read.eq.1)then
    1163. 

  1163.            fp%fmc_g(i,j)=fuelmc_g
    1164. 

  1164.         endif
    1165. 

  1165.         ! print *,'fmc_g:',fire_fmc_read,i,j,fp%fmc_g(i,j)
    1166. 

    1167. 

  1167.         ! end jm addition
    1168. 

    1169. 

  1169.         !
    1170. 

  1170.         !*** rest copied from wf2_janice/fire_startup.m4 with minimal changes
    1171. 

  1171.         !
    1172. 

    1173. 

  1173.         !     ...Settings of fire spread parameters from Rothermel follows. These
    1174. 

  1174.         !        don't need to be recalculated later.
    1175. 

  1175.         
    1176. 

  1176.         bmst     = fp%fmc_g(i,j) / (1.+fp%fmc_g(i,j))
    1177. 

  1177.         fuelloadm= (1.-bmst) * fgi(k)  !  fuelload without moisture
    1178. 

  1178.         fuelload = fuelloadm * (.3048)**2 * 2.205    ! to lb/ft^2
    1179. 

  1179.         fueldepth = fueldepthm(k)/0.3048               ! to ft
    1180. 

  1180.         betafl = fuelload/(fueldepth * fueldens(k))! packing ratio
    1181. 

  1181.         betaop = 3.348 * savr(k)**(-0.8189)     ! optimum packing ratio
    1182. 

  1182.         qig = 250. + 1116.*fp%fmc_g(i,j)            ! heat of preignition, btu/lb
    1183. 

  1183.         epsilon = exp(-138./savr(k) )    ! effective heating number
    1184. 

  1184.         rhob = fuelload/fueldepth    ! ovendry bulk density, lb/ft^3
    1185. 

    1186. 

  1186.         c = 7.47 * exp( -0.133 * savr(k)**0.55)    ! const in wind coef
    1187. 

  1187.         fp%bbb(i,j) = 0.02526 * savr(k)**0.54      ! const in wind coef
    1188. 

  1188.         !if(fire_wind_log_interp .eq. 4 .or. fire_use_windrf .eq. 1) then
    1189. 

  1189.         !     c = c * windrf(k)**fp%bbb(i,j)             ! jm: multiply wind by reduction factor
    1190. 

  1190.         !endif
    1191. 

  1191.         e = 0.715 * exp( -3.59e-4 * savr(k))       ! const in wind coef
    1192. 

  1192.         fp%phiwc(i,j) = c * (betafl/betaop)**(-e)
    1193. 

    1194. 

  1194.         ! phis = 5.275 *(fp%betafl(i,j))**(-0.3) *tanphim**2   ! slope factor
    1195. 

  1195.         fp%phisc(i,j) = 5.275 *(betafl)**(-0.3)    ! const in slope coeff
    1196. 

    1197. 

  1197.         rtemp2 = savr(k)**1.5
    1198. 

  1198.         gammax = rtemp2/(495. + 0.0594*rtemp2)              ! maximum rxn vel, 1/min
    1199. 

  1199.         a = 1./(4.774 * savr(k)**0.1 - 7.27)   ! coef for optimum…
    1200.

Large files files are truncated, but you can click here to view the full file