-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbalan.F90
273 lines (222 loc) · 9.64 KB
/
balan.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
!...end-of-timestep calculation to determine if water and energy
!...balance is maintained
subroutine balan(sib,nsl_old,wwwliq_old,wwwice_old,capac_old,cas_q)
use kinds
use sibtype
use sib_const_module, only: &
nsoil, &
snofac, &
dtt, &
dti, &
tau
use physical_parameters, only: &
hltm
implicit none
!-----------------------------------------------------------------------
!...WATER BALANCE...
!
! WATER IN + CHANGE IN STORAGE = WATER OUT
!
! inputs:
! precipitation (sib_prog%cupr,sib_prog%lspr) (mm/sec)
!
! output2:
! CAS-mixed layer latent heat flux (sib_diag%fws) (W/m^2)
!
! runoff
! surface runoff (sib_diag%roffo) (mm/hr)
! subsurface runoff (sib_diag%roff) (mm/hr)
!
! storage
! interception reservoirs
! canopy interception (sib_prog%capac(1)) (kg/m^2)
! ground interception (sib_prog%capac(2)) (kg/m^2)
! soil
! soil water (sib_prog%www_liq) (kg/m^2)
! soil ice (sib_prog%www_ice) (kg/m^2)
! snow
! snow water (sib_prog%www_liq) (kg/m^2)
! snow ice (sib_prog%www_ice) (kg/m^2)
!
! canopy airspace (cas_q) (kg/m^2)
!
! ** note: indices for soil water and ice arrays
! will be 1 to 10 for soil (1 at top,
! 10 at bottom of soil column), and from
! -4 (top) to 0 (bottom) for up to 5
! snow layers.
!
!----------------------------------------------------------------------
!...ENERGY BALANCE
!
! ENERGY IN + CHANGE IN STORAGE = ENERGY OUT
!
! inputs:
! radiation absorbed by canopy (sib_diag%radt(1)) (W/m^2)
! radiation absorbed by ground (sib_diag%radt(2)) (W/m^2)
!
! outputs:
! latent heat
! transpiration (sib_diag%ect) (J/m^2)
! evaporation
! from canopy interception storage (sib_diag%eci) (J/m^2)
! from ground interception storage (sib_diag%egi) (J/m^2)
! from surface soil layer (sib_diag%egs) (J/m^2)
! from snow (sib_diag%ess) (J/m^2)
! sensible heat
! canopy (sib_diag%hc) (J/m^2)
! ground (sib_diag%hg) (J/m^2)
! snow (sib_diag%hs) (J/m^2)
!
! storage
! canopy storage flux (sib_diag%chf) (W/m^2)
! soil storage flux (sib_diag%shf) (W/m^2)
!
!
!----------------------------------------------------------------------
type(sib_t), intent(inout) :: sib
!-----------------------------------------------------------------------
! input variables
!-----------------------------------------------------------------------
integer(kind=int_kind),intent(in) :: nsl_old
real(kind=dbl_kind),intent(in), &
dimension(nsl_old+1:nsoil) :: wwwliq_old ! (kg/m^2)
real(kind=dbl_kind),intent(in), &
dimension(nsl_old+1:nsoil) :: wwwice_old ! (kg/m^2)
real(kind=dbl_kind),intent(in) :: capac_old(2) ! (kg/m^2)
real(kind=dbl_kind),intent(in) :: cas_q ! (kg/m^2)
!-----------------------------------------------------------------------
! end input variables
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
real(kind=dbl_kind) :: dstor ! change in storage, canopy and
! surface interception (kg/m^2)
real(kind=dbl_kind) :: dqsoil ! change in storage, soil only
! liquid and ice (kg/m^2)
real(kind=dbl_kind) :: dqsnow ! change in storage, snow (kg/m^2)
real(kind=dbl_kind) :: snownew ! amount of snow, end-of-timestep (kg/m^2)
real(kind=dbl_kind) :: snowold ! amount of snow, beg-of-timestep (kg/m^2)
real(kind=dbl_kind) :: evap ! evaporation, from CAS to atmosphere.
! this encompasses CAS storage,
! ground interception stores, and
! evap from top soil layer as well
! as snow (kg/m^2)
real(kind=dbl_kind) :: runoff ! total runoff: surface + subsurface
! (kg/m^2)
real(kind=dbl_kind) :: sbeg,send! sum holders
real(kind=dbl_kind) :: precip ! precip (kg/m^2)
real(kind=dbl_kind) :: cas_q_new ! CAS water stored as vapor (kg/m^2)
real(kind=dbl_kind) :: vbal ! vapor balance (kg/m^2)
real(kind=dbl_kind) :: wbal ! water balance
real(kind=dbl_kind) :: rhs,lhs ! right and left hand sides of energy
! balance equation (W/m^2)
real(kind=dbl_kind) :: ebal ! energy imbalance (W/m^2)
integer(kind=int_kind) :: i ! loop index
!-----------------------------------------------------------------------
! end local variables
!-----------------------------------------------------------------------
wbal = 0.0
ebal = 0.0
!...NOTE...unless otherwise noted, all units will be kg/m^2 water
!...change in canopy and surface interception storage
dstor = (sib%prog%capac(1) - capac_old(1)) + & ! canopy interception
(sib%prog%capac(2) - capac_old(2)) ! surface interception
!...change in soil water
dqsoil = 0.0
do i=1,nsoil
dqsoil = dqsoil + (sib%prog%www_liq(i) - wwwliq_old(i)) + &
(sib%prog%www_ice(i) - wwwice_old(i))
enddo
!...beginning of timestep snow water (liquid + ice)
snowold = 0.0
do i=nsl_old+1,0
snowold = snowold + wwwliq_old(i) + wwwice_old(i)
enddo
!...end of timestep snow water (liquid + ice)
snownew = 0.0
do i=sib%prog%nsl+1,0
snownew = snownew + sib%prog%www_liq(i) + sib%prog%www_ice(i)
enddo
dqsnow = snownew - snowold
!...evaporation, from canopy and ground surface interception
!...stores (sib%diag%eci, sib%diag%egi),
!...as well as top soil layer (sib%diag%egs) and snow evaporation
evap = sib%diag%fws * dtt / hltm ! total CAS water vapor flux
!...runoff, surface + subsurface
runoff = (sib%diag%roffo + sib%diag%roff)*dtt/3600.0_dbl_kind
!...precip
precip = (sib%prog%lspr + sib%prog%cupr) * dtt
!...vapor
vbal = cas_q_new - cas_q
!...WATER BALANCE...
wbal = precip - (evap + runoff) - &
(dqsoil + dqsnow + dstor)
!...setting water balance error at 0.1 kg/m^2 (or 0.1mm)
if(abs(wbal) > 1.0_dbl_kind) then
! print*,'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
! print'(2(a,g16.6))','WATER IMBALANCE: hour=',tau,'imbalance=' &
! ,wbal,' kg/m^2'
! print'(a)','TERMS: input - output - storage (all units kg/m^2)'
! print'(3(a,g14.6))','input=',precip,' output=',evap+runoff &
! ,' storage=',dqsoil + dqsnow + dstor
! print'(2(a,g16.6))','evap(fws)=',evap,' runoff=',runoff
! print*,' '
! print'(a,4g14.6)','output (eci,egi,egs,ess):',sib%diag%eci*dtt/hltm &
! ,sib%diag%egi*dtt/hltm,sib%diag%egs*dtt/hltm &
! ,sib%diag%ess*dtt*snofac/hltm
! print'(a,4f14.6)','storage (interception, soil, snow,vapor):' &
! ,dstor,dqsoil,dqsnow,vbal
! print'(a,2f14.6)','runoff (overland, subsfc):' &
! ,sib%diag%roffo*dtt/3600.0_dbl_kind &
! ,sib%diag%qqq*dtt/3600.0_dbl_kind
! print*,'soil moisture'
sbeg = 0.0
send = 0.0
do i=1,nsoil
! print'(i4,2(a,g14.6))',i,' beg liq=',wwwliq_old(i),' end liq=' &
! ,sib%prog%www_liq(i)
sbeg = sbeg + wwwliq_old(i)
send = send + sib%prog%www_liq(i)
enddo
! print'(2(a,g14.6))','sum beg=',sbeg,' sum end=',send
sbeg = 0.0
send = 0.0
! print*,'soil ice'
do i=1,nsoil
! print'(i4,2(a,g14.6))',i,' beg ice=',wwwice_old(i),' end ice=' &
! ,sib%prog%www_ice(i)
sbeg = sbeg + wwwice_old(i)
send = send + sib%prog%www_ice(i)
enddo
! print'(2(a,g14.6))','sum beg=',sbeg,' sum end=',send
! print*,'WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW'
endif
!----------------- END OF WATER BALANCE CODE ----------------------------
!----------------- ENERGY BALANCE ---------------------------------------
rhs = (sib%diag%ect + sib%diag%eci + sib%diag%egi + sib%diag%egs + &
sib%diag%ess + sib%diag%hc + sib%diag%hg + sib%diag%hs) * dti
lhs = sib%diag%radt(1) + sib%diag%radt(2) - sib%diag%chf - sib%diag%shf
ebal = lhs - rhs
if (ebal > 1.0E10) then
! if (ebal /= 0.0) then
print*,'EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE'
print*, 'point',sib%stat%pt_num
print'(4(a,g14.6))','ENERGY IMBALANCE: hour =',tau,'rhs=' &
,rhs,' lhs=',lhs,' imbalance=',ebal
print*,'imbalance amount=',ebal,' (W/m^2)'
print'(2(a,g14.6))','sib%diag%ect=',sib%diag%ect*dti,' sib%diag%eci=' &
,sib%diag%eci*dti
print'(2(a,g14.6))','sib%diag%egi=',sib%diag%egi*dti,' sib%diag%egs=' &
,sib%diag%egs*dti
print'(2(a,g14.6))','sib%diag%ess=',sib%diag%ess*dti,' sib%diag%hc=' &
,sib%diag%hc*dti
print'(2(a,g14.6))','sib%diag%hg=',sib%diag%hg*dti,' sib%diag%hs=' &
,sib%diag%hs*dti
print'(2(a,g14.6))','radt(1)=',sib%diag%radt(1),' radt(2)=', &
sib%diag%radt(2)
print'(2(a,g14.6))','sib%diag%chf=',sib%diag%chf,' sib%diag%shf=', &
sib%diag%shf
print*,'EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE'
endif
end subroutine balan