-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcarbon_zhang2.f90
819 lines (766 loc) · 34.9 KB
/
carbon_zhang2.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
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
subroutine carbon_zhang2
use parm
!!============================================
!!Input variables
!! sol_bd(:,:) |Mg/m**3 |bulk density of the soil
!! sol_st(:,:) |mm H2O |amount of water stored in the soil layer on
!! |current day
!! sol_fc(:,:) |mm H2O |amount of water available to plants in soil
!! |layer at field capacity (fc - wp),Index:(layer,HRU)
!! sol_wp(:,:) |mm H20/mm soil|water content of soil at -1.5 MPa (wilting
!! |point)
!! sol_wpmm(:,:) |mm H20 |water content of soil at -1.5 MPa (wilting
!! |point)
!! sol_percc(k,j)
!! sol_latc(k,j)
!!==============================================
!!Transput variables;
!! sol_HSC(:,:) : mass of C present in slow humus (kg ha-1)
!! sol_HSN(:,:) : mass of N present in slow humus (kg ha-1)
!! sol_HPC(:,:) : mass of C present in passive humus (kg ha-1)
!! sol_HPN(:,:) : mass of N present in passive humus (kg ha-1)
!! sol_LM(:,:) : mass of metabolic litter (kg ha-1)
!! sol_LMC(:,: : mass of C in metabolic litter (kg ha-1)
!! sol_LMN(:,:) : mass of N in metabolic litter (kg ha-1)
!! sol_LS(:,:) : mass of structural litter (kg ha-1)
!! sol_LSC(:,:) : mass of C in structural litter (kg ha-1)
!! sol_LSL(:,:) : mass of lignin in structural litter (kg ha-1)
!! sol_LSN(:,:) : mass of N in structural litter (kg ha-1)
!! STD(:) : standing dead (kg ha-1) (Not used)
!! STDL(:) : mass of lignin in standing dead (kg ha-1) (Not used)
!! STDN(:) : mass of N in standing dead (dead plants + sorbed from soil; kg ha-1) (Not used)
!! STDNEl(:) : standing dead N after enrichment with sorbed N in a soil layer (kg ha-1)
!! sol_NO3(:,:) : weight of NO3-N in soil layer (kg ha-1)
!! sol_NH3(:,:) : weight of NH3-N in soil layer (kg ha-1)
!!==============================================
!!read in parameters
!HPR : rate of transformation of passive humus under optimal conditions (subsurface
!layers = 0.000012 day-1) (Parton et al.,1993, 1994)
!HSR : rate of transformation of slow humus under optimal conditions (all layers
!= 0.0005 day-1) (Parton et al., 1993, 1994; Vitousek et al., 1993)
!KOC : liquid–solid partition coefficient for microbial biomass (10^3 m^3 Mg-1)
!BMR : rate of transformation of microbial biomass and associated products under optimal
! conditions (surface = 0.0164 day-1; all other layers = 0.02 day-1) (Parton et al., 1993, 1994)
!LMR : rate of transformation of metabolic litter under optimal conditions (surface =
!0.0405 day-1; all other layers = 0.0507 day-1) (Parton et al., 1994)
!Sf : fraction of mineral N sorbed to litter: 0.05 for surface litter, 0.1 for belowground litter
!Cf : carbon fraction of organic materials 0.42; from data of Pinck et al., 1950)
!LSR : rate of potential transformation of structural litter under optimal conditions
!(surface = 0.0107 day-1; all other layers= 0.0132 day-1) (Parton et al., 1994)
!PRMT_51 !COEF ADJUSTS MICROBIAL ACTIVITY FUNCTION IN TOP SOIL LAYER (0.1_1.),
!PRMT_45 !COEF IN CENTURY EQ ALLOCATING SLOW TO PASSIVE HUMUS(0.001_0.05) ORIGINAL VALUE = 0.003,
!!==============================================
!! local variables
!ABCO2 : allocation from biomass to CO2; 0.6 (surface litter), 0.85–0.68*(CLAF + SILF) (all other layers) (Parton et al., 1993, 1994)
!ABL : carbon allocation from biomass to leaching; ABL = (1-exp(-f/(0.01* SW+ 0.1*(KdBM)*DB)) (Williams, 1995)
!ABP : allocation from biomass to passive humus; 0 (surface litter), 0.003 + 0.032*CLAF (all other layers) (Parton et al., 1993, 1994)
!ALMCO2 : allocation from metabolic litter to CO2; 0.6 (surface litter), 0.55 (all other layers) (Parton et al., 1993, 1994)
!ALSLCO2 : allocation from lignin of structural litter to CO2; 0.3 (Parton et al., 1993, 1994)
!ALSLNCO2: allocation from non-lignin of structural litter to CO2; 0.6 (surface litter), 0.55 (all other layers) (Parton et al., 1993, 1994)
!APCO2 : allocation from passive humus to CO2; 0.55 (Parton et al., 1993, 1994)
!ASCO2 : allocation from slow humus to CO2; 0.55 (Parton et al., 1993, 1994)
!ASP : allocation from slow humus to passive; 0 (surface litter), 0.003-0.009*CLAF (all other layers) (Parton et al., 1993, 1994)
!BMC : mass of C in soil microbial biomass and associated products (kg ha-1)
!BMCTP : potential transformation of C in microbial biomass (kg ha-1 day-1)
!BMN : mass of N in soil microbial biomass and associated products (kg ha-1)
!BMNTP : potential transformation of N in microbial biomass (kg ha-1 day-1)
!CDG : soil temperature control on biological processes
!CNR : C/N ratio of standing dead
!CS : combined factor controlling biological processes [CS = sqrt(CDG×SUT)* 0.8*OX*X1), CS < 10; CS = 10, CS>=10 (Williams, 1995)]
!DBp : soil bulk density of plow layer (Mg m-3) (Not used)
!HSCTP : potential transformation of C in slow humus (kg ha-1 day-1)
!HSNTP : potential transformation of N in slow humus (kg ha-1 day-1)
!HPCTP : potential transformation of C in passive humus (kg ha-1 day-1)
!HPNTP : potential transformation of N in passive humus (kg ha-1 day-1)
!LMF : fraction of the litter that is metabolic
!LMNF : fraction of metabolic litter that is N (kg kg-1)
!LMCTP : potential transformation of C in metabolic litter (kg ha-1 day-1)
!LMNTP : potential transformation of N in metabolic litter (kg ha-1 day-1)
!LSCTP : potential transformation of C in structural litter (kg ha-1 day-1)
!LSF : fraction of the litter that is structural
!LSLF : fraction of structural litter that is lignin (kg kg-1)
!LSNF : fraction of structural litter that is N (kg kg-1)
!LSLCTP : potential transformation of C in lignin of structural litter (kg ha-1 day-1)
!LSLNCTP : potential transformation of C in nonlignin structural litter (kg ha-1 day-1)
!LSNTP : potential transformation of N in structural litter (kg ha-1 day-1)
!NCBM : N/C ratio of biomass
!NCHP : N/C ratio passive humus
!NCHS : N/C ratio of the slow humus
!OX : oxygen control on biological processes with soil depth
!SUT : soil water control on biological processes
!X1 : tillage control on residue decomposition (Not used)
!XBMT : control on transformation of microbial biomass by soil texture and structure.
!Its values: surface litter layer = 1; all other layers = 1-0.75*(SILF + CLAF) (Parton et al., 1993, 1994)
!XLSLF : control on potential transformation of structural litter by lignin fraction
!of structural litter [XLSLF = exp(-3* LSLF) (Parton et al., 1993, 1994)]
integer :: j, k, kk
real :: sol_mass, sol_min_n
real :: fc, wc, sat, void, sut, cdg, OX, CS
real :: X1,X3, XX
real :: LMF, LSF, LSLF, XLSLF, LSR, BMR, XBMT, HSR, HPR
real :: LSCTA, LSLCTA, LSLNCTA,LSNTA, LMCTA, LMNTA, BMCTA, BMNTA, HSCTA, HSNTA, HPCTA, HPNTA
real :: LSCTP, LSLCTP, LSLNCTP, LSNTP, LMR, LMCTP, LMNTP, BMCTP,HSCTP, HSNTP, HPCTP, HPNTP
real :: NCHP, Nf, NCBM, NCHS, ALSLCO2, ALSLNCO2,ALMCO2,ABCO2, A1CO2, APCO2, ASCO2, ABP, ASP, A1, ASX, APX
real :: PRMT_51, PRMT_45
real :: DF1, DF2, SNMN, DF3, DF4, DF5, DF6, ADD, ADF1, ADF2, ADF3, ADF4, ADF5
real :: TOT
real :: PN1, PN2, PN3, PN4, PN5, PN6, PN7, PN8, PN9
real :: SUM, CPN1, CPN2, CPN3, CPN4, CPN5
real :: WMIN,DMDN, wdn, Delta_BMC, DeltaWN
!! initilize local variables
DeltaWN = 0.
DeltaBMC = 0.
wdn = 0.
X1 = 0.
X3 = 0.
XX = 0.
fc = 0.
wc = 0.
sat = 0.
void = 0.
sut = 0.
cdg = 0.
OX = 0.
CS = 0.
LMF = 0.
LSF = 0.
LSLF = 0.
XLSLF = 0.
LSR = 0.
LSCTP = 0.
LSCTA = 0.
LSLCTA = 0.
LSLNCTA = 0.
SNTA = 0.
LMCTA = 0.
LMNTA = 0.
BMCTA = 0.
BMNTA = 0.
HSCTA = 0.
HSNTA= 0.
HPCTA= 0.
HPNTA= 0.
LSLCTP= 0.
LSLNCTP= 0.
LSNTP= 0.
LMR= 0.
LMCTP= 0.
LMNTP= 0.
BMR= 0.
XBMT= 0.
BMCTP= 0.
HSR= 0.
HSCTP= 0.
HSNTP= 0.
HPR= 0.
HPCTP= 0.
HPNTP= 0.
NCHP= 0.
Nf= 0.
NCBM= 0.
NCHS= 0.
ALSLCO2= 0.
ALSLNCO2= 0.
ALMCO2= 0.
ABCO2= 0.
A1CO2= 0.
APCO2= 0.
ASCO2= 0.
ABP= 0.
ASP= 0.
A1= 0.
ASX= 0.
APX= 0.
DF1= 0.
DF2= 0.
SNMN= 0.
DF3= 0.
DF4= 0.
DF5= 0.
DF6= 0.
ADD= 0.
ADF1= 0.
ADF2= 0.
ADF3= 0.
ADF4= 0.
ADF5= 0.
PN1= 0.
PN2= 0.
PN3= 0.
PN4= 0.
PN5= 0.
PN6= 0.
PN7= 0.
PN8= 0.
PN9= 0.
TOT= 0.
SUM= 0.
CPN1= 0.
CPN2= 0.
CPN3= 0.
CPN4= 0.
CPN5= 0.
WMIN= 0.
DMDN= 0.
j=0
j=ihru
do k = 1,sol_nly(j)
! a simple equation to calculate Bulk Density from DSSAT (Not Used)
!XZ = sol_cbn(k,j)
!sol_BDM(k)=ZZ/(1./BD(J)-XZ/.224)
end do
!!for debug purpose by zhang
if (iyr == 1941 .and. i==134) then
!write(*,*) 'stop'
end if
!calculate tillage factor using DSSAT
if (tillage_switch(j) .eq. 1 .and. tillage_days(j) .le. 30) then
tillage_factor(j) = 1.6
else
tillage_factor(j) = 1.0
end if
!calculate tillage factor using DSSAT
!!calculate C/N dynamics for each soil layer
!!===========================================
do k = 1, sol_nly(j)
if (k == 1) then
!10 cm / 1000 = 0.01m; 1 ha = 10000 m2; ton/m3; * 1000 --> final unit is kg/ha; rock fraction is considered
sol_mass = (10) / 1000.* 10000. * sol_bd(k,j)* 1000. *(1- sol_rock(k,j) / 100.)
else
sol_mass = (sol_z(k,j) - sol_z(k-1,j)) / 1000.* 10000. * sol_bd(k,j)* 1000. *(1- sol_rock(k,j) / 100.)
end if
!!If k = 1, then using temperature, soil moisture in layer 2 to calculate decomposition factor
!!Not
kk =0
if (k == 1) then
kk = 2
else
kk = k
end if
!! mineralization can occur only if temp above 0 deg
!check sol_st soil water content in each soil ayer mm H2O
if (sol_tmp(k,j) > 0. .AND. sol_st(k,j) > 0.) then
!!from Armen
!!compute soil water factor - sut
fc = sol_fc(k,j) + sol_wpmm(k,j) ! units mm
wc = sol_st(k,j) + sol_wpmm(k,j) ! units mm
sat = sol_ul(k,j) + sol_wpmm(k,j) ! units mm
void = sol_por(k,j) * (1. - wc / sat) ! fraction
!!from Armen
sut = 0.
!sut = .1 + .9 * Sqrt(sol_st(kk,j) / sol_fc(kk,j))
!sut = .1 + .9 * Sqrt(wc / fc)
X1=wc-sol_wpmm(k,j)
IF(X1<0.)THEN
SUT=.1*(sol_st(kk,j) /sol_wpmm(k,j))**2
ELSE
SUT = .1 + .9 * Sqrt(sol_st(k,j) / sol_fc(k,j))
END IF
sut = Min(1., sut)
sut = Max(.05, sut)
!check X1, FC, S15
!from Armen
!wf = fwf(fc,wc,sol_wpmm(kk,j))
!of = fof(void,sol_por(kk,j))
!sut = wf * of
!from Armen
!!compute tillage factor (X1)
!use the tillfactor module from Armen
!X1 = ftilf(tillagef(kk,j), wc, sat)
X1 = 1.0
!calculate tillage factor using DSSAT
if (tillage_switch(j) .eq. 1 .and. tillage_days(j) .le. 30) then
if (k == 1) then
X1 = 1.6
else
if (sol_z(k,j) .le. tillage_depth(j)) then
X1 = 1.6
elseif (sol_z(k-1,j) .lt. tillage_depth(j)) then
X1 = 1.0 + 0.6*(tillage_depth(j) - sol_z(k-1,j))/(sol_z(k,j) - sol_z(k-1,j))
end if
end if
else
X1 = 1.0
end if
!calculate tillage factor using DSSAT
!!compute soil temperature factor
!!When sol_tep is larger than 35, cdg is negative?
cdg = 0.
!if (sol_tmp(kk,j) <= 35.0) then
!cdg = sol_tmp(kk,j) / (sol_tmp(kk,j) + exp(5.058 - 0.2504 * sol_tmp(kk,j)))
cdg = sol_tmp(k,j) / (sol_tmp(k,j) + exp(5.058459 - 0.2503591 * sol_tmp(k,j)))
!end if
!!from Armen
cdg = fcgd(sol_tmp(k,j))
!!from Armen
!!compute oxygen (OX)
OX = 0.
!OX = 1 - (0.9* sol_z(k,j)/1000.) / (sol_z(k,j)/1000.+ exp(1.50-3.99*sol_z(k,j)/1000.))
!OX = 1 - (0.8* sol_z(k,j)) / (sol_z(k,j)+ exp(1.50-3.99*sol_z(k,j)))
OX=1.-0.8*((sol_z(kk,j)+sol_z(kk-1,j))/2)/(((sol_z(kk,j)+sol_z(kk-1,j))/2)+EXP(18.40961-0.023683632*((sol_z(kk,j)+sol_z(kk-1,j))/2)))
!! compute combined factor
CS = 0.
CS=MIN(10.,SQRT(cdg*sut)*0.9*OX*X1)
!! call denitrification (to use void and cdg factor)
wdn = 0.
cdg = fcgd(sol_tmp(k,j))
if (cdg > 0. .and. void <= 0.1) then
call ndenit(k,j,cdg,wdn,void)
end if
wshd_dnit = wshd_dnit + wdn * hru_dafr(j)
wdntl = wdntl + wdn
sol_min_n = sol_no3(k,j) + sol_nh3(k,j)
!lignin content in structural litter (fraction)
RLR = 0.
RLR = min(0.8,sol_LSL(k,j)/(sol_LS(k,j) + 1.E-5))
!HSR=PRMT(47) !CENTURY SLOW HUMUS TRANSFORMATION RATE D^-1(0.00041_0.00068) ORIGINAL VALUE = 0.000548,
HSR = 5.4799998E-04
!HPR=PRMT(48) !CENTURY PASSIVE HUMUS TRANSFORMATION RATE D^-1(0.0000082_0.000015) ORIGINAL VALUE = 0.000012
HPR = 1.2000000E-05
APCO2=.55
ASCO2=.60
PRMT_51 =0. !COEF ADJUSTS MICROBIAL ACTIVITY FUNCTION IN TOP SOIL LAYER (0.1_1.),
PRMT_51 = 1.
!!The following codes are clculating of the N:C ration in the newly formed SOM for each pool
!!please note that in the surface layer, no new materials enter Passive pool, therefore, no NCHP is
!!calculated for the first layer.
IF(k==1)THEN
CS=CS*PRMT_51
ABCO2=.55
A1CO2=.55
BMR=.0164
LMR=.0405
LSR=.0107
NCHP=.1
XBM=1.
! COMPUTE N/C RATIOS
!X1=.1*(WLMN(LD1)+WLSN(LD1))/(RSD(LD1)+1.E-5)
X1 = 0.1*(sol_LSN(k,j)+sol_LMN(k,j))/(sol_rsd(k,j)/1000+1.E-5) !relative notrogen content in residue (%)
IF(X1>2.)THEN
NCBM=.1
GO TO 6
END IF
IF(X1>.01)THEN
NCBM=1./(20.05-5.0251*X1)
ELSE
NCBM=.05
END IF
6 NCHS=NCBM/(5.*NCBM+1.)
GO TO 2
END IF
!ABCO2=.17+.0068*SAN(ISL)
ABCO2 = 0.17 + 0.0068 * sol_sand(k,j)
A1CO2=.55
BMR=.02
LMR=.0507
LSR=.0132
XBM=.25+.0075*sol_sand(k,j)
!X1=1000.*(WNH3(ISL)+WNO3(ISL))/WT(ISL)
!WT is soil mass in tons
X1 = 1000. * sol_min_n/(sol_mass/1000)
IF(X1>7.15)THEN
NCBM=.33
NCHS=.083
NCHP=.143
ELSE
NCBM=1./(15.-1.678*X1)
NCHS=1./(20.-1.119*X1)
NCHP=1./(10.-.42*X1)
END IF
2 ABP=.003+.00032*sol_clay(k,j)
!SMS(3,ISL)=SMS(3,ISL)+CS
PRMT_45 = 0. !COEF IN CENTURY EQ ALLOCATING SLOW TO PASSIVE HUMUS(0.001_0.05) ORIGINAL VALUE = 0.003,
PRMT_45 = 5.0000001E-02
ASP=MAX(.001,PRMT_45-.00009*sol_clay(k,j))
! POTENTIAL TRANSFORMATIONS STRUCTURAL LITTER
X1=LSR*CS*EXP(-3.*RLR)
LSCTP=X1*sol_LSC(k,j)
LSLCTP=LSCTP*RLR
LSLNCTP=LSCTP*(1.-RLR)
LSNTP=X1*sol_LSN(k,j)
! POTENTIAL TRANSFORMATIONS METABOLIC LITTER
X1=LMR*CS
LMCTP=sol_LMC(k,j)*X1
LMNTP=sol_LMN(k,j)*X1
! POTENTIAL TRANSFORMATIONS MICROBIAL BIOMASS
X1=BMR*CS*XBM
BMCTP=sol_BMC(k,j)*X1
BMNTP=sol_BMN(k,j)*X1
! POTENTIAL TRANSFORMATIONS SLOW HUMUS
X1=HSR*CS
HSCTP=sol_HSC(k,j)*X1
HSNTP=sol_HSN(k,j)*X1
! POTENTIAL TRANSFORMATIONS PASSIVE HUMUS
X1=CS*HPR
HPCTP=sol_HPC(k,j)*X1
HPNTP=sol_HPN(k,j)*X1
! ESTIMATE N DEMAND
A1=1.-A1CO2
ASX=1.-ASCO2-ASP
APX=1.-APCO2
PN1=LSLNCTP*A1*NCBM !Structural Litter to Biomass
PN2=.7*LSLCTP*NCHS !Structural Litter to Slow
PN3=LMCTP*A1*NCBM !Metabolic Litter to Biomass
!PN4=BMCTP*ABL*NCBM !Biomass to Leaching (calculated in NCsed_leach)
PN5=BMCTP*ABP*NCHP !Biomass to Passive
PN6=BMCTP*(1.-ABP-ABCO2)*NCHS !Biomass to Slow
PN7=HSCTP*ASX*NCBM !Slow to Biomass
PN8=HSCTP*ASP*NCHP !Slow to Passive
PN9=HPCTP*APX*NCBM !Passive to Biomass
!PN1=LSLNCTP*A1*NCBM
!PN2=.7*LSLCTP*NCHS
!PN3=LMCTP*A1*NCBM
!PN5=BMCTP*ABP*NCHP
!PN6=BMCTP*(1.-ABP-ABCO2)*NCHS
!PN7=HSCTP*ASX*NCBM
!PN8=HSCTP*ASP*NCHP
!PN9=HPCTP*APX*NCBM
! COMPARE SUPPLY AND DEMAND FOR N
SUM=0.
CPN1=0.
CPN2=0.
CPN3=0.
CPN4=0.
CPN5=0.
X1=PN1+PN2
IF(LSNTP<X1)THEN
CPN1=X1-LSNTP
ELSE
SUM=SUM+LSNTP-X1
END IF
IF(LMNTP<PN3)THEN
CPN2=PN3-LMNTP
ELSE
SUM=SUM+LMNTP-PN3
END IF
X1=PN5+PN6
IF(BMNTP<X1)THEN
CPN3=X1-BMNTP
ELSE
SUM=SUM+BMNTP-X1
END IF
X1=PN7+PN8
IF(HSNTP<X1)THEN
CPN4=X1-HSNTP
ELSE
SUM=SUM+HSNTP-X1
END IF
IF(HPNTP<PN9)THEN
CPN5=PN9-HPNTP
ELSE
SUM=SUM+HPNTP-PN9
END IF
! WNH3(ISL)=WNH3(ISL)+SUM
!total available N
WMIN=MAX(1.E-5,sol_NO3(k,j) + sol_NH3(k,j)+SUM)
!WMIN=MAX(1.E-5,sol_NO3(k,j) +SUM)
!total demand for potential tranformaiton of SOM
DMDN=CPN1+CPN2+CPN3+CPN4+CPN5
X3=1.
! REDUCE DEMAND IF SUPPLY LIMITS
IF(WMIN<DMDN) then
X3=WMIN/DMDN
end if
!SMS(5,ISL)=SMS(5,ISL)+X3
! ACTUAL TRANSFORMATIONS
IF(CPN1>0.)THEN
LSCTA=LSCTP*X3
LSNTA=LSNTP*X3
LSLCTA=LSLCTP*X3
LSLNCTA=LSLNCTP*X3
ELSE
LSCTA=LSCTP
LSNTA=LSNTP
LSLCTA=LSLCTP
LSLNCAT=LSLNCTP
END IF
IF(CPN2>0.)THEN
LMCTA=LMCTP*X3
LMNTA=LMNTP*X3
ELSE
LMCTA=LMCTP
LMNTA=LMNTP
END IF
IF(CPN3>0.)THEN
BMCTA=BMCTP*X3
BMNTA=BMNTP*X3
ELSE
BMCTA=BMCTP
BMNTA=BMNTP
END IF
IF(CPN4>0.)THEN
HSCTA=HSCTP*X3
HSNTA=HSNTP*X3
ELSE
HSCTA=HSCTP
HSNTA=HSNTP
END IF
IF(CPN5>0.)THEN
HPCTA=HPCTP*X3
HPNTA=HPNTP*X3
ELSE
HPCTA=HPCTP
HPNTA=HPNTP
END IF
!Recalculate demand using actural transformations
!revised from EPIC code by Zhang
PN1=LSLNCTA*A1*NCBM !Structural Litter to Biomass
PN2=.7*LSLCTA*NCHS !Structural Litter to Slow
PN3=LMCTA*A1*NCBM !Metabolic Litter to Biomass
!PN4=BMCTP*ABL*NCBM !Biomass to Leaching (calculated in NCsed_leach)
PN5=BMCTA*ABP*NCHP !Biomass to Passive
PN6=BMCTA*(1.-ABP-ABCO2)*NCHS !Biomass to Slow
PN7=HSCTA*ASX*NCBM !Slow to Biomass
PN8=HSCTA*ASP*NCHP !Slow to Passive
PN9=HPCTA*APX*NCBM !Passive to Biomass
! COMPARE SUPPLY AND DEMAND FOR N
SUM=0.
CPN1=0.
CPN2=0.
CPN3=0.
CPN4=0.
CPN5=0.
X1=PN1+PN2
IF(LSNTA<X1)THEN
CPN1=X1-LSNTA
ELSE
SUM=SUM+LSNTA-X1
END IF
IF(LMNTA<PN3)THEN
CPN2=PN3-LMNTA
ELSE
SUM=SUM+LMNTA-PN3
END IF
X1=PN5+PN6
IF(BMNTA<X1)THEN
CPN3=X1-BMNTA
ELSE
SUM=SUM+BMNTA-X1
END IF
X1=PN7+PN8
IF(HSNTA<X1)THEN
CPN4=X1-HSNTA
ELSE
SUM=SUM+HSNTA-X1
END IF
IF(HPNTA<PN9)THEN
CPN5=PN9-HPNTA
ELSE
SUM=SUM+HPNTA-PN9
END IF
! WNH3(ISL)=WNH3(ISL)+SUM
!total available N
WMIN=MAX(1.E-5,sol_NO3(k,j) + sol_NH3(k,j)+SUM)
!WMIN=MAX(1.E-5,sol_NO3(k,j) +SUM)
!total demand for potential tranformaiton of SOM
DMDN=CPN1+CPN2+CPN3+CPN4+CPN5
!DMDN=DMDN*X3
!SGMN=SGMN+SUM
!supply - demand
sol_RNMN(k,j)=SUM-DMDN
! UPDATE
IF(sol_RNMN(k,j)>0.)THEN
sol_NH3(k,j)=sol_NH3(k,j)+sol_RNMN(k,j)
! WNO3(ISL)=WNO3(ISL)+RNMN(ISL)
GO TO 21
END IF
X1=sol_NO3(k,j)+sol_RNMN(k,j)
IF(X1<0.)THEN
sol_RNMN(k,j)=-sol_NO3(k,j)
sol_NO3(k,j)=1.E-10
ELSE
sol_NO3(k,j)=X1
END IF
21 DF1=LSNTA
DF2=LMNTA
!!DF represents Demand from
!SNMN=SNMN+sol_RNMN(k,j)
!calculate P flows
!! compute humus mineralization on active organic p
hmp = 0.
hmp_rate = 0.
hmp_rate = 1.4* (HSNTA + HPNTA)/(sol_HSN(k,j) + sol_HPN(k,j) + 1.e-6)
!hmp_rate = 1.4* (HSNTA )/(sol_HSN(k,j) + sol_HPN(k,j) + 1.e-6)
hmp = hmp_rate*sol_orgp(k,j)
hmp = Min(hmp, sol_orgp(k,j))
sol_orgp(k,j) = sol_orgp(k,j) - hmp
sol_solp(k,j) = sol_solp(k,j) + hmp
!! compute residue decomp and mineralization of
!! fresh organic n and p (upper two layers only)
rmp = 0.
decr = 0.
decr = (LSCTA + LMCTA)/(sol_LSC(k,j) + sol_LMC(k,j) + 1.e-6)
decr = min(1., decr)
rmp = decr * sol_fop(k,j)
sol_fop(k,j) = sol_fop(k,j) - rmp
sol_solp(k,j) = sol_solp(k,j) + .8 * rmp
sol_orgp(k,j) = sol_orgp(k,j) + .2 * rmp
!calculate P flows
!SMS(9,ISL)=SMS(9,ISL)+RNMN(ISL)
LSCTA = Min(sol_LSC(k,j),LSCTA)
sol_LSC(k,j)=MAX(1.E-10,sol_LSC(k,j)-LSCTA)
LSLCTA = min(sol_LSLC(k,j),LSLCTA)
sol_LSLC(k,j)=MAX(1.E-10,sol_LSLC(k,j)-LSLCTA)
sol_LSLNC(k,j)=MAX(1.E-10,sol_LSLNC(k,j)-LSLNCTA)
LMCTA=MIN(sol_LMC(k,j),LMCTA)
IF (sol_LM(k,j) > 0.) THEN
RTO = MAX(0.42,sol_LMC(k,j)/sol_LM(k,j))
sol_LM(k,j) = sol_LM(k,j) - LMCTA/RTO
sol_LMC(k,j) = sol_LMC(k,j) - LMCTA
END IF
!sol_LMC(k,j)=MAX(1.E-10,sol_LMC(k,j)-LMCTA)
!sol_LM(k,j)=MAX(1.E-10,sol_LM(k,j)-LMCTA/.42)
sol_LSL(k,j)=MAX(1.E-10,sol_LSL(k,j)-LSLCTA/.42)
sol_LS(k,j)=MAX(1.E-10,sol_LS(k,j)-LSCTA/.42)
X3=APX*HPCTA+ASX*HSCTA+A1*(LMCTA+LSLNCTA)
sol_BMC(k,j)=sol_BMC(k,j)-BMCTA+X3
!DeltaBMC = DeltaBMC -BMCTA+X3
DF3=BMNTA-NCBM*X3
!!DF3 is the supply of BMNTA - demand of N to meet the Passive, Slow, Metabolic, and Non-lignin Structural
!! C pools transformaitons into microbiomass pool
X1=.7*LSLCTA+BMCTA*(1.-ABP-ABCO2)
sol_HSC(k,j)=sol_HSC(k,j)-HSCTA+X1
DF4=HSNTA-NCHS*X1
!!DF4 Slow pool supply of N - N demand for microbiomass C transformed into slow pool
X1=HSCTA*ASP+BMCTA*ABP
sol_HPC(k,j)=sol_HPC(k,j)-HPCTA+X1
DF5=HPNTA-NCHP*X1
!!DF5 Passive pool demand of N - N demand for microbiomass C transformed into passive pool
DF6=sol_min_n-sol_NO3(k,j)-sol_NH3(k,j)
!!DF6 Supply of mineral N - available mineral N = N demanded from mineral pool
!SMS(10,ISL)=SMS(10,ISL)-DF6
ADD=DF1+DF2+DF3+DF4+DF5+DF6
ADF1=ABS(DF1)
ADF2=ABS(DF2)
ADF3=ABS(DF3)
ADF4=ABS(DF4)
ADF5=ABS(DF5)
TOT=ADF1+ADF2+ADF3+ADF4+ADF5
XX=ADD/(TOT+1.E-10)
sol_LSN(k,j)=MAX(.001,sol_LSN(k,j)-DF1+XX*ADF1)
sol_LMN(k,j)=MAX(.001,sol_LMN(k,j)-DF2+XX*ADF2)
sol_BMN(k,j)=sol_BMN(k,j)-DF3+XX*ADF3
sol_HSN(k,j)=sol_HSN(k,j)-DF4+XX*ADF4
sol_HPN(k,j)=sol_HPN(k,j)-DF5+XX*ADF5
sol_RSPC(k,j)=.3*LSLCTA+A1CO2*(LSLNCTA+LMCTA)+ABCO2*BMCTA+ASCO2*HSCTA+APCO2*HPCTA
rspc_d(j) = rspc_d(j) + sol_RSPC(k,j)
!SMM(74,MO)=SMM(74,MO)+RSPC(ISL)
!SMS(8,ISL)=SMS(8,ISL)+RSPC(ISL)
!TRSP=TRSP+RSPC(ISL)
!VAR(74)=VAR(74)+RSPC(ISL)
!RSD(ISL)=.001*(WLS(ISL)+WLM(ISL))
sol_rsd(k,j)= sol_LS(k,j)+sol_LM(k,j)
sol_orgn(k,j) = sol_HPN(k,j)
sol_aorgn(k,j) = sol_HSN(k,j)
sol_fon(k,j) = sol_LMN(k,j) + sol_LSN(k,j)
sol_cbn(k,j) = 100*(sol_LSC(k,j)+sol_LMC(k,j) +sol_HSC(k,j) + sol_HPC(k,j) + sol_BMC(k,j))/sol_mass
!! septic changes 1/28/09 gsm
!! compute denitrification while simulating septic tank
!wdn = 0.
!if (i_sep(j) /= k .and. ipop_sep(j) > 0) then
!! compute soil water factor
! sut = 0.
!! change for domain error 1/29/09 gsm check with Jeff !!!
! if (sol_st(kk,j) < 0.) sol_st(kk,j) = .0000001
! sut = .1 + .9 * Sqrt(sol_st(kk,j) / sol_fc(kk,j))
! sut = Min(1., sut)
! sut = Max(.05, sut)
! !!compute soil temperature factor
! xx = 0.
! cdg = 0.
! xx = sol_tmp(kk,j)
! cdg = .9 * xx / (xx + Exp(9.93 - .312 * xx)) + .1
! cdg = Max(.1, cdg)
! if (sut >= sdnco) then
! wdn = sol_no3(k,j) * (1. - Exp(-cdn * cdg * sol_cbn(k,j)))
! else
! wdn = 0.
! endif
! sol_no3(k,j) = sol_no3(k,j) - wdn
!end if
! septic changes 1/28/09 gsm
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! hmntl |kg N/ha |amount of nitrogen moving from active
!! |organic to nitrate pool in soil profile
!! |on current day in HRU
!! hmptl |kg P/ha |amount of phosphorus moving from the
!! |organic to labile pool in soil profile
!! |on current day in HRU
!! rmn2tl |kg N/ha |amount of nitrogen moving from the fresh
!! |organic (residue) to the nitrate(80%) and
!! |active organic(20%) pools in soil profile
!! |on current day in HRU
!! rmptl |kg P/ha |amount of phosphorus moving from the
!! |fresh organic (residue) to the labile(80%)
!! |and organic(20%) pools in soil profile
!! |on current day in HRU
!! rwntl |kg N/ha |amount of nitrogen moving from active
!! |organic to stable organic pool in soil
!! |profile on current day in HRU
!! sol_aorgn(:,:)|kg N/ha |amount of nitrogen stored in the active
!! |organic (humic) nitrogen pool
!! sol_fon(:,:) |kg N/ha |amount of nitrogen stored in the fresh
!! |organic (residue) pool
!! sol_fop(:,:) |kg P/ha |amount of phosphorus stored in the fresh
!! |organic (residue) pool
!! sol_no3(:,:) |kg N/ha |amount of nitrogen stored in the
!! |nitrate pool in soil layer
!! sol_orgn(:,:) |kg N/ha |amount of nitrogen stored in the stable
!! |organic N pool
!! sol_orgp(:,:) |kg P/ha |amount of phosphorus stored in the organic
!! |P pool in soil layer
!! sol_rsd(:,:) |kg/ha |amount of organic matter in the soil
!! |classified as residue
!! sol_solp(:,:) |kg P/ha |amount of phosohorus stored in solution
!! wdntl |kg N/ha |amount of nitrogen lost from nitrate pool
!! |by denitrification in soil profile on
!! |current day in HRU
!! wshd_dnit |kg N/ha |average annual amount of nitrogen lost from
!! |nitrate pool due to denitrification in
!! |watershed
!! wshd_hmn |kg N/ha |average annual amount of nitrogen moving
!! |from active organic to nitrate pool in
!! |watershed
!! wshd_hmp |kg P/ha |average annual amount of phosphorus moving
!! |from organic to labile pool in watershed
!! wshd_rmn |kg N/ha |average annual amount of nitrogen moving
!! |from fresh organic (residue) to nitrate
!! |and active organic pools in watershed
!! wshd_rmp |kg P/ha |average annual amount of phosphorus moving
!! |from fresh organic (residue) to labile
!! |and organic pools in watershed
!! wshd_rwn |kg N/ha |average annual amount of nitrogen moving
!! |from active organic to stable organic pool
!! |in watershed
! call ndenit(k,j,cdg,wdn,0.05)
!! end if
!! summary calculations
!! calculations are based on century model, and not alighned with SWAT old algorithm yet.
if (curyr > nyskip) then
hmn = 0.
hmn = sol_RNMN(k,j)
wshd_hmn = wshd_hmn + hmn * hru_dafr(j)
rwn = 0.
rwn = HSNTA
wshd_rwn = wshd_rwn + rwn * hru_dafr(j)
wshd_hmp = wshd_hmp + hmp * hru_dafr(j)
rmn1 = 0.
rmn1 = (LSNTA+LMNTA)
wshd_rmn = wshd_rmn + rmn1 * hru_dafr(j)
wshd_rmp = wshd_rmp + rmp * hru_dafr(j)
wshd_dnit = wshd_dnit + wdn * hru_dafr(j)
hmntl = hmntl + hmn
rwntl = rwntl + rwn
hmptl = hmptl + hmp
rmn2tl = rmn2tl + rmn1
rmptl = rmptl + rmp
wdntl = wdntl + wdn
end if
end if
end do
!write (*,*) iyr,i, DeltaBMC
!write (*,*) iyr,i, wdntl
return
end