-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathNCsed_leach.f90
225 lines (194 loc) · 8.41 KB
/
NCsed_leach.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
subroutine orgncswat2(iwave)
!! ~ ~ ~ PURPOSE ~ ~ ~
!! this subroutine calculates the amount of organic nitrogen removed in
!! surface runoff - when using CSWAT==2 it
!! ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! da_ha |ha |area of watershed in hectares
!! enratio |none |enrichment ratio calculated for day in HRU
!! erorgn(:) |none |organic N enrichment ratio, if left blank
!! |the model will calculate for every event
!! ihru |none |HRU number
!! iwave |none |flag to differentiate calculation of HRU and
!! |subbasin sediment calculation
!! |iwave = 0 for HRU
!! |iwave = subbasin # for subbasin
!! sedyld(:) |metric tons |daily soil loss caused by water erosion in
!! |HRU
!! sol_bd(:,:) |Mg/m**3 |bulk density of the soil
!! sol_z(:,:) |mm |depth to bottom of soil layer
!! sub_bd(:) |Mg/m^3 |bulk density in subbasin first soil layer
!! sub_fr(:) |none |fraction of watershed area in subbasin
!! sub_orgn(:) |kg N/ha |amount of nitrogen stored in all organic
!! sedc_d(:) |kg C/ha |amount of C lost with sediment
!!
!!
!!
!! |pools
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! sedorgn(:) |kg N/ha |amount of organic nitrogen in surface runoff
!! |in HRU for the day
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! conc | |concentration of organic N in soil
!! er |none |enrichment ratio
!! j |none |HRU number
!! wt1 |none |conversion factor (mg/kg => kg/ha)
!! xx |kg N/ha |amount of organic N in first soil layer
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
use parm
integer, intent (in) :: iwave
integer :: j
real :: xx, wt1, er, conc
real :: sol_mass, QBC, VBC, YBC, YOC, YW, TOT, YEW, X1, PRMT_21, PRMT_44
real :: DK, V, X3, CO, CS, perc_clyr, latc_clyr
integer :: k
latc_clyr = 0.
j = 0
j = ihru
!!for debug purpose by zhang
!if (iyr == 1991 .and. i==235) then
! write(*,*) 'stop'
!end if
xx = 0.
wt1 = 0. !! conversion factor
er = 0. !! enrichment ratio
if (iwave <= 0) then
!! HRU calculations
!xx = sol_n(1,j) + sol_fon(1,j) + sol_mn(1,j)
xx = sol_LSN(1,j)+sol_LMN(1,j)+sol_HPN(1,j)+sol_HSN(1,j) !+sol_BMN(1,j)
!wt = sol_bd(1,j) * sol_z(1,j) * 10. (tons/ha)
!wt1 = wt/1000
wt1 = sol_bd(1,j) * sol_z(1,j) / 100.
if (erorgn(j) > .001) then
er = erorgn(j)
else
er = enratio
end if
else
!! subbasin calculations
xx = sub_orgn(iwave)
wt1 = sub_bd(iwave) * sol_z(1,j) / 100.
er = enratio
end if
conc = 0.
conc = xx * er / wt1
if (iwave <= 0) then
!! HRU calculations
sedorgn(j) = .001 * conc * sedyld(j) / hru_ha(j)
else
!! subbasin calculations
sedorgn(j) = .001 * conc * sedyld(j) / (da_ha * sub_fr(iwave))
end if
!! update soil nitrogen pools only for HRU calculations
if (iwave <= 0 .and. xx > 1.e-6) then
xx1 = (1. - sedorgn(j) / xx)
!!add by zhang to update soil nitrogen pools
sol_LSN(1,j) = sol_LSN(1,j) * xx1
sol_LMN(1,j) = sol_LMN(1,j) * xx1
sol_HPN(1,j) = sol_HPN(1,j) * xx1
sol_HSN(1,j) = sol_HSN(1,j) * xx1
!sol_BMN(1,j) = sol_BMN(1,j) * xx1
end if
!return
!Calculate runoff and leached C&N from micro-biomass
latc_clyr = 0.
sol_mass = 0.
!kg/ha
sol_mass = (sol_z(1,j) / 1000.) * 10000. * sol_bd(1,j)* 1000. * (1- sol_rock(1,j) / 100.)
QBC=0. !c loss with runoff or lateral flow
VBC=0. !c los with vertical flow
YBC=0. !BMC LOSS WITH SEDIMENT
YOC=0. !Organic C loss with sediment
YW=0. !YW = WIND EROSION (T/HA)
TOT=sol_HPC(1,j)+sol_HSC(1,j)+sol_LMC(1,j)+sol_LSC(1,j) !Total organic carbon in layer 1
!YEW = MIN(er*(sedyld(j)/hru_ha(j)+YW/hru_ha(j))/(sol_mass/1000.),.9)
! Not sure whether should consider enrichment ratio or not!
YEW = MIN((sedyld(j)/hru_ha(j)+YW/hru_ha(j))/(sol_mass/1000.),.9) !fraction of soil erosion of total soil mass
X1=1.-YEW
!YEW=MIN(ER*(YSD(NDRV)+YW)/WT(LD1),.9)
!ER enrichment ratio
!YSD water erosion
!YW wind erosion
YOC=YEW*TOT
sol_HSC(1,j)=sol_HSC(1,j)*X1
sol_HPC(1,j)=sol_HPC(1,j)*X1
sol_LS(1,j)=sol_LS(1,j)*X1
sol_LM(1,j)=sol_LM(1,j)*X1
sol_LSL(1,j)=sol_LSL(1,j)*X1
sol_LSC(1,j)=sol_LSC(1,j)*X1
sol_LMC(1,j)=sol_LMC(1,j)*X1
sol_LSLC(1,j)=sol_LSLC(1,j)*X1
sol_LSLNC(1,j)=sol_LSC(1,j)-sol_LSLC(1,j)
if (surfq(j) > 0) then
!write(*,*) 'stop'
end if
IF(sol_BMC(1,j)>.01) THEN
PRMT_21 = 0. !KOC FOR CARBON LOSS IN WATER AND SEDIMENT(500._1500.) KD = KOC * C
PRMT_21 = 1000.
sol_WOC(1,j) = sol_LSC(1,j)+sol_LMC(1,j)+sol_HPC(1,j)+sol_HSC(1,j)+sol_BMC(1,j)
DK=.0001*PRMT_21*sol_WOC(1,j)
!X1=PO(LD1)-S15(LD1)
X1 = sol_por(1,j)*sol_z(1,j)-sol_wpmm(1,j) !mm
IF (X1 <= 0.) THEN
X1 = 0.01
END IF
XX=X1+DK
!V=QD+Y4
V = surfq(j) + sol_prk(1,j) + flat(1,j)
!QD surface runoff
X3=0.
IF(V>1.E-10)THEN
X3=sol_BMC(1,j)*(1.-EXP(-V/XX)) !loss of biomass C
PRMT_44 = 0. !RATIO OF SOLUBLE C CONCENTRATION IN RUNOFF TO PERCOLATE(0.1_1.)
PRMT_44 = .5
CO=X3/(sol_prk(1,j) + PRMT_44*(surfq(j)+flat(1,j))) !CS is the horizontal concentration
CS=PRMT_44*CO !CO is the vertical concentration
VBC=CO*(sol_prk(1,j)) !!! sol_prk(:,:) |mm H2O |percolation from soil layer on current day
sol_BMC(1,j)=sol_BMC(1,j)-X3
QBC=CS*(surfq(j)+flat(1,j))
! COMPUTE WBMC LOSS WITH SEDIMENT
IF(YEW>0.)THEN
CS=DK*sol_BMC(1,j)/XX
YBC=YEW*CS
END IF
END IF
END IF
sol_BMC(1,j)=sol_BMC(1,j)-YBC
surfqc_d(j) = QBC*(surfq(j)/(surfq(j)+flat(1,j)+1.e-6))
sol_latc(1,j) = QBC*(flat(1,j)/(surfq(j)+flat(1,j)+1.e-6))
sol_percc(1,j) = VBC
sedc_d(j) = YOC + YBC
latc_clyr = latc_clyr + sol_latc(1,j)
DO k=2,sol_nly(j)
if (sol_prk(k,j) > 0 .and. k == sol_nly(j)) then
!write (*,*) 'stop'
end if
sol_thick = 0.
sol_thick = sol_z(k,j)-sol_z(k-1,j)
sol_WOC(k,j) = sol_LSC(k,j)+sol_LMC(k,j)+sol_HPC(k,j)+sol_HSC(k,j)
Y1=sol_BMC(k,j)+VBC
VBC=0.
IF(Y1>=.01)THEN
V=sol_prk(k,j) + flat(k,j)
IF(V>0.)VBC=Y1*(1.-EXP(-V/(sol_por(k,j)*sol_thick-sol_wpmm(k,j)+.0001*PRMT_21*sol_WOC(k,j))))
END IF
sol_latc(k,j) = VBC*(flat(k,j)/(sol_prk(k,j) + flat(k,j)+1.e-6))
sol_percc(k,j) = VBC-sol_latc(k,j)
sol_BMC(k,j)=Y1-VBC
!! calculate nitrate in percolate
!perc_clyr = 0.
perc_clyr = perc_clyr + sol_percc(k,j)
latc_clyr = latc_clyr + sol_latc(k,j)
END DO
latc_d(j) = latc_clyr
percc_d(j) = perc_clyr
return
end