-
Notifications
You must be signed in to change notification settings - Fork 102
/
check_icond.f90
executable file
·288 lines (252 loc) · 19.3 KB
/
check_icond.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
! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
module check_icond_module
USE nrtype
! access missing values
USE globalData,only:integerMissing ! missing integer
USE globalData,only:realMissing ! missing double precision number
! define modeling decisions
USE mDecisions_module,only: &
moisture, & ! moisture-based form of Richards' equation
mixdform ! mixed form of Richards' equation
implicit none
private
public::check_icond
contains
! ************************************************************************************************
! public subroutine check_icond: read model initial conditions
! ************************************************************************************************
subroutine check_icond(nGRU, & ! number of GRUs and HRUs
progData, & ! model prognostic (state) variables
mparData, & ! model parameters
indxData, & ! layer index data
err,message) ! error control
! --------------------------------------------------------------------------------------------------------
! modules
USE nrtype
USE var_lookup,only:iLookParam ! variable lookup structure
USE var_lookup,only:iLookProg ! variable lookup structure
USE var_lookup,only:iLookIndex ! variable lookup structure
USE globalData,only:gru_struc ! gru-hru mapping structures
USE data_types,only:gru_hru_doubleVec ! actual data
USE data_types,only:gru_hru_intVec ! actual data
USE globaldata,only:iname_soil,iname_snow ! named variables to describe the type of layer
USE multiconst,only:&
LH_fus, & ! latent heat of fusion (J kg-1)
iden_ice, & ! intrinsic density of ice (kg m-3)
iden_water,& ! intrinsic density of liquid water (kg m-3)
gravity, & ! gravitational acceleration (m s-2)
Tfreeze ! freezing point of pure water (K)
USE snow_utils_module,only:fracliquid ! compute volumetric fraction of liquid water in snow based on temperature
USE updatState_module,only:updateSnow ! update snow states
USE updatState_module,only:updateSoil ! update soil states
implicit none
! --------------------------------------------------------------------------------------------------------
! variable declarations
! dummies
integer(i4b) ,intent(in) :: nGRU ! number of grouped response units
type(gru_hru_doubleVec),intent(inout) :: progData ! prognostic vars
type(gru_hru_doubleVec),intent(in) :: mparData ! parameters
type(gru_hru_intVec) ,intent(in) :: indxData ! layer indexes
integer(i4b) ,intent(out) :: err ! error code
character(*) ,intent(out) :: message ! returned error message
! locals
character(len=256) :: cmessage ! downstream error message
integer(i4b) :: iGRU ! loop index
integer(i4b) :: iHRU ! loop index
! temporary variables for realism checks
integer(i4b) :: iLayer ! index of model layer
integer(i4b) :: iSoil ! index of soil layer
real(rkind) :: fLiq ! fraction of liquid water on the vegetation canopy (-)
real(rkind) :: vGn_m ! van Genutchen "m" parameter (-)
real(rkind) :: tWat ! total water on the vegetation canopy (kg m-2)
real(rkind) :: scalarTheta ! liquid water equivalent of total water [liquid water + ice] (-)
real(rkind) :: h1,h2 ! used to check depth and height are consistent
integer(i4b) :: nLayers ! total number of layers
real(rkind) :: kappa ! constant in the freezing curve function (m K-1)
integer(i4b) :: nSnow ! number of snow layers
real(rkind),parameter :: xTol=1.e-10_rkind ! small tolerance to address precision issues
real(rkind),parameter :: canIceTol=1.e-3_rkind ! small tolerance to allow existence of canopy ice for above-freezing temperatures (kg m-2)
! --------------------------------------------------------------------------------------------------------
! Start procedure here
err=0; message="check_icond/"
! --------------------------------------------------------------------------------------------------------
! Check that the initial conditions do not conflict with parameters, structure, etc.
! --------------------------------------------------------------------------------------------------------
do iGRU = 1,nGRU
do iHRU = 1,gru_struc(iGRU)%hruCount
! ensure the spectral average albedo is realistic
if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) > mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMax)%dat(1)) &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMax)%dat(1)
if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) < mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinWinter)%dat(1)) &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinWinter)%dat(1)
! ensure the visible albedo is realistic
if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) > mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxVisible)%dat(1)) &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxVisible)%dat(1)
if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) < mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinVisible)%dat(1)) &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinVisible)%dat(1)
! ensure the nearIR albedo is realistic
if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) > mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxNearIR)%dat(1)) &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMaxNearIR)%dat(1)
if(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) < mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinNearIR)%dat(1)) &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) = mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%albedoMinNearIR)%dat(1)
end do
end do
! ensure the initial conditions are consistent with the constitutive functions
do iGRU = 1,nGRU
do iHRU = 1,gru_struc(iGRU)%hruCount
! associate local variables with variables in the data structures
associate(&
! state variables in the vegetation canopy
scalarCanopyTemp => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarCanopyTemp)%dat(1) , & ! canopy temperature
scalarCanopyIce => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarCanopyIce)%dat(1) , & ! mass of ice on the vegetation canopy (kg m-2)
scalarCanopyLiq => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarCanopyLiq)%dat(1) , & ! mass of liquid water on the vegetation canopy (kg m-2)
! state variables in the snow+soil domain
mLayerTemp => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerTemp)%dat , & ! temperature (K)
mLayerVolFracLiq => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracLiq)%dat , & ! volumetric fraction of liquid water in each snow layer (-)
mLayerVolFracIce => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracIce)%dat , & ! volumetric fraction of ice in each snow layer (-)
mLayerMatricHead => progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerMatricHead)%dat , & ! matric head (m)
mLayerLayerType => indxData%gru(iGRU)%hru(iHRU)%var(iLookINDEX%layerType)%dat , & ! type of layer (ix_soil or ix_snow)
! depth varying soil properties
vGn_alpha => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%vGn_alpha)%dat , & ! van Genutchen "alpha" parameter (m-1)
vGn_n => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%vGn_n)%dat , & ! van Genutchen "n" parameter (-)
theta_sat => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%theta_sat)%dat , & ! soil porosity (-)
theta_res => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%theta_res)%dat , & ! soil residual volumetric water content (-)
! snow parameters
snowfrz_scale => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%snowfrz_scale)%dat(1) , & ! scaling parameter for the snow freezing curve (K-1)
FCapil => mparData%gru(iGRU)%hru(iHRU)%var(iLookPARAM%FCapil)%dat(1) & ! fraction of pore space in tension storage (-)
) ! (associate local variables with model parameters)
! compute the constant in the freezing curve function (m K-1)
kappa = (iden_ice/iden_water)*(LH_fus/(gravity*Tfreeze)) ! NOTE: J = kg m2 s-2
! check canopy ice content for unrealistic situations
if(scalarCanopyIce > canIceTol .and. scalarCanopyTemp > Tfreeze)then
! ice content > threshold, terminate run
write(message,'(A,E22.16,A,E9.3,A,F7.3,A,F7.3,A)') trim(message)//'canopy ice (=',scalarCanopyIce,') > canIceTol (=',canIceTol,') when canopy temperature (=',scalarCanopyTemp,') > Tfreeze (=',Tfreeze,')'
err=20; return
else if(scalarCanopyIce > 0._rkind .and. scalarCanopyTemp > Tfreeze)then
! if here, ice content < threshold. Could be sublimation on previous timestep or simply wrong input. Print a warning
write(*,'(A,E22.16,A,F7.3,A,F7.3,A)') 'Warning: canopy ice content in restart file (=',scalarCanopyIce,') > 0 when canopy temperature (=',scalarCanopyTemp,') > Tfreeze (=',Tfreeze,'). Continuing.',NEW_LINE('a')
end if
! number of layers
nLayers = gru_struc(iGRU)%hruInfo(iHRU)%nSnow + gru_struc(iGRU)%hruInfo(iHRU)%nSoil
nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow
! loop through all layers
do iLayer=1,nLayers
! compute liquid water equivalent of total water (liquid plus ice)
if (iLayer>nSnow) then ! soil layer = no volume expansion
iSoil = iLayer - nSnow
vGn_m = 1._rkind - 1._rkind/vGn_n(iSoil)
scalarTheta = mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)
else ! snow layer = volume expansion allowed
iSoil = integerMissing
vGn_m = realMissing
scalarTheta = mLayerVolFracIce(iLayer)*(iden_ice/iden_water) + mLayerVolFracLiq(iLayer)
end if
! *****
! * check that the initial volumetric fraction of liquid water and ice is reasonable...
! *************************************************************************************
select case(mlayerLayerType(iLayer))
! ***** snow
case(iname_snow)
! (check liquid water)
if(mLayerVolFracLiq(iLayer) < 0._rkind)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water < 0: layer = ',iLayer; err=20; return; end if
if(mLayerVolFracLiq(iLayer) > 1._rkind)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water > 1: layer = ',iLayer; err=20; return; end if
! (check ice)
if(mLayerVolFracIce(iLayer) > 0.80_rkind)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice > 0.80: layer = ',iLayer; err=20; return; end if
if(mLayerVolFracIce(iLayer) < 0.05_rkind)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice < 0.05: layer = ',iLayer; err=20; return; end if
! check total water
if(scalarTheta > 0.80_rkind)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > 0.80: layer = ',iLayer; err=20; return; end if
if(scalarTheta < 0.05_rkind)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] < 0.05: layer = ',iLayer; err=20; return; end if
! ***** soil
case(iname_soil)
! (check liquid water)
if(mLayerVolFracLiq(iLayer) < theta_res(iSoil)-xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water < theta_res: layer = ',iLayer; err=20; return; end if
if(mLayerVolFracLiq(iLayer) > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water > theta_sat: layer = ',iLayer; err=20; return; end if
! (check ice)
if(mLayerVolFracIce(iLayer) < 0._rkind )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice < 0: layer = ' ,iLayer; err=20; return; end if
if(mLayerVolFracIce(iLayer) > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice > theta_sat: layer = ',iLayer; err=20; return; end if
! check total water
if(scalarTheta < theta_res(iSoil)-xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] < theta_res: layer = ',iLayer; err=20; return; end if
if(scalarTheta > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > theta_sat: layer = ',iLayer; err=20; return; end if
case default
write(*,*) 'Cannot recognize case in initial vol water/ice check: type=', mlayerLayerType(iLayer)
err=20; message=trim(message)//'cannot identify layer type'; return
end select
! *****
! * check that the initial conditions are consistent with the constitutive functions...
! *************************************************************************************
select case(mLayerLayerType(iLayer))
! ** snow
case(iname_snow)
! check that snow temperature is less than freezing
if(mLayerTemp(iLayer) > Tfreeze)then
message=trim(message)//'initial snow temperature is greater than freezing'
err=20; return
end if
! ensure consistency among state variables
call updateSnow(&
! input
mLayerTemp(iLayer), & ! intent(in): temperature (K)
scalarTheta, & ! intent(in): mass fraction of total water (-)
snowfrz_scale, & ! intent(in): scaling parameter for the snow freezing curve (K-1)
! output
mLayerVolFracLiq(iLayer), & ! intent(out): volumetric fraction of liquid water (-)
mLayerVolFracIce(iLayer), & ! intent(out): volumetric fraction of ice (-)
fLiq, & ! intent(out): fraction of liquid water (-)
err,cmessage) ! intent(out): error control
if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
! ** soil
case(iname_soil)
! ensure consistency among state variables
call updateSoil(&
! input
mLayerTemp(iLayer), & ! intent(in): layer temperature (K)
mLayerMatricHead(iLayer-nSnow), & ! intent(in): matric head (m)
vGn_alpha(iSoil),vGn_n(iSoil),theta_sat(iSoil),theta_res(iSoil),vGn_m, & ! intent(in): van Genutchen soil parameters
! output
scalarTheta, & ! intent(out): volumetric fraction of total water (-)
mLayerVolFracLiq(iLayer), & ! intent(out): volumetric fraction of liquid water (-)
mLayerVolFracIce(iLayer), & ! intent(out): volumetric fraction of ice (-)
err,cmessage) ! intent(out): error control
if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors)
case default; err=10; message=trim(message)//'unknown case for model layer'; return
end select
end do ! (looping through layers)
! end association to variables in the data structures
end associate
! if snow layers exist, compute snow depth and SWE
if(nSnow > 0)then
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%scalarSWE)%dat(1) = sum( (progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)*iden_water + &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow)*iden_ice) * &
progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerDepth)%dat(1:nSnow) )
end if ! if snow layers exist
! check that the layering is consistent
do iLayer=1,nLayers
h1 = sum(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerDepth)%dat(1:iLayer)) ! sum of the depths up to the current layer
h2 = progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(iLayer) - progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(0) ! difference between snow-atm interface and bottom of layer
if(abs(h1 - h2) > 1.e-6_rkind)then
write(message,'(a,1x,i0,a,f5.3,a,f5.3)') trim(message)//'mis-match between layer depth and layer height; layer = ', iLayer, '; sum depths = ',h1,'; height = ',h2
err=20; return
end if
end do
end do ! iHRU
end do ! iGRU
end subroutine check_icond
end module check_icond_module