-
Notifications
You must be signed in to change notification settings - Fork 0
/
watbal.f90
2107 lines (1683 loc) · 74.9 KB
/
watbal.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
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!**********************************************************************
! WATBAL.FOR
! RAD November 2008.
!
! This file contains all the subroutines for calculating the water
! balance. Based on SPA water balance routines, with various modifications.
!
! Modifications mostly apply to soil water potential, hydraulic and thermal
! conductivity functions, and format of input parameters. The code is also organized
! very differently into subroutines and functions than the original SPA code.
!
! The subroutines are:
!
! INITWATBAL - initializes water balance variables.
! CALCSOILPARS - calculates relative water uptake from soil layers
! WATBALLAY - calls water balance routines, for layered soil.
! HEATBALANCE - redistributes heat, and calls the soil T profile calculation.
! FINDSOILTK - estimates soil surface temperature from heat balance.
! ENERGYFUN - returns the soil surface energy balance closure given soil T.
! ENERGYCALC - calculates the components of the soil surface energy balance.
! SOILWPFUN - soil water potential function.
! SOILCONDFUN - soil hydraulic conductivity function.
! SOILRESCALC - soil to root hydraulic conductance function.
! WATERUPTAKELAYER - called by CALCSOILPARS.
! CANOPY_BALANCE - canopy interception, canopy drainage and wet ET.
! CANSTOR - canopy storage function, integrated by CANOPY_BALANCE
! WETTINGLAYERS - surface wetting and drying, calculates dry layer thickness.
! QEFLUX - soil evaporation (latent heat flux, QE).
! SOIL_BALANCE - integrates soil gravitational drainage.
! SOILSTOR - soil storage function, integrated by SOIL_BALANCE.
! INFILTRATE - infiltration of rainfall into top layers of soil.
! WATERTHERMAL - redistributes heat according to water movement.
! CRANKNICHOLS - calculates soil T profile.
! THERMCONDSUB - soil thermal conductivity.
! VOLHCFUN - Calculate volumetric heat capacity.
! SCALEUP - scales single-tree estimates of ET and radiation absorbed
! to stand-average values.
! ASSIGNSOILWATER - Assign soil water content, depending on parameters
! IWATBALSIM, WSOILMETHOD and USEMEASSW
!
!**********************************************************************
!**********************************************************************
SUBROUTINE INITWATBAL(LAYTHICK,WETTINGBOT,WETTINGTOP, &
POREFRAC,WATERGAIN,WATERLOSS,PPTGAIN, &
INITWATER,DRYTHICKMIN,DRYTHICK, &
CANOPY_STORE, SURFACE_WATERMM, &
FRACWATER,WSOIL,WSOILROOT, &
NLAYER,NROOTLAYER,ICEPROP, &
QE,RUNOFF,OUTFLOW,SOILDEPTH, &
SOILDATA,USEMEASSW)
! Initializes various water balance related variables.
! RAD, May 2008
!**********************************************************************
USE maestcom
USE switches
IMPLICIT NONE
INTEGER SOILDATA,USEMEASSW,IOERROR
INTEGER NLAYER,I,NROOTLAYER,IZEROVAR
REAL LAYTHICK(MAXSOILLAY),WETTINGBOT(10),WETTINGTOP(10)
REAL POREFRAC(MAXSOILLAY),WATERGAIN(MAXSOILLAY)
REAL WATERLOSS(MAXSOILLAY),INITWATER(MAXSOILLAY)
REAL FRACWATER(MAXSOILLAY),PPTGAIN(MAXSOILLAY)
REAL ICEPROP(MAXSOILLAY)
REAL DRYTHICK,DRYTHICKMIN,QE,RUNOFF,OUTFLOW
REAL CANOPY_STORE,SURFACE_WATERMM,SOILDEPTH
REAL WSOIL,WSOILROOT,ZEROVAR
FRACWATER = INITWATER
ICEPROP = 0.
DRYTHICK = DRYTHICKMIN
WETTINGTOP = 0.
WETTINGBOT = 0.
WETTINGBOT(1) = LAYTHICK(1) ! See io.for in SPA
QE = 0. ! Soil evaporation (as latent heat flux)
WATERGAIN = 0. ! array of water gains per layer (m)
WATERLOSS = 0. ! array of water losses per layer (m)
PPTGAIN = 0.
RUNOFF = 0. ! Cumulative Runoff
OUTFLOW = 0.
CANOPY_STORE = 0.
SURFACE_WATERMM = 0.
! Soil properties in the layer below the 'lowest layer'. This is the 'soil core'
! that water drains into. In SPA, the 'core' counter is used for this layer.
POREFRAC(NLAYER+1) = POREFRAC(NLAYER)
FRACWATER(NLAYER+1) = FRACWATER(NLAYER)
LAYTHICK(NLAYER+1) = LAYTHICK(NLAYER)
! Total depth of rooted soil (m), for conversion of WSOIL to water content.
SOILDEPTH = SUM(LAYTHICK(1:NROOTLAYER))
! Initial total soil water storage in the whole soil, and in the rooted profile.
WSOIL = 0.
WSOILROOT = 0.
DO I = 1,NLAYER
WSOIL = WSOIL + INITWATER(I)*LAYTHICK(I)*1000
END DO
DO I = 1,NROOTLAYER
WSOILROOT = WSOILROOT + INITWATER(I)*LAYTHICK(I)*1000
END DO
! If there is soil data in the met.dat file, but the user has set usemeassw=0
! as well as iwatbalsim=0, set usemeassw to 1, and print a warning.
IF(SOILDATA.GT.0.AND.USEMEASSW.EQ.0)THEN
USEMEASSW = 1
CALL SUBERROR('WARNING: USEMEASSW is set to 1, measured soil water is used.',IWARN,IOERROR)
ENDIF
! Write initial water storage and fraction to output file, all others set to 0
! Move this to inout.for, use same routine as for output somehow...
ZEROVAR = -999.00
IZEROVAR = 0
IF (IOFORMAT .EQ. 1) THEN
WRITE (UWATBAL) REAL(IZEROVAR),REAL(IZEROVAR),WSOIL,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR
ELSE IF (IOFORMAT .EQ. 0) THEN
WRITE (UWATBAL,520)IZEROVAR,IZEROVAR,WSOIL,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR, &
ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR,ZEROVAR
520 FORMAT (I7,I7,51(F14.4,1X))
END IF
RETURN
END SUBROUTINE INITWATBAL
!**********************************************************************
SUBROUTINE CALCSOILPARS(NLAYER,NROOTLAYER,ISPEC,SOILWP,FRACWATER, &
FRACORGANIC,POREFRAC,SOILCOND,THERMCOND, &
ROOTMASS,ROOTLEN,LAYTHICK,ICEPROP, &
EQUALUPTAKE,RETFUNCTION, &
USEMEASSW, SOILDATA, SOILMOISTURE, &
PSIE,BPAR,KSAT,ROOTRESIST, &
ROOTRESFRAC, &
ROOTRAD,MINROOTWP,TOTLAI, &
WIND, ZHT, Z0HT, GAMSOIL, &
WEIGHTEDSWP,TOTESTEVAP, &
FRACUPTAKE,TOTSOILRES,ALPHARET,WS,WR,NRET)
! Calculates soil water potential, hydraulic conductivity,
! soil-to-root hydraulic conductance (leaf specific), thermal
! conductivity, and fraction water uptake from each soil layer
! (that has roots). No uptake takes place yet (see WATBALLAY subroutine).
! RAD, May 2008
!**********************************************************************
USE maestcom
USE metcom
IMPLICIT NONE
INTEGER I, NLAYER, NROOTLAYER, EQUALUPTAKE, RETFUNCTION, SOILDATA
INTEGER USEMEASSW,ISPEC
REAL BPAR(MAXSOILLAY),PSIE(MAXSOILLAY),KSAT(MAXSOILLAY)
REAL ALPHARET(MAXSOILLAY),WS(MAXSOILLAY),WR(MAXSOILLAY),NRET(MAXSOILLAY)
REAL SOILWP(MAXSOILLAY),FRACWATER(MAXSOILLAY),POREFRAC(MAXSOILLAY)
REAL SOILCOND(MAXSOILLAY),ROOTMASS(MAXSOILLAY),ROOTLEN(MAXSOILLAY)
REAL LAYTHICK(MAXSOILLAY), ICEPROP(MAXSOILLAY)
REAl SOILRRES1(MAXSOILLAY), SOILRRES2(MAXSOILLAY)
REAL FRACUPTAKE(MAXSOILLAY),SOILRRES(MAXSOILLAY)
REAL THERMCOND(MAXSOILLAY), FRACORGANIC(MAXSOILLAY)
REAL MINROOTWP,KTOT,SOILZPD,GAMSOIL,TOTESTEVAP
REAL SOILMOISTURE,ROOTRESIST,ROOTRESFRAC,ROOTRAD,TOTLAI
REAL WIND,ZHT,Z0HT,WEIGHTEDSWP,TOTSOILRES
REAL, EXTERNAL :: SOILCONDFUN
REAL, EXTERNAL :: THERMCONDFUN
REAL, EXTERNAL :: SOILWPFUN
REAL, EXTERNAL :: SOILCONDFUN2
REAL, EXTERNAL :: GBCANMS
! Update soil water potential,soil conductivity, soil thermal conductivity and
! soil to root conductance for each soil layer.
DO I=1,NLAYER
SOILWP(I) = SOILWPFUN(FRACWATER(I),PSIE(I),BPAR(I), &
POREFRAC(I),ALPHARET(I),WS(I),WR(I),NRET(I),RETFUNCTION)
SOILCOND(I) = SOILCONDFUN(FRACWATER(I),KSAT(I), &
BPAR(I),POREFRAC(I),WS(I),WR(I),NRET(I),RETFUNCTION)
THERMCOND(I) = THERMCONDFUN(I, SOILWP(I), FRACWATER(I), &
POREFRAC(I), BPAR(I), &
FRACORGANIC(I),RETFUNCTION)
ENDDO
! Calculate soil conductivity and soil-to-root resistance if using measured soil water content.
IF(USEMEASSW.EQ.1)THEN
! Only supported (sofar) for soil water potential
IF(SOILDATA.EQ.POTENTIAL)THEN
SOILCOND(1) = SOILCONDFUN2(SOILMOISTURE,KSAT(1),BPAR(1), &
PSIE(1))
ENDIF
ENDIF
CALL SOILRESCALC(USEMEASSW, SOILCOND,ROOTRESIST, &
ROOTMASS,ROOTLEN,LAYTHICK,ROOTRAD, &
NROOTLAYER,SOILRRES1,SOILRRES2)
! Fractional water uptake from each soil layer
! Note that water is not actually taken up yet.
CALL WATERUPTAKELAYER(SOILWP,SOILRRES1,SOILRRES2, &
ROOTRESFRAC, &
MINROOTWP,TOTLAI, &
ICEPROP,EQUALUPTAKE, &
USEMEASSW, SOILDATA, SOILMOISTURE, &
ROOTLEN,NROOTLAYER,WEIGHTEDSWP, &
FRACUPTAKE,TOTSOILRES,LAYTHICK, &
TOTESTEVAP)
! Aerodynamic conductance between soil surface and air,
! assuming turbulent transfer (so that conductance is the same for
! momentum, heat and mass transfer (Jones 1992)).
SOILZPD = 0.0 ! Reference height is the soil surface.
GAMSOIL = GBCANMS(WIND,ZHT,Z0HT,SOILZPD)
RETURN
END
!**********************************************************************
SUBROUTINE WATBALLAY(IDAY,IHOUR,PPT,RUTTERB,RUTTERD,MAXSTORAGE, &
THROUGHFALL,RADINTERC, &
CANOPY_STORE, EVAPSTORE, DRAINSTORE, &
SURFACE_WATERMM, &
POREFRAC,WETTINGBOT,WETTINGTOP,NLAYER, &
NROOTLAYER,LAYTHICK,SOILTK,QE, &
TAIRK,VPDPA,WIND, &
ZHT,Z0HT,ZPD,PRESS,ETMM,ETMMSPEC,NOSPEC, &
USEMEASET,ETMEAS,FRACUPTAKESPEC, &
ICEPROP,FRACWATER,DRAINLIMIT, &
KSAT,BPAR,WSOIL,WSOILROOT,DISCHARGE, &
DRYTHICKMIN,DRYTHICK,QEMM,OVERFLOW, &
WATERGAIN,WATERLOSS,PPTGAIN,KEEPWET, &
EXPINF,WS,WR,NRET,RETFUNCTION)
! Do water balance for layered soil.
! Replaces subroutines WATERFLUXES and WATERTHERMAL in SPA
! No calculation of heat balance yet, but should be placed in separate routine.
!
! RAD, July 2008
!**********************************************************************
USE maestcom
IMPLICIT NONE
INTEGER I,J,RR,NLAYER,NROOTLAYER,USEMEASET,IDAY,IHOUR
INTEGER KEEPWET
INTEGER RETFUNCTION,NOSPEC,ISPEC
REAL POREFRAC(MAXSOILLAY),WETTINGBOT(10),WETTINGTOP(10)
REAL LAYTHICK(MAXSOILLAY)
REAL WATERGAIN(MAXSOILLAY),WATERLOSS(MAXSOILLAY)
REAL FRACUPTAKESPEC(MAXSOILLAY,MAXSP),FRACWATER(MAXSOILLAY)
REAL FRACUPTAKE(MAXSOILLAY)
REAL ICEPROP(MAXSOILLAY),PPTGAIN(MAXSOILLAY)
REAL BPAR(MAXSOILLAY), KSAT(MAXSOILLAY)
REAL DRAINLIMIT(MAXSOILLAY)
real check1, check2, check3, check4
REAL MAXSTORAGE,LAMBDASOIL,ETMMSPEC(MAXSP)
REAL SOILTK,TAIRK,VPDPA,VPDKPA,QEM,QEMM,SURFACE_WATERMM
REAL PPT,EVAPSTORE,DRAINSTORE,WSOILROOT,TOTESTEVAPMM
REAL SOILTC,TAIRC,WSOIL,SNOW,QE,THROUGHFALL,WIND,ZHT,Z0HT
REAL ZPD,PRESS,RADINTERC,RUTTERB,RUTTERD,CANOPY_STORE,DRYTHICKMIN
REAL DRYTHICK,ETLOSS,ETMEAS,ETMM,DISCHARGE,EXPINF,OVERFLOW
REAL WATERCONTENT,ICECONTENT
REAL, EXTERNAL :: HEATEVAP
REAL WS(MAXSOILLAY),WR(MAXSOILLAY),NRET(MAXSOILLAY)
! Conversions
SOILTC = SOILTK - FREEZE
TAIRC = TAIRK - FREEZE
VPDKPA = VPDPA / 1000
! Latent heat of soil evaporation.
LAMBDASOIL = HEATEVAP(SOILTC)
! Init.
WSOIL = 0.
WATERLOSS = 0.
WATERGAIN = 0.
SNOW = 0.
! mm t-1 (mm per timestep)
QEMM = (-QE/LAMBDASOIL)*SPERHR
! m t-1
QEM = QEMM / 1000
! Assume that water infiltrates in one timestep.
SURFACE_WATERMM = 0.
! Determine throughfall.
IF(THROUGHFALL.EQ.1.0.OR.MAXSTORAGE.EQ.0.0)THEN
SURFACE_WATERMM = PPT
EVAPSTORE = 0.0
DRAINSTORE = 0.0
ELSE
CALL CANOPY_BALANCE(PPT,WIND,ZHT,Z0HT,ZPD, &
PRESS,TAIRC,RADINTERC, &
VPDPA,THROUGHFALL, &
RUTTERB,RUTTERD,MAXSTORAGE, &
CANOPY_STORE, SURFACE_WATERMM, &
EVAPSTORE, DRAINSTORE)
ENDIF
! Calculates the thickness of the top dry layer (if any).
CALL WETTINGLAYERS(POREFRAC,WETTINGBOT,WETTINGTOP, &
SURFACE_WATERMM,SNOW,TAIRK,QE, &
NLAYER,LAYTHICK,DRYTHICKMIN,DRYTHICK)
! Option to keep soil wet, also keep DRYTHICK = 0
IF(KEEPWET.EQ.1)DRYTHICK = 0.001
! From which layer is soil evaporation withdrawn?
! Note that it currently can only come from 1st or 2nd layer!
! This probably puts constraints on the thickness of soil layers (default = 10cm)
IF(DRYTHICK.LT.LAYTHICK(1))THEN
RR=1 !the dry zone does not extend beneath the top layer
ELSE
RR=2 !it does
ENDIF
! Evaporation
IF(QE.LT.0.)THEN
WATERLOSS(RR) = WATERLOSS(RR) + QEM
! Dew formation (positive values of QE are now actually prevented in QEFLUX, see there).
ELSE
WATERGAIN(1) = WATERGAIN(1) - QEM
ENDIF
! Use measured or modelled ET for water balance calculations:
IF(USEMEASET.EQ.1)THEN
ETLOSS = MAX(0.0, ETMEAS / 1000)
DO I=1,NROOTLAYER
WATERLOSS(I) = WATERLOSS(I) + ETLOSS*FRACUPTAKE(I)
ENDDO
! Use modelled; allow for multiple species.
ELSE
DO ISPEC=1,NOSPEC
! Water loss from each rooted layer (i.e. *root water uptake)
ETLOSS = ETMMSPEC(ISPEC) / 1000
FRACUPTAKE(1:MAXSOILLAY) = FRACUPTAKESPEC(1:MAXSOILLAY, ISPEC)
DO I=1,NROOTLAYER
WATERLOSS(I) = WATERLOSS(I) + ETLOSS*FRACUPTAKE(I)
ENDDO
ENDDO
ENDIF
check1 = sum(fracuptake(1:nrootlayer))
check2 = etmm/1000 + qem
check3 = sum(etmmspec(1:nospec))/1000 + qem
check4 = sum(waterloss(1:nrootlayer))
! Calculate drainage for each soil layer
DO J = 1,NLAYER
CALL SOIL_BALANCE(J, POREFRAC, ICEPROP, FRACWATER, &
LAYTHICK, DRAINLIMIT(J), WATERLOSS, WATERGAIN, &
KSAT(J), BPAR(J),WS(J),WR(J),NRET(J),RETFUNCTION)
ENDDO
! Loss of water at lowest soil layer (discharge, or deep drainage) (mm).
DISCHARGE = 1000 * WATERGAIN(NLAYER + 1)
! Infiltration of water reaching the surface
! NOTE: WATERGAIN and WATERLOSS are not updated here.
! Instead, PPTGAIN is an output array that is added to water balance.
! This is because of heat balance: temperature of rain equals air temperature,
! temperature of soil water equals soil temperature.
CALL INFILTRATE(SURFACE_WATERMM,NLAYER,POREFRAC,FRACWATER, &
LAYTHICK,WATERGAIN,WATERLOSS, &
EXPINF,PPTGAIN,OVERFLOW)
! Add and subtract gains and losses for each layer
DO J = 1,NLAYER
! Option to not change initial soil water.
IF(KEEPWET.NE.1)THEN
! M of water
WATERCONTENT = (FRACWATER(J)*(1. - ICEPROP(J)))*LAYTHICK(J)
ICECONTENT = (FRACWATER(J)*ICEPROP(J))*LAYTHICK(J)
WATERCONTENT = MAX(0., WATERCONTENT + WATERGAIN(J) + &
PPTGAIN(J) - WATERLOSS(J))
! Volumetric water content
FRACWATER(J) = (WATERCONTENT+ICECONTENT) / LAYTHICK(J)
IF(WATERCONTENT.EQ.0.0)THEN
ICEPROP(J) = 0.0
ELSE
ICEPROP(J) = ICECONTENT / (WATERCONTENT+ICECONTENT)
ENDIF
ENDIF
! Total soil water storage (mm)
WSOIL = WSOIL + FRACWATER(J)*LAYTHICK(J)*1000
! Do error checking here...
ENDDO
! Total soil water storage in rooted zone
WSOILROOT = 0.
DO J = 1,NROOTLAYER
WSOILROOT = WSOILROOT + FRACWATER(J)*LAYTHICK(J)*1000
ENDDO
RETURN
END
!**********************************************************************
SUBROUTINE HEATBALANCE(NLAYER,FRACWATER,POREFRAC,TAIRK,SOILTK, &
SOILTEMP,LAYTHICK,WATERGAIN,WATERLOSS, &
PPTGAIN,THERMCOND)
! Does the heat balance calculation: soil surface temperature,
! which is based on energy balance closure (net radiation,
! sensible and latent heat flux, soil heat flux),
! soil T profile, redistribution of heat due to water movement.
! Based on various bits of code from SPA, in subroutines
! SOILDAY, WATERTHERMAL, ENERGY, ETC.
!**********************************************************************
USE maestcom
IMPLICIT NONE
INTEGER NLAYER
REAL THERMCOND(MAXSOILLAY),SOILTEMP(MAXSOILLAY)
REAL FRACWATER(MAXSOILLAY),POREFRAC(MAXSOILLAY)
REAL LAYTHICK(MAXSOILLAY),WATERGAIN(MAXSOILLAY)
REAL WATERLOSS(MAXSOILLAY),PPTGAIN(MAXSOILLAY)
REAL VOLHC(MAXSOILLAY)
REAL TAIRK,SOILTK
! Update soil temperature based on movement (and uptake) of soil water.
! This should be a small component of the heat balance.
CALL WATERTHERMAL(NLAYER, FRACWATER, POREFRAC, &
SOILTEMP, LAYTHICK, WATERGAIN, &
TAIRK, WATERLOSS, PPTGAIN, VOLHC)
! Calculate soil temperature profile with Crank-Nicholson scheme.
CALL CRANKNICHOLS(NLAYER, LAYTHICK, SOILTK, SOILTEMP, &
VOLHC, THERMCOND)
RETURN
END
!**********************************************************************
SUBROUTINE FINDSOILTK(IDAY,TAIRK, GAMSOIL, &
PRESSPA, SOILTK, SOILTK2, VPDKPA, RGLOB, &
THERMCOND1, LAYTHICK1, POREFRAC1, &
SOILWP1,DRYTHICK,TORTPAR,VIEWFACTOR)
! Finds the soil surface temperature from energy balance, by
! calling the ZBRENT routine, which uses the ENERGY function.
! Based on SPA (bits in SOILDAY mostly). (June 2008, RAD).
!**********************************************************************
USE maestcom
IMPLICIT NONE
INTEGER IDAY
INTEGER EXTRAINT(10)
REAL EXTRAPARS(EXTRAPARDIM)
REAL LAYTHICK1,VIEWFACTOR
REAL TAIRK,GAMSOIL,PRESSPA,SOILTK,SOILTK2,VPDKPA
REAL RGLOB,THERMCOND1,POREFRAC1,SOILWP1
REAL DRYTHICK,TORTPAR
REAL T1,T2,XACC,TEST
REAL, EXTERNAL :: ENERGYFUN
REAL, EXTERNAL :: ZBRENT
! Put parameters into EXTRAPARS array:
EXTRAPARS(1) = GAMSOIL
EXTRAPARS(2) = PRESSPA
EXTRAPARS(3) = SOILTK2
EXTRAPARS(4) = VPDKPA
EXTRAPARS(5) = RGLOB
EXTRAPARS(6) = TAIRK
EXTRAPARS(7) = THERMCOND1
EXTRAPARS(8) = LAYTHICK1
EXTRAPARS(9) = POREFRAC1
EXTRAPARS(10) = SOILWP1
EXTRAPARS(11) = DRYTHICK
EXTRAPARS(12) = TORTPAR
EXTRAPARS(13) = VIEWFACTOR
! Set bounds for root-finding (quite liberal bounds!).
T1 = TAIRK - 50.
T2 = TAIRK + 50.
! Error tolerance for soil temperature numerical solution (degrees C).
XACC = 0.0001
! Find the soil surface temperature that gives closure in energy balance:
SOILTK = ZBRENT(ENERGYFUN,T1,T2,XACC,EXTRAPARS, EXTRAINT)
! Print warning if the solution is at the bounds.
IF(SOILTK.EQ.T1.OR.SOILTK.EQ.T2) &
CALL SUBERROR('WARNING: SOIL T SOLUTION AT BOUNDS.', &
IWARN,0)
RETURN
END
!**********************************************************************
REAL FUNCTION ENERGYFUN(SOILTK,EXTRAPARS,EXTRAINT)
! Calculates the energy balance closure (based on Hinzmann et al. 1998).
! Used to find the soil surface temperature that achieves closure
! (i.e. zero net balance),in the function FINDSOILTK.
! Taken from SPA (June 2008, RAD).
!**********************************************************************
USE maestcom
IMPLICIT NONE
REAL EXTRAPARS(EXTRAPARDIM)
INTEGER EXTRAINT(10)
REAL LAYTHICK1,VIEWFACTOR,GAMSOIL,PRESSPA,SOILTK2,VPDKPA
REAL RGLOB,TAIRK,THERMCOND1,POREFRAC1,SOILWP1,DRYTHICK
REAL TORTPAR,QE,QH,QN,QC,ESOIL,SOILTK
! Get parameters from EXTRAPARS array:
GAMSOIL = EXTRAPARS(1)
PRESSPA = EXTRAPARS(2)
SOILTK2 = EXTRAPARS(3)
VPDKPA = EXTRAPARS(4)
RGLOB = EXTRAPARS(5)
TAIRK = EXTRAPARS(6)
THERMCOND1 = EXTRAPARS(7)
LAYTHICK1 = EXTRAPARS(8)
POREFRAC1 = EXTRAPARS(9)
SOILWP1 = EXTRAPARS(10)
DRYTHICK = EXTRAPARS(11)
TORTPAR = EXTRAPARS(12)
VIEWFACTOR = EXTRAPARS(13)
! Subroutine that actually does all the work.
CALL ENERGYCALC(SOILTK,GAMSOIL,PRESSPA,SOILTK2, &
VPDKPA,RGLOB,TAIRK,THERMCOND1,LAYTHICK1, &
POREFRAC1,SOILWP1,DRYTHICK,TORTPAR, &
VIEWFACTOR, &
QH,QE,QN,QC,ESOIL)
! Energy balance:
ENERGYFUN = QH + QE + QN + QC
RETURN
END
!**********************************************************************
SUBROUTINE ENERGYCALC(SOILTK,GAMSOIL,PRESSPA,SOILTK2, &
VPDKPA,RGLOB,TAIRK,THERMCOND1,LAYTHICK1, &
POREFRAC1,SOILWP1, &
DRYTHICK,TORTPAR, &
VIEWFACTOR,QH,QE,QN,QC,ESOIL)
! Calculate components of the soil energy balance.
! For arguments, subscript (1 or 2) refers to soil layer (1 = surface).
! RAD, June 2008.
!**********************************************************************
USE maestcom
IMPLICIT NONE
REAL LAYTHICK1,VIEWFACTOR,SOILTK,GAMSOIL,PRESSPA,SOILTK2
REAL VPDKPA,RGLOB,TAIRK,THERMCOND1,POREFRAC1,SOILWP1
REAL DRYTHICK,TORTPAR,QH,QE,QN,QC,ESOIL,TAIRC
REAL RHO,QEFLUX
REAL, EXTERNAL :: RHOFUN
REAL, EXTERNAL :: ESOILFUN
! Conversions
TAIRC = TAIRK - FREEZE
! Density of air (kg m-3)
RHO = RHOFUN(TAIRK)
! Sensible heat flux (W m-2)
QH = CPAIR * RHO * GAMSOIL * (TAIRK - SOILTK)
! Latent heat flux (W m-2)
QE = QEFLUX(SOILTK,TAIRK,VPDKPA,POREFRAC1,SOILWP1, &
GAMSOIL,PRESSPA,DRYTHICK,TORTPAR)
! No soil evap if surface is frozen
IF(SOILTK.LE.FREEZE)QE = 0.
! Net radiation - emitted longwave varies with surface temp.
ESOIL = ESOILFUN(SOILTK)
QN = (1-SOILALBEDO)*RGLOB - ESOIL
! Soil heat flux (flux of heat out of layer 1 into layer 2).
QC = -THERMCOND1 * (SOILTK - SOILTK2)/LAYTHICK1
! Was 0.5*LAYTHICK1; but actually has to transport to middle of next layer (otherwise we get a mismatch!)
! This actually should be the depth from the middle of layer 1 to the middle of layer 1, which
! only is LAYTHICK1 when 1 and 2 have the same thickness!!!
RETURN
END
!**********************************************************************
REAL FUNCTION SOILWPFUN(WATERCONTENT,PSIE,BPAR, &
POROSITY,ALPHARET,WS,WR,NRET,RETFUNCTION)
! Soil water potential (MPa), using the Campbell (1974) equation.
! PSIE (MPa),BPAR(-),POROSITY(m3 m-3),WATERCONTENT (m3 m-3)
! Output in MPa
! Or using van Genutchen (1980) equation, alpha (MPa-1),WS and WR (m3 m-3)
! RAD, April 2008.
!**********************************************************************
USE maestcom
IMPLICIT NONE
INTEGER RETFUNCTION,IOFATAL
REAL WATERCONTENT,PSIE,BPAR,POROSITY
REAL ALPHARET,WS,WR,NRET
SOILWPFUN = 1.0 ! INIT WITH WRONG VALUE.
IF(WATERCONTENT.EQ.0.) THEN
SOILWPFUN = -999.0 ! bug mathias décembre 2012 pour éviter valeure -Infinity qui crée NA dans WEIGHTEDSMP
ELSE
IF(RETFUNCTION.EQ.1)THEN
SOILWPFUN = PSIE*(WATERCONTENT/POROSITY)**(-BPAR)
ENDIF
IF(RETFUNCTION.EQ.2)THEN
SOILWPFUN = PSIE*(WATERCONTENT)**(-BPAR)
ENDIF
IF(RETFUNCTION.EQ.3) THEN
IF(WATERCONTENT.LT.WR)THEN
SOILWPFUN = -999.0
ELSE IF (WATERCONTENT.GT.WS) THEN
SOILWPFUN = -0.0001
ELSE
SOILWPFUN = - 1/ALPHARET * (((WS-WR)/(WATERCONTENT-WR))**(NRET/(NRET-1)) - 1)**(1/NRET)
END IF
END IF
ENDIF
IF(SOILWPFUN.GT.0) THEN
SOILWPFUN = 0.0
CALL SUBERROR('ERROR IN SOIL WATER POTENTIAL CALCULATION', &
IOFATAL,0)
ENDIF
IF(SOILWPFUN.LT.-999) THEN ! modification mathias decembre 2012 pour eviter bug
SOILWPFUN = -999.0
ENDIF
RETURN
END ! SOILWPFUN
!**********************************************************************
REAL FUNCTION SOILCONDFUN(WATERCONTENT,KSAT,BPAR,POROSITY,WS,WR,N,RETFUNCTION)
! Soil hydraulic conductivity, using the Campbell (1974) equation. or the van Genutchen equation
! Inputs:
! KSAT (ms-1),BPAR(-),POROSITY(m3 m-3),WATERCONTENT (m3 m-3)
! Output in m s-1
! RAD, April 2008
!**********************************************************************
USE maestcom
IMPLICIT NONE
INTEGER RETFUNCTION
REAL KSAT,WATERCONTENT,BPAR,POROSITY
REAL WS,WR,N,ALPHAPOT
IF (RETFUNCTION.EQ.1.OR.RETFUNCTION.EQ.2) THEN
! Campbell 1974
SOILCONDFUN = KSAT*(WATERCONTENT/POROSITY)**(2*BPAR + 3)
ELSE IF (RETFUNCTION.EQ.3) THEN
! van genuchten
IF (WATERCONTENT.LT.WR) THEN
SOILCONDFUN = 0
ELSE IF (WATERCONTENT.GT.WS) THEN
SOILCONDFUN = KSAT
ELSE
ALPHAPOT = (( (WS - WR)/(WATERCONTENT-WR) )**(N/(N-1)) - 1) ** (1/N)
SOILCONDFUN = KSAT * (1 - ALPHAPOT**(N-1) * ( 1+ ALPHAPOT**N )**(-1+1/N) ) **2 &
/ (1 + ALPHAPOT**N)**(1/2 - 1/(2*N))
END IF
END IF
! Avoid underflow (avoiding very small numbers).
IF(SOILCONDFUN.LT.1E-30) SOILCONDFUN = 1E-30
END !SOILCONDFUN
!**********************************************************************
REAL FUNCTION SOILCONDFUN2(SOILWP,KSAT,BPAR,PSIE)
! Soil hydraulic conductivity, using the Campbell (1974) equation.
! This is an alternatve to SOILCONDFUN, which uses water content directly.
! Inputs:
! KSAT (ms-1),BPAR(-),PSIE(MPa),SOILWP(MPa)
! Output in m s-1
! RAD, April 2008
!**********************************************************************
USE maestcom
IMPLICIT NONE
REAL KSAT, SOILWP, PSIE, BPAR
! Campbell 1974
SOILCONDFUN2 = KSAT*(PSIE/SOILWP)**(2+3/BPAR)
! Avoid underflow (avoiding very small numbers).
IF(SOILCONDFUN2.LT.1E-30) SOILCONDFUN2 = 1E-30
END !SOILCONDFUN2
!**********************************************************************
SUBROUTINE SOILRESCALC(USEMEASSW,SOILCOND,ROOTRESIST,ROOTMASS, &
ROOTLEN,LAYTHICK,ROOTRAD,NROOTLAYER, &
SOILR1,SOILR2)
! Calculate Soil-to-root hydraulic resistance
! Output is SOILR in MPa s m2 mmol-1
! Modified from SPA (RAD, April 2008).
!**********************************************************************
USE maestcom
IMPLICIT NONE
INTEGER NROOTLAYER,USEMEASSW,NITER,I
REAL LSOIL,MPAM,LA
REAL SOILCOND(MAXSOILLAY),ROOTMASS(MAXSOILLAY),ROOTLEN(MAXSOILLAY)
REAL LAYTHICK(MAXSOILLAY), SOILRRES(MAXSOILLAY)
REAL SOILR1(MAXSOILLAY), SOILR2(MAXSOILLAY)
REAL ROOTRESIST,ROOTRAD,DEPTH
REAL ROOTRESCONS,RS,RS2,SPAROOTRESIST
REAL KSOIL,KS,LOGRR,TOTROOT,MEANROOTLEN
REAL, EXTERNAL :: SOILCONDFUN
! Hydraulic head for conversion (MPa / m)
MPAM = 0.009807
! New representation of root resistance following Hacke et al.
! Root resistance in a layer is proportional to 1/rootmass and proportional
! to root length (i.e. depth of layer)
LA = 0.0
DO I = 1,NROOTLAYER
DEPTH = SUM(LAYTHICK(1:I)) - LAYTHICK(I)/2
IF(ROOTLEN(I).GT.0.0)THEN
LA = LA + DEPTH/ROOTLEN(I)
ENDIF
ENDDO
IF(LA.GT.0.0)THEN
ROOTRESCONS = ROOTRESIST / LA
ELSE
ROOTRESCONS = 1E09 ! Arbitrarily large number
ENDIF
DO I=1,NROOTLAYER
! Depth to middle of layer (ca. root path length)
DEPTH = SUM(LAYTHICK(1:I)) - LAYTHICK(I)/2
! Soil hydraulic conductivity in m2 s-1 MPa-1
!LSOIL = SOILCOND(I)/MPAM !converts from ms-1 to m2 s-1 MPa-1
! ... in mol m-1 s-1 MPa-1. Note that original KSAT was given in same units.
KSOIL = SOILCOND(I) / (H2OVW * GRAV * 1E-03)
IF(KSOIL.LT.1E-35)THEN !prevent floating point error
SOILRRES(I) = 1e35
ELSE
! Reformulated to match Duursma et al. 2008.
IF(ROOTLEN(I).GT.0.0)THEN
! Radius of soil cylinder around root in single-root model
RS = SQRT(1./(ROOTLEN(I)*PI))
LOGRR = LOG(RS/ROOTRAD)
IF(LOGRR.LT.0)CALL SUBERROR( &
'Root radius larger than soil-root radius - fine root density too high!', &
IWARN,0)
KS = ROOTLEN(I)*LAYTHICK(I)*2.0*pi*KSOIL/LOGRR
SOILR1(I) = 1/KS
! As in SPA:
!RS2 = LOG(RS/ROOTRAD)/(2.0*PI*ROOTLEN(I)*LAYTHICK(I)*LSOIL)
! convert from MPa s m2 m-3 to MPa s m2 mmol-1
!SOILR1(I) = RS2*1E-6*18*0.001
! Note : this component is calculated but not used (see wateruptakelayer proc). More research needed!
! Second component of below ground resistance related to root hydraulics.
SOILR2(I) = ROOTRESCONS * DEPTH / ROOTLEN(I)
ELSE
SOILR2(I) = 0.0
SOILR1(I) = 0.0
ENDIF
SOILRRES(I) = SOILR1(I) + SOILR2(I)
ENDIF
ENDDO
! When using measured soil water, use water content of top layer,
! but all fine roots to calculate total soil conductance
IF(USEMEASSW.EQ.1)THEN
KSOIL = SOILCOND(1) / (H2OVW * GRAV * 1E-03)
! Total average fine root density
TOTROOT = 0
DO I=1,NROOTLAYER
TOTROOT = TOTROOT + ROOTLEN(I)*LAYTHICK(I)
ENDDO
MEANROOTLEN = TOTROOT / SUM(LAYTHICK(1:NROOTLAYER))
! Average soil cylinder around roots.
RS = SQRT(1./(MEANROOTLEN*PI))
LOGRR = LOG(RS/ROOTRAD)
! Total soil conductance
KS = TOTROOT*2.0*pi*KSOIL/LOGRR
SOILR1(1) = 1/KS
SOILR2(1) = ROOTRESCONS * (SUM(LAYTHICK(1:NROOTLAYER))/2) / MEANROOTLEN
SOILRRES(1) = SOILR1(1) + SOILR2(1)
ENDIF
RETURN
END !SOILRESCALC
!**********************************************************************
SUBROUTINE WATERUPTAKELAYER(SOILWP,SOILRRES1,SOILRRES2, &
ROOTRESFRAC, &
MINROOTWP,TOTLAI, &
ICEPROP,EQUALUPTAKE, &
USEMEASSW, SOILDATA, SOILMOISTURE, &
ROOTLEN,NROOTLAYER,WEIGHTEDSWP, &
FRACUPTAKE,TOTSOILRES,LAYTHICK, &
TOTESTEVAP) ! rajout laythick mathias décembre 2012
! Function for deciding from which layer water is withdrawn,
! and to calculate soil water potential weighted by root fraction.
! Taken from SPA, April 2008 (RAD), but heavily modified.
!**********************************************************************
USE maestcom
USE metcom
IMPLICIT NONE
INTEGER I,NROOTLAYER,EQUALUPTAKE,SOILDATA,USEMEASSW
REAL ESTEVAP(MAXSOILLAY),FRACUPTAKE(MAXSOILLAY),RSUM
REAL FRACUPTAKE2(MAXSOILLAY), LAYTHICK(MAXSOILLAY)
REAL ICEPROP(MAXSOILLAY), SOILWP(MAXSOILLAY)
REAL SOILRRES1(MAXSOILLAY),SOILRRES2(MAXSOILLAY)
REAL SOILRRES(MAXSOILLAY),RESTOT(MAXSOILLAY)
REAL ROOTLEN(MAXSOILLAY),SWPDIFF
REAL MINROOTWP,KTOT,KTOT2,TOTESTEVAP,TOTESTEVAPMM
REAL WEIGHTEDSWP,TOTSOILRES,SOILMOISTURE
REAL EMAXLEAF,TOTLAI,FRACSUM,ROOTRESFRAC,TOTESTEVAPWET
!Reset.
TOTESTEVAP = 0.
TOTESTEVAPWET = 0.
WEIGHTEDSWP = 0.
ESTEVAP = 0.
FRACUPTAKE = 0.
! Total soil conductance at the leaf level; including soil and root component.
! SOILRRES2 is commented out : don't use root-component of resistance (is part of plant resistance)
SOILRRES = SOILRRES1 !+ SOILRRES2
! Estimated max transpiration from gradient-gravity / soil resis
IF(EQUALUPTAKE.EQ.0)THEN
DO I=1,NROOTLAYER
! Estimated maximum uptake rate with the current soil watyer
! if aboveground resistance is zero.
IF(SOILRRES(I).GT.0.0)THEN
ESTEVAP(I)=(SOILWP(I)-MINROOTWP)/SOILRRES(I)
ELSE
ESTEVAP(I)=0.0 ! When no roots present
ENDIF
! No negative uptake.
ESTEVAP(I)=MAX(0.,ESTEVAP(I))
! Soil water potential weighted by layer Emax (from SPA)
WEIGHTEDSWP = WEIGHTEDSWP + SOILWP(I)*ESTEVAP(I)
! No uptake from frozen soils
IF(ICEPROP(I).GT.0.)ESTEVAP(I)=0.
ENDDO
! Option 1 (not currently used) : SPA method to figure out relative water uptake.
! Fraction uptake in each layer by Emax in each layer (SPA)
IF(SUM(ESTEVAP(1:NROOTLAYER)).EQ.0.0)THEN
FRACUPTAKE = 0.0
ELSE
DO I=1,NROOTLAYER
FRACUPTAKE(I) = ESTEVAP(I)/SUM(ESTEVAP(1:NROOTLAYER))
ENDDO
ENDIF
! sum of EMAX for each layer
TOTESTEVAP = SUM(ESTEVAP(1:NROOTLAYER)) ! modification mathias décembre 2012
! Soil water potential is weighted by ESTEVAP.
IF(TOTESTEVAP.GT.0.)THEN
WEIGHTEDSWP = WEIGHTEDSWP / TOTESTEVAP
ELSE
WEIGHTEDSWP = 0.0
DO I=1,NROOTLAYER
WEIGHTEDSWP = WEIGHTEDSWP + SOILWP(I) * LAYTHICK(I)
ENDDO
WEIGHTEDSWP = WEIGHTEDSWP / SUM(LAYTHICK(1:NROOTLAYER))
ENDIF
! Resistances are in parallel. This variable is also not used.
RSUM = 0.0
DO I=1,NROOTLAYER
IF(SOILRRES(I).GT.0.0)THEN
RSUM = RSUM + 1/SOILRRES(I)
ENDIF
ENDDO
TOTSOILRES = 1/RSUM
ENDIF
! Debugging option; equal uptake in all layers.
! Also used (for now) when USEMEASET=1.
IF(EQUALUPTAKE.EQ.1)THEN
FRACUPTAKE = 1.0 / REAL(NROOTLAYER)
TOTSOILRES = SOILRRES1(1)
ENDIF
! Option 2 (currently used) :
! Taylor and Keppler: relative water uptake is
! proportional to root length density and Psi difference.
! See : Taylor, H.M. and B. Keppler. 1975. Water uptake by cotton root systems:
! an examination of assumptions in the single root model. Soil Science. 120:57-67.
DO I=1,NROOTLAYER
IF(SUM(ESTEVAP(1:NROOTLAYER)).GT.0.)THEN
SWPDIFF = MAX(0., (SOILWP(I)-MINROOTWP))
FRACUPTAKE2(I) = ROOTLEN(I)*SWPDIFF
ELSE
! No water uptake possible.
FRACUPTAKE2(I) = 0.
ENDIF
ENDDO
IF(SUM(FRACUPTAKE2(1:NROOTLAYER)).GT.0)THEN
FRACUPTAKE2 = FRACUPTAKE2 / SUM(FRACUPTAKE2(1:NROOTLAYER)) ! Make sure that it sums to 1.
ELSE
FRACUPTAKE2 = 0.0
ENDIF
! If using measured soil water, replace with measurement.
IF(USEMEASSW.EQ.1)THEN
IF(SOILDATA.EQ.POTENTIAL)THEN
WEIGHTEDSWP = SOILMOISTURE