-
Notifications
You must be signed in to change notification settings - Fork 1
/
S3M_Module_Phys_Snow_Apps_Hydraulics.f90
165 lines (125 loc) · 8.13 KB
/
S3M_Module_Phys_Snow_Apps_Hydraulics.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
!------------------------------------------------------------------------------------
! File: S3M_Module_Phys_Snow_Apps_Hydraulics.f90
! Author: Francesco Avanzi
!
! Created on June 10, 2020 5:00 PM
! Last update on October 27, 2020 10:45 AM
!
! Snow-hydraulics application module
!------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------
! Module Header
module S3M_Module_Phys_Snow_Apps_Hydraulics
!------------------------------------------------------------------------------------
! External module(s)
use S3M_Module_Namelist, only: oS3M_Namelist
use S3M_Module_Vars_Loader, only: oS3M_Vars
use S3M_Module_Tools_Debug
! Implicit none for all subroutines in this module
implicit none
!------------------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------------------
! Subroutine to compute snowpack outflow based on Darcy's
!(Avanzi et al. 2015, https://doi.org/10.1016/j.advwatres.2015.09.021)
subroutine S3M_Phys_Snow_Apps_Outflow(iID, iRows, iCols, iDtForcing, &
dVarRhoW, &
a2dVarDem, &
a2dVarSWE_D, a2dVarRho_D, a2dVarSWE_W, &
a2dVarRainfall, a2dVarMeltingS, a2dVarRefreezingS, &
a2dVarOutflow_K, a2dVarH_D, a2dVarH_S, a2dVarSSA, a2dVarPerm, a2dVarTheta_W)
!------------------------------------------------------------------------------------------
! Variable(s) declaration
integer(kind = 4) :: iID, iRows, iCols, iDtForcing
real(kind = 4) :: dVarRhoW
real(kind = 4), dimension(iRows, iCols) :: a2dVarDem
real(kind = 4), dimension(iRows, iCols) :: a2dVarTheta_W, a2dVarSWE_D, a2dVarRho_D, a2dVarSWE_W
real(kind = 4), dimension(iRows, iCols) :: a2dVarOutflow_K, a2dVarRainfall, a2dVarMeltingS, a2dVarRefreezingS
real(kind = 4), dimension(iRows, iCols) :: a2dVarPorosity, a2dVarH_D, a2dVarH_S
real(kind = 4), dimension(iRows, iCols) :: a2dVarSr, a2dVarSr_irr, a2dVarSr_star
real(kind = 4), dimension(iRows, iCols) :: a2dVarSSA, a2dVarr_e, a2dVarPerm, a2dVarCond
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Initialize variable(s)
a2dVarPorosity = -9999; a2dVarSr = -9999; a2dVarSr_irr = -9999; a2dVarSr_star = -9999; a2dVarr_e = -9999;
a2dVarCond = -9999;
! Info start
call mprintf(.true., iINFO_Extra, ' Phys :: Snow :: Outflow ... ' )
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Compute porosity
where( (a2dVarDem.ge.0.0) .and. (a2dVarH_D.gt.0.0))
a2dVarPorosity = 1 - a2dVarRho_D/917
elsewhere(a2dVarDem.ge.0.0)
a2dVarPorosity = 0.0
endwhere
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Compute dry-snow depth
where( (a2dVarDem.ge.0.0) .and. (a2dVarRho_D.gt.0.0) )
a2dVarH_D = ((a2dVarSWE_D/1000)*dVarRhoW)/a2dVarRho_D
elsewhere
a2dVarH_D = 0.0
endwhere
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Compute control volume and saturation degree
where( (a2dVarDem.ge.0.0) .and. (a2dVarH_D.gt.0.0) &
.and. ((a2dVarSWE_W/1000 - (a2dVarPorosity*a2dVarH_D)) .ge. 0.0) ) !we divide by 1000 to convert mm to m
!and so compare a2dVarSWE_W with a2dVarH_D
a2dVarH_S = a2dVarH_D + ((a2dVarSWE_W/1000) - (a2dVarPorosity*a2dVarH_D))
a2dVarSr = 1.0
elsewhere( (a2dVarDem.ge.0.0) .and. (a2dVarH_D.gt.0.0) )
a2dVarH_S = a2dVarH_D
a2dVarSr = (a2dVarSWE_W/1000)/(a2dVarPorosity*a2dVarH_D)
elsewhere( (a2dVarDem.ge.0.0) )
a2dVarH_S = 0.0
a2dVarSr = 0.0
endwhere
where( (a2dVarDem.ge.0.0) .and. (a2dVarPorosity.gt.0.0))
a2dVarSr_irr = 0.02*((a2dVarRho_D/dVarRhoW)/a2dVarPorosity)
a2dVarSr_star = (a2dVarSr - a2dVarSr_irr)/(1 - a2dVarSr_irr)
elsewhere
a2dVarSr_irr = 0.0
a2dVarSr_star = 0.0
endwhere
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Compute SSA, r_e, permeability, and conductivity
where( (a2dVarDem.ge.0.0) .and. (a2dVarRho_D.gt.0.0))
a2dVarSSA = -308.2*LOG(a2dVarRho_D/1000) - 206 ! SSA here will be in cm2/g as in Domine et al 07.
! a2dVarRho_D must be in g/cm3
a2dVarSSA = a2dVarSSA/10 !now a2dVarSSA is in m2/kg
a2dVarr_e = 3/(a2dVarSSA*917)
a2dVarPerm = 3*(a2dVarr_e**2)*EXP((-0.013*a2dVarRho_D))
a2dVarCond = a2dVarPerm*(a2dVarSr_star**3)
endwhere
where(a2dVarSr .lt. a2dVarSr_irr) a2dVarCond = 0.0
! Here we force a2dVarCond to 0 where a2dVarSr lt a2dVarSr_irr; the implication is the a2dVarOutflow will be 0
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Compute a2dVarOutflow_K, but some conditions apply....
! Note that a2dVarRainfall has been already added to SWE_W at this point.
where( (a2dVarSr.ge.0.5) .or. (a2dVarSWE_D .lt. 10) )
a2dVarOutflow_K = a2dVarSWE_W
! This is a pragmatic solution that is particularly helpful given that S3M is employed for flood forecasting:
! wherever saturation is extremely high OR the snowpack is very shallow, we assume snow has no retention capability
! and all wet SWE (here possibly including a2dVarRainfall already, although in principle it should
! have been allocated to a2dVarOutflow_ExcessRain) is immediately discharged.
elsewhere(a2dVarSWE_W .le. (5.47*(10**5)*a2dVarCond*(iDtForcing)*1000))
! 5.47*(10**5) includes viscosity and density.
a2dVarOutflow_K = a2dVarSWE_W
elsewhere
a2dVarOutflow_K = 5.47*(10**5)*a2dVarCond*(iDtForcing)*1000
! 5.47*(10**5) includes viscosity and density.
endwhere
where( a2dVarOutflow_K.le.0.0 ) a2dVarOutflow_K = 0.0 !just in case....
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
! Info end
call mprintf(.true., iINFO_Extra, ' Phys :: Snow :: Phase Part OUtflow... OK' )
!------------------------------------------------------------------------------------------
end subroutine S3M_Phys_Snow_Apps_Outflow
!------------------------------------------------------------------------------------------
end module S3M_Module_Phys_Snow_Apps_Hydraulics
!------------------------------------------------------------------------------------