forked from wrf-model/WRF
/
module_fr_sfire_atm.F
180 lines (130 loc) · 5.31 KB
/
module_fr_sfire_atm.F
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
!WRF:MEDIATION_LAYER:FIRE_MODEL
!*** Jan Mandel August 2007 - February 2008
!*** email: jmandel@ucar.edu or Jan.Mandel@gmail.com or Jan.Mandel@cudenver.edu
! Routines dealing with the atmosphere
module module_fr_sfire_atm
use module_model_constants, only: cp,xlv
use module_fr_sfire_util
contains
SUBROUTINE fire_tendency( &
ids,ide, kds,kde, jds,jde, & ! dimensions
ims,ime, kms,kme, jms,jme, &
its,ite, kts,kte, jts,jte, &
grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to atm grid
alfg,alfc,z1can, & ! coeffients, properties, geometry
zs,z_at_w,dz8w,mu,rho, &
rthfrten,rqvfrten) ! theta and Qv tendencies
! This routine is atmospheric physics
! it does NOT go into module_fr_sfire_phys because it is not fire physics
! taken from the code by Ned Patton, only order of arguments change to the convention here
! --- this routine takes fire generated heat and moisture fluxes and
! calculates their influence on the theta and water vapor
! --- note that these tendencies are valid at the Arakawa-A location
IMPLICIT NONE
! --- incoming variables
INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
ims,ime, kms,kme, jms,jme, &
its,ite, kts,kte, jts,jte
REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: grnhfx,grnqfx ! W/m^2
REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: canhfx,canqfx ! W/m^2
REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: zs ! topography (m abv sealvl)
REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: mu ! dry air mass (Pa)
REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: z_at_w ! m abv sealvl
REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: dz8w ! dz across w-lvl
REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: rho ! density
REAL, INTENT(in) :: alfg ! extinction depth ground fire heat (m)
REAL, INTENT(in) :: alfc ! extinction depth crown fire heat (m)
REAL, INTENT(in) :: z1can ! height of crown fire heat release (m)
! --- outgoing variables
REAL, INTENT(out), DIMENSION( ims:ime,kms:kme,jms:jme ) :: &
rthfrten, & ! theta tendency from fire (in mass units)
rqvfrten ! Qv tendency from fire (in mass units)
! --- local variables
INTEGER :: i,j,k
INTEGER :: i_st,i_en, j_st,j_en, k_st,k_en
REAL :: cp_i
REAL :: rho_i
REAL :: xlv_i
REAL :: z_w
REAL :: fact_g, fact_c
REAL :: alfg_i, alfc_i
REAL, DIMENSION( its:ite,kts:kte,jts:jte ) :: hfx,qfx
!! character(len=128)::msg
do j=jts,jte
do k=kts,min(kte+1,kde)
do i=its,ite
rthfrten(i,k,j)=0.
rqvfrten(i,k,j)=0.
enddo
enddo
enddo
! --- set some local constants
cp_i = 1./cp ! inverse of specific heat
xlv_i = 1./xlv ! inverse of latent heat
alfg_i = 1./alfg
alfc_i = 1./alfc
!!write(msg,'(8e11.3)')cp,cp_i,xlv,xlv_i,alfg,alfc,z1can
!!call message(msg)
call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_tendency:grnhfx')
call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_tendency:grnqfx')
! --- set loop indicies : note that
i_st = MAX(its,ids+1)
i_en = MIN(ite,ide-1)
k_st = kts
k_en = MIN(kte,kde-1)
j_st = MAX(jts,jds+1)
j_en = MIN(jte,jde-1)
! --- distribute fluxes
DO j = j_st,j_en
DO k = k_st,k_en
DO i = i_st,i_en
! --- set z (in meters above ground)
z_w = z_at_w(i,k,j) - zs(i,j) ! should be zero when k=k_st
! --- heat flux
fact_g = cp_i * EXP( - alfg_i * z_w )
IF ( z_w < z1can ) THEN
fact_c = cp_i
ELSE
fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) )
END IF
hfx(i,k,j) = fact_g * grnhfx(i,j) + fact_c * canhfx(i,j)
!! write(msg,2)i,k,j,z_w,grnhfx(i,j),hfx(i,k,j)
!!2 format('hfx:',3i4,6e11.3)
!! call message(msg)
! --- vapor flux
fact_g = xlv_i * EXP( - alfg_i * z_w )
IF (z_w < z1can) THEN
fact_c = xlv_i
ELSE
fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) )
END IF
qfx(i,k,j) = fact_g * grnqfx(i,j) + fact_c * canqfx(i,j)
!! if(hfx(i,k,j).ne.0. .or. qfx(i,k,j) .ne. 0.)then
!! write(msg,1)i,k,j,hfx(i,k,j),qfx(i,k,j)
!!1 format('tend:',3i6,2e11.3)
!! call message(msg)
! endif
END DO
END DO
END DO
! --- add flux divergence to tendencies
!
! multiply by dry air mass (mu) to eliminate the need to
! call sr. calculate_phy_tend (in dyn_em/module_em.F)
DO j = j_st,j_en
DO k = k_st,k_en-1
DO i = i_st,i_en
rho_i = 1./rho(i,k,j)
rthfrten(i,k,j) = - mu(i,j) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j)
rqvfrten(i,k,j) = - mu(i,j) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j)
END DO
END DO
END DO
call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rthfrten,'fire_tendency:rthfrten')
call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rqvfrten,'fire_tendency:rqvfrten')
RETURN
END SUBROUTINE fire_tendency
!
!***
!
end module module_fr_sfire_atm