-
Notifications
You must be signed in to change notification settings - Fork 92
/
FatesPlantHydraulicsMod.F90
5344 lines (4559 loc) · 269 KB
/
FatesPlantHydraulicsMod.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
module FatesPlantHydraulicsMod
! ==============================================================================================
! This module contains the relevant code for plant hydraulics. Currently, one hydraulics module
! is available. Other methods of estimating plant hydraulics may become available in future
! releases. For now, please cite the following reference if this module is used to generate
! published research:
!
! Christoffersen, B.O., Gloor, M., Fauset, S., Fyllas, N. M., Galbraith, D. R., Baker,
! T. R., Kruijt, B., Rowland, L., Fisher, R. A., Binks, O. J., Sevanto, S., Xu, C., Jansen,
! S., Choat, B., Mencuccini, M., McDowell, N. G., Meir, P. Linking hydraulic traits to
! tropical forest function in a size-structured and trait-driven model (TFS~v.1-Hydro).
! Geoscientific Model Development, 9(11), 2016, pp: 4227-4255,
! https://www.geosci-model-dev.net/9/4227/2016/, DOI = 10.5194/gmd-9-4227-2016.
!
! WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
!
! PLANT HYDRAULICS IS AN EXPERIMENTAL OPTION THAT IS STILL UNDERGOING TESTING.
!
! WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
!
! ==============================================================================================
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!99
! (TODO: THE ROW WIDTH ON THIS MODULE ARE TOO LARGE. NAG COMPILERS
! WILL FREAK IF LINES ARE TOO LONG. BEFORE SUBMITTING THIS TO
! MASTER WE NEED TO GO THROUGH AND GET THESE LINES BELOW
! 100 spaces (for readability), or 130 (for NAG)
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!99
use FatesGlobals, only : endrun => fates_endrun
use FatesGlobals, only : fates_log
use FatesConstantsMod, only : r8 => fates_r8
use FatesConstantsMod, only : fates_huge
use FatesConstantsMod, only : denh2o => dens_fresh_liquid_water
use FatesConstantsMod, only : grav => grav_earth
use FatesConstantsMod, only : ifalse, itrue
use FatesConstantsMod, only : pi_const
use FatesConstantsMod, only : cm2_per_m2
use FatesConstantsMod, only : g_per_kg
use FatesConstantsMod, only : nearzero
use EDParamsMod , only : hydr_kmax_rsurf1
use EDParamsMod , only : hydr_kmax_rsurf2
use EDTypesMod , only : ed_site_type
use EDTypesMod , only : ed_patch_type
use EDTypesMod , only : ed_cohort_type
use FatesInterfaceMod , only : bc_in_type
use FatesInterfaceMod , only : bc_out_type
use FatesInterfaceMod , only : hlm_use_planthydro
use FatesAllometryMod, only : bsap_allom
use FatesAllometryMod, only : CrownDepth
use FatesAllometryMod , only : set_root_fraction
use FatesAllometryMod , only : i_hydro_rootprof_context
use FatesHydraulicsMemMod, only: ed_site_hydr_type
use FatesHydraulicsMemMod, only: ed_cohort_hydr_type
use FatesHydraulicsMemMod, only: n_hypool_leaf
use FatesHydraulicsMemMod, only: n_hypool_tot
use FatesHydraulicsMemMod, only: n_hypool_stem
use FatesHydraulicsMemMod, only: numLWPmem
use FatesHydraulicsMemMod, only: n_hypool_troot
use FatesHydraulicsMemMod, only: n_hypool_aroot
use FatesHydraulicsMemMod, only: n_porous_media
use FatesHydraulicsMemMod, only: nshell
use FatesHydraulicsMemMod, only: n_hypool_ag
use FatesHydraulicsMemMod, only: porous_media
use FatesHydraulicsMemMod, only: cap_slp
use FatesHydraulicsMemMod, only: cap_int
use FatesHydraulicsMemMod, only: cap_corr
use FatesHydraulicsMemMod, only: rwcft
use FatesHydraulicsMemMod, only: C2B
use FatesHydraulicsMemMod, only: InitHydraulicsDerived
use FatesHydraulicsMemMod, only: nlevsoi_hyd_max
use FatesHydraulicsMemMod, only: cohort_recruit_water_layer
use FatesHydraulicsMemMod, only: recruit_water_avail_layer
use PRTGenericMod, only : all_carbon_elements
use PRTGenericMod, only : leaf_organ, fnrt_organ, sapw_organ
use PRTGenericMod, only : store_organ, repro_organ, struct_organ
use clm_time_manager , only : get_step_size, get_nstep
use EDPftvarcon, only : EDPftvarcon_inst
! CIME Globals
use shr_log_mod , only : errMsg => shr_log_errMsg
use shr_infnan_mod , only : isnan => shr_infnan_isnan
implicit none
private
integer, parameter :: van_genuchten = 1
integer, parameter :: campbell = 2
integer :: iswc = campbell
! 1=leaf, 2=stem, 3=troot, 4=aroot
! Several of these may be better transferred to the parameter file in due time (RGK)
integer, public :: use_ed_planthydraulics = 1 ! 0 => use vanilla btran
! 1 => use BC hydraulics;
! 2 => use CX hydraulics
logical, public :: do_dqtopdth_leaf = .false. ! should a nonzero dqtopdth_leaf
! term be applied to the plant
! hydraulics numerical solution?
logical, public :: do_dyn_xylemrefill = .false. ! should the dynamics of xylem refilling
! (i.e., non-instantaneous) be considered
! within plant hydraulics?
logical, public :: do_kbound_upstream = .true. ! should the hydraulic conductance at the
! boundary between nodes be taken to be a
! function of the upstream loss of
! conductivity (flc)?
logical, public :: do_growthrecruiteffects = .true. ! should size- or root length-dependent
! hydraulic properties and states be
! updated every day when trees grow or
! when recruitment happens?
logical,parameter :: debug = .false. !flag to report warning in hydro
character(len=*), parameter, private :: sourcefile = &
__FILE__
! We use this parameter as the value for which we set un-initialized values
real(r8), parameter :: un_initialized = -9.9e32_r8
!
! !PUBLIC MEMBER FUNCTIONS:
public :: AccumulateMortalityWaterStorage
public :: RecruitWaterStorage
public :: hydraulics_drive
public :: InitHydrSites
public :: HydrSiteColdStart
public :: BTranForHLMDiagnosticsFromCohortHydr
public :: InitHydrCohort
public :: DeallocateHydrCohort
public :: UpdateH2OVeg
public :: CopyCohortHydraulics
public :: FuseCohortHydraulics
public :: updateSizeDepTreeHydProps
public :: updateWaterDepTreeHydProps
public :: updateSizeDepTreeHydStates
public :: initTreeHydStates
public :: updateSizeDepRhizHydProps
public :: updateSizeDepRhizHydStates
public :: RestartHydrStates
public :: SavePreviousCompartmentVolumes
public :: SavePreviousRhizVolumes
public :: UpdateTreeHydrNodes
public :: UpdateTreeHydrLenVolCond
public :: UpdateWaterDepTreeHydrCond
public :: ConstrainRecruitNumber
!------------------------------------------------------------------------------
! 01/18/16: Created by Brad Christoffersen
! 02/xx/17: Refactoring by Ryan Knox and Brad Christoffersen
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
subroutine hydraulics_drive( nsites, sites, bc_in,bc_out,dtime )
! ARGUMENTS:
! -----------------------------------------------------------------------------------
integer,intent(in) :: nsites
type(ed_site_type),intent(inout),target :: sites(nsites)
type(bc_in_type),intent(in) :: bc_in(nsites)
type(bc_out_type),intent(inout) :: bc_out(nsites)
real(r8),intent(in) :: dtime
select case (use_ed_planthydraulics)
case (1)
call FillDrainRhizShells(nsites, sites, bc_in, bc_out )
call hydraulics_BC(nsites, sites,bc_in,bc_out,dtime )
case (2)
!call Hydraulics_CX()
case DEFAULT
end select
end subroutine Hydraulics_Drive
! =====================================================================================
subroutine RestartHydrStates(sites,nsites,bc_in,bc_out)
! It is assumed that the following state variables have been read in by
! the restart machinery.
!
! co_hydr%th_ag
! co_hydr%th_troot
! co_hydr%th_aroot
! si_hydr%r_node_shell
! si_hydr%v_shell
! si_hydr%h2osoi_liqvol_shell
! si_hydr%l_aroot_layer
!
! The goal of this subroutine is to call
! the correct sequence of hydraulics initializations to repopulate
! information that relies on these key states, as well as other vegetation
! states such as carbon pools and plant geometry.
integer , intent(in) :: nsites
type(ed_site_type) , intent(inout), target :: sites(nsites)
type(bc_in_type) , intent(in) :: bc_in(nsites)
type(bc_out_type) , intent(inout) :: bc_out(nsites)
! locals
! ----------------------------------------------------------------------------------
! LL pointers
type(ed_patch_type),pointer :: cpatch ! current patch
type(ed_cohort_type),pointer :: ccohort ! current cohort
type(ed_cohort_hydr_type),pointer :: ccohort_hydr
integer :: s ! site loop counter
do s = 1,nsites
cpatch => sites(s)%oldest_patch
do while(associated(cpatch))
ccohort => cpatch%shortest
do while(associated(ccohort))
ccohort_hydr => ccohort%co_hydr
! This calculates node heights
call UpdateTreeHydrNodes(ccohort_hydr,ccohort%pft,ccohort%hite, &
sites(s)%si_hydr%nlevsoi_hyd,bc_in(s))
! This calculates volumes, lengths and max conductances
call UpdateTreeHydrLenVolCond(ccohort,sites(s)%si_hydr%nlevsoi_hyd,bc_in(s))
! Since this is a newly initialized plant, we set the previous compartment-size
! equal to the ones we just calculated.
call SavePreviousCompartmentVolumes(ccohort_hydr)
! Set some generic initial values
ccohort_hydr%refill_days = 3.0_r8
ccohort_hydr%lwp_mem(:) = 0.0_r8
ccohort_hydr%lwp_stable = 0.0_r8
ccohort_hydr%lwp_is_unstable = .false.
ccohort_hydr%flc_ag(:) = 1.0_r8
ccohort_hydr%flc_troot(:) = 1.0_r8
ccohort_hydr%flc_aroot(:) = 1.0_r8
ccohort_hydr%flc_min_ag(:) = 1.0_r8
ccohort_hydr%flc_min_troot(:) = 1.0_r8
ccohort_hydr%flc_min_aroot(:) = 1.0_r8
ccohort_hydr%refill_thresh = -0.01_r8
ccohort_hydr%refill_days = 3.0_r8
ccohort => ccohort%taller
enddo
cpatch => cpatch%younger
end do
sites(s)%si_hydr%l_aroot_layer_init(:) = un_initialized
sites(s)%si_hydr%r_node_shell_init(:,:) = un_initialized
sites(s)%si_hydr%v_shell_init(:,:) = un_initialized
! Update static quantities related to the rhizosphere
call UpdateSizeDepRhizVolLenCon(sites(s), bc_in(s))
! We update the "initial" values of the rhizosphere after
! the previous call to make sure that the conductances are updated
! Now we set the prevous to the current so that the water states
! are not perturbed
call SavePreviousRhizVolumes(sites(s), bc_in(s))
end do
call UpdateH2OVeg(nsites,sites,bc_out)
return
end subroutine RestartHydrStates
! ====================================================================================
subroutine initTreeHydStates(site_p, cc_p, bc_in)
! REQUIRED INPUTS:
!
! csite%si_hydr%psisoi_liq_innershell(:)
! ccohort_hydr%z_node_troot(:)
! ccohort_hydr%z_node_aroot
! ccohort_hydr%z_node_ag
!
! !DESCRIPTION:
!
! !USES:
! !ARGUMENTS:
type(ed_site_type), intent(inout), target :: site_p ! current cohort pointer
type(ed_cohort_type), intent(inout), target :: cc_p ! current cohort pointer
type(bc_in_type) , intent(in) :: bc_in
!
! !LOCAL VARIABLES:
type(ed_cohort_type), pointer :: cCohort
type(ed_site_type), pointer :: csite
type(ed_cohort_hydr_type), pointer :: ccohort_hydr
integer :: j,k,FT ! indices
real(r8) :: dz
real(r8) :: smp
cCohort => cc_p
ccohort_hydr => cCohort%co_hydr
csite => site_p
FT = cCohort%pft
!convert soil water contents to water potential in each soil layer and
!assign it to the absorbing root (assume absorbing root water potential
!in equlibrium w/ surrounding soil)
do j=1, site_p%si_hydr%nlevsoi_hyd
!call swcVG_psi_from_th(waterstate_inst%h2osoi_liqvol_shell(c,j,1), &
! watsat(c,j), watres(c,j), alpha_VG(c,j), n_VG(c,j), m_VG(c,j), l_VG(c,j), &
! smp)
!ccohort_hydr%psi_aroot(j) = smp
!ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j)
ccohort_hydr%psi_aroot(j) = -0.2_r8 !do not assume the equalibrium between soil and root
call th_from_psi(ft, 4, ccohort_hydr%psi_aroot(j), ccohort_hydr%th_aroot(j), csite%si_hydr, bc_in )
end do
!initialize plant water potentials at hydrostatic equilibrium (dh/dz = 0)
!the assumption is made here that initial conditions for soil water will
!be in (or at least close to) hydrostatic equilibrium as well, so that
!it doesn't matter which absorbing root layer the transporting root water
!potential is referenced to.
do k=1, n_hypool_troot
dz = ccohort_hydr%z_node_troot(k) - ccohort_hydr%z_node_aroot(1)
ccohort_hydr%psi_troot(k) = ccohort_hydr%psi_aroot(1) - 1.e-6_r8*denh2o*grav*dz
if (ccohort_hydr%psi_troot(k)>0.0_r8) ccohort_hydr%psi_troot(k) = -0.01_r8
call th_from_psi(ft, 3, ccohort_hydr%psi_troot(k), ccohort_hydr%th_troot(k), csite%si_hydr, bc_in)
end do
!working our way up a tree, assigning water potentials that are in
!hydrostatic equilibrium with the water potential immediately below
dz = ccohort_hydr%z_node_ag(n_hypool_ag) - ccohort_hydr%z_node_troot(1)
ccohort_hydr%psi_ag(n_hypool_ag) = ccohort_hydr%psi_troot(1) - 1.e-6_r8*denh2o*grav*dz
if (ccohort_hydr%psi_ag(n_hypool_ag)>0.0_r8) ccohort_hydr%psi_ag(n_hypool_ag) = -0.01_r8
call th_from_psi(ft, 2, ccohort_hydr%psi_ag(n_hypool_ag), ccohort_hydr%th_ag(n_hypool_ag), csite%si_hydr, bc_in)
do k=n_hypool_ag-1, 1, -1
dz = ccohort_hydr%z_node_ag(k) - ccohort_hydr%z_node_ag(k+1)
ccohort_hydr%psi_ag(k) = ccohort_hydr%psi_ag(k+1) - 1.e-6_r8*denh2o*grav*dz
if(ccohort_hydr%psi_ag(k)>0.0_r8) ccohort_hydr%psi_ag(k)= -0.01_r8
call th_from_psi(ft, porous_media(k), ccohort_hydr%psi_ag(k), ccohort_hydr%th_ag(k), csite%si_hydr, bc_in)
end do
ccohort_hydr%lwp_mem(:) = ccohort_hydr%psi_ag(1) ! initializes the leaf water potential memory
ccohort_hydr%lwp_stable = ccohort_hydr%psi_ag(1)
ccohort_hydr%lwp_is_unstable = .false. ! inital value for leaf water potential stability flag
ccohort_hydr%flc_ag(:) = 1.0_r8
ccohort_hydr%flc_troot(:) = 1.0_r8
ccohort_hydr%flc_aroot(:) = 1.0_r8
ccohort_hydr%flc_min_ag(:) = 1.0_r8
ccohort_hydr%flc_min_troot(:) = 1.0_r8
ccohort_hydr%flc_min_aroot(:) = 1.0_r8
ccohort_hydr%refill_thresh = -0.01_r8
ccohort_hydr%refill_days = 3.0_r8
ccohort_hydr%errh2o_growturn_ag(:) = 0.0_r8
ccohort_hydr%errh2o_growturn_troot(:) = 0.0_r8
ccohort_hydr%errh2o_growturn_aroot(:) = 0.0_r8
ccohort_hydr%errh2o_pheno_ag(:) = 0.0_r8
ccohort_hydr%errh2o_pheno_troot(:) = 0.0_r8
ccohort_hydr%errh2o_pheno_aroot(:) = 0.0_r8
!ccohort_hydr%th_aroot_prev(:) = 0.0_r8
!ccohort_hydr%th_aroot_prev_ucnorr(:)= 0.0_r8
!initialize cohort-level btran
call flc_gs_from_psi(cCohort, ccohort_hydr%psi_ag(1))
end subroutine initTreeHydStates
! =====================================================================================
subroutine UpdateTreeHydrNodes(ccohort_hydr,ft,plant_height,nlevsoi_hyd,bc_in)
! --------------------------------------------------------------------------------
! This subroutine calculates the nodal heights critical to hydraulics in the plant
!
! Inputs: Plant height
! Plant functional type
! Number of soil hydraulic layers
!
! Outputs: cohort_hydr%z_node_ag(:)
! %z_lower_ag(:)
! %z_upper_ag(:)
! %z_node_troot(:)
! %z_lower_troot(:)
! %z_upper_troot(:)
! %z_node_aroot(:)
! --------------------------------------------------------------------------------
! Arguments
type(ed_cohort_hydr_type), intent(inout) :: ccohort_hydr
integer,intent(in) :: ft ! plant functional type index
real(r8), intent(in) :: plant_height ! [m]
integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions
! Locals
real(r8) :: roota ! root profile parameter a zeng2001_crootfr
real(r8) :: rootb ! root profile parameter b zeng2001_crootfr
real(r8) :: crown_depth ! crown depth for the plant [m]
real(r8) :: dz_canopy ! discrete crown depth intervals [m]
real(r8) :: z_stem ! the height of the plants stem below crown [m]
real(r8) :: dz_stem ! vertical stem discretization [m]
real(r8) :: dcumul_rf ! cumulative root distribution discretization [-]
real(r8) :: cumul_rf ! cumulative root distribution where depth is determined [-]
real(r8) :: z_cumul_rf ! depth at which cumul_rf occurs [m]
integer :: k ! Loop counter for compartments
! Crown Nodes
! in special case where n_hypool_leaf = 1, the node height of the canopy
! water pool is 1/2 the distance from the bottom of the canopy to the top of the tree
roota = EDPftvarcon_inst%roota_par(ft)
rootb = EDPftvarcon_inst%rootb_par(ft)
call CrownDepth(plant_height,crown_depth)
dz_canopy = crown_depth / real(n_hypool_leaf,r8)
do k=1,n_hypool_leaf
ccohort_hydr%z_lower_ag(k) = plant_height - dz_canopy*real(k,r8)
ccohort_hydr%z_node_ag(k) = ccohort_hydr%z_lower_ag(k) + 0.5_r8*dz_canopy
ccohort_hydr%z_upper_ag(k) = ccohort_hydr%z_lower_ag(k) + dz_canopy
enddo
! Stem Nodes
! in special case where n_hypool_stem = 1, the node height of the stem water pool is
! 1/2 the height from the ground to the bottom of the canopy
z_stem = plant_height - crown_depth
dz_stem = z_stem / real(n_hypool_stem,r8)
do k=n_hypool_leaf+1,n_hypool_ag
ccohort_hydr%z_upper_ag(k) = real(n_hypool_stem - (k - 1 - n_hypool_leaf),r8)*dz_stem
ccohort_hydr%z_node_ag(k) = ccohort_hydr%z_upper_ag(k) - 0.5_r8*dz_stem
ccohort_hydr%z_lower_ag(k) = ccohort_hydr%z_upper_ag(k) - dz_stem
enddo
! Transporting Root Nodes
! in special case where n_hypool_troot = 1, the node depth of the single troot pool
! is the depth at which 50% total root distribution is attained
dcumul_rf = 1._r8/real(n_hypool_troot,r8)
do k=1,n_hypool_troot
cumul_rf = dcumul_rf*real(k,r8)
call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, &
0.001_r8, 0.001_r8, cumul_rf, z_cumul_rf)
z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd)))
ccohort_hydr%z_lower_troot(k) = -z_cumul_rf
call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, &
0.001_r8, 0.001_r8, cumul_rf-0.5_r8*dcumul_rf, z_cumul_rf)
z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd)))
ccohort_hydr%z_node_troot(k) = -z_cumul_rf
call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, &
0.001_r8, 0.001_r8, cumul_rf-1.0_r8*dcumul_rf+1.E-10_r8, z_cumul_rf)
z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd)))
ccohort_hydr%z_upper_troot(k) = -z_cumul_rf
enddo
! Absorbing root depth
ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd)
! Shouldn't this be updating the upper and lower values as well?
! (RGK 12-2018)
if(nlevsoi_hyd == 1) then
ccohort_hydr%z_node_troot(:) = ccohort_hydr%z_node_aroot(nlevsoi_hyd)
end if
return
end subroutine UpdateTreeHydrNodes
! =====================================================================================
subroutine SavePreviousCompartmentVolumes(ccohort_hydr)
type(ed_cohort_hydr_type),intent(inout) :: ccohort_hydr
! Saving the current compartment volumes into an "initial" save-space
! allows us to see how the compartments change size when plants
! change size and effect water contents
ccohort_hydr%v_ag_init(:) = ccohort_hydr%v_ag(:)
ccohort_hydr%v_troot_init(:) = ccohort_hydr%v_troot(:)
ccohort_hydr%v_aroot_layer_init(:) = ccohort_hydr%v_aroot_layer(:)
return
end subroutine SavePreviousCompartmentVolumes
! =====================================================================================
subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in)
! DESCRIPTION: Updates absorbing root length (total and its vertical distribution)
! as well as the consequential change in the size of the 'representative' rhizosphere
! shell radii, volumes, and compartment volumes of plant tissues
! !USES:
use shr_sys_mod , only : shr_sys_abort
! ARGUMENTS:
type(ed_site_type) , intent(in) :: currentSite ! Site stuff
type(ed_cohort_type) , intent(inout) :: ccohort ! current cohort pointer
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions
! Locals
integer :: nlevsoi_hyd ! Number of total soil layers
type(ed_cohort_hydr_type), pointer :: ccohort_hydr
integer :: ft
nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd
ccohort_hydr => ccohort%co_hydr
ft = ccohort%pft
! Save the current vegetation compartment volumes into
! a save space so that it can be compared with the updated quantity.
call SavePreviousCompartmentVolumes(ccohort_hydr)
! This updates all of the z_node positions
call UpdateTreeHydrNodes(ccohort_hydr,ft,ccohort%hite,nlevsoi_hyd,bc_in)
! This updates plant compartment volumes, lengths and
! maximum conductances. Make sure for already
! initialized vegetation, that SavePreviousCompartment
! volumes, and UpdateTreeHydrNodes is called prior to this.
call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
end subroutine updateSizeDepTreeHydProps
! =====================================================================================
subroutine updateWaterDepTreeHydProps(currentSite,ccohort,bc_in)
! DESCRIPTION: Updates absorbing root length (total and its vertical distribution)
! as well as the consequential change in the size of the 'representative' rhizosphere
! shell radii, volumes, and compartment volumes of plant tissues
! !USES:
use shr_sys_mod , only : shr_sys_abort
! ARGUMENTS:
type(ed_site_type) , intent(in) :: currentSite ! Site stuff
type(ed_cohort_type) , intent(inout) :: ccohort ! current cohort pointer
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions
! Locals
integer :: nlevsoi_hyd ! Number of total soil layers
type(ed_cohort_hydr_type), pointer :: ccohort_hydr
integer :: ft
nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd
ccohort_hydr => ccohort%co_hydr
ft = ccohort%pft
! This updates plant compartment volumes, lengths and
! maximum conductances. Make sure for already
! initialized vegetation, that SavePreviousCompartment
! volumes, and UpdateTreeHydrNodes is called prior to this.
call UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in)
end subroutine updateWaterDepTreeHydProps
! =====================================================================================
subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
! -----------------------------------------------------------------------------------
! This subroutine calculates three attributes of a plant:
! 1) the volumes of storage compartments in the plants
! 2) the lenghts of the organs
! 3) the conductances
! These and are not dependent on the hydraulic state of the
! plant, it is more about the structural characteristics and how much biomass
! is present in the different tissues.
!
! Inputs, plant geometries, plant carbon pools, z_node values
!
! -----------------------------------------------------------------------------------
! Arguments
type(ed_cohort_type),intent(inout) :: ccohort
integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions
type(ed_cohort_hydr_type),pointer :: ccohort_hydr ! Plant hydraulics structure
integer :: j,k
integer :: ft ! Plant functional type index
real(r8) :: roota ! root profile parameter a zeng2001_crootfr
real(r8) :: rootb ! root profile parameter b zeng2001_crootfr
real(r8) :: leaf_c ! Current amount of leaf carbon in the plant [kg]
real(r8) :: fnrt_c ! Current amount of fine-root carbon in the plant [kg]
real(r8) :: sapw_c ! Current amount of sapwood carbon in the plant [kg]
real(r8) :: struct_c ! Current amount of structural carbon in the plant [kg]
real(r8) :: b_canopy_carb ! total leaf (canopy) biomass in carbon units [kgC/indiv]
real(r8) :: b_canopy_biom ! total leaf (canopy) biomass in dry wt units [kg/indiv]
real(r8) :: b_woody_carb ! total woody biomass in carbon units [kgC/indiv]
real(r8) :: b_woody_bg_carb ! belowground woody biomass in carbon units [kgC/indiv]
real(r8) :: b_stem_carb ! aboveground stem biomass in carbon units [kgC/indiv]
real(r8) :: b_stem_biom ! aboveground stem biomass in dry wt units [kg/indiv]
real(r8) :: b_bg_carb ! belowground biomass (coarse + fine roots) in carbon units [kgC/indiv]
real(r8) :: b_tot_carb ! total individual biomass in carbon units [kgC/indiv]
real(r8) :: v_stem ! aboveground stem volume [m3/indiv]
real(r8) :: z_stem ! the height of the plants stem below crown [m]
real(r8) :: sla ! specific leaf area [cm2/g]
real(r8) :: v_canopy ! total leaf (canopy) volume [m3/indiv]
real(r8) :: denleaf ! leaf dry mass per unit fresh leaf volume [kg/m3]
real(r8) :: a_sapwood ! sapwood area [m2]
real(r8) :: a_sapwood_target ! sapwood cross-section area at reference height, at target biomass [m2]
real(r8) :: bsw_target ! sapwood carbon, at target [kgC]
real(r8) :: v_sapwood ! sapwood volume [m3]
real(r8) :: b_troot_carb ! transporting root biomass in carbon units [kgC/indiv]
real(r8) :: b_troot_biom ! transporting root biomass in dry wt units [kg/indiv]
real(r8) :: v_troot ! transporting root volume [m3/indiv]
real(r8) :: rootfr ! mass fraction of roots in each layer [kg/kg]
real(r8), allocatable :: rootfrs(:) ! Vector of root fractions (only used in 1 layer case) [kg/kg]
real(r8) :: crown_depth ! Depth of the plant's crown [m]
real(r8) :: kmax_node1_nodekplus1(n_hypool_ag) ! cumulative kmax, petiole to node k+1,
! conduit taper effects excluded [kg s-1 MPa-1]
real(r8) :: kmax_node1_lowerk(n_hypool_ag) ! cumulative kmax, petiole to upper boundary of node k,
! conduit taper effects excluded [kg s-1 MPa-1]
real(r8) :: chi_node1_nodekplus1(n_hypool_ag) ! ratio of cumulative kmax with taper effects
! included to that without [-]
real(r8) :: chi_node1_lowerk(n_hypool_ag) ! ratio of cumulative kmax with taper effects
! included to that without [-]
real(r8) :: dz_node1_nodekplus1 ! cumulative distance between canopy
! node and node k + 1 [m]
real(r8) :: dz_node1_lowerk ! cumulative distance between canopy
! node and upper boundary of node k [m]
real(r8) :: kmax_treeag_tot ! total stem (petiole to transporting root node)
! hydraulic conductance [kg s-1 MPa-1]
real(r8) :: kmax_tot ! total tree (leaf to root tip)
! hydraulic conductance [kg s-1 MPa-1]
real(r8),parameter :: taper_exponent = 1._r8/3._r8 ! Savage et al. (2010) xylem taper exponent [-]
ccohort_hydr => ccohort%co_hydr
ft = ccohort%pft
leaf_c = ccohort%prt%GetState(leaf_organ, all_carbon_elements)
sapw_c = ccohort%prt%GetState(sapw_organ, all_carbon_elements)
fnrt_c = ccohort%prt%GetState(fnrt_organ, all_carbon_elements)
struct_c = ccohort%prt%GetState(struct_organ, all_carbon_elements)
roota = EDPftvarcon_inst%roota_par(ft)
rootb = EDPftvarcon_inst%rootb_par(ft)
!roota = 4.372_r8 ! TESTING: deep (see Zeng 2001 Table 1)
!rootb = 0.978_r8 ! TESTING: deep (see Zeng 2001 Table 1)
!roota = 8.992_r8 ! TESTING: shallow (see Zeng 2001 Table 1)
!rootb = 8.992_r8 ! TESTING: shallow (see Zeng 2001 Table 1)
if(leaf_c > 0._r8) then
! ------------------------------------------------------------------------------
! Part 1. Set the volumes of the leaf, stem and root compartments
! and lenghts of the roots
! ------------------------------------------------------------------------------
b_woody_carb = sapw_c + struct_c
b_woody_bg_carb = (1.0_r8-EDPftvarcon_inst%allom_agb_frac(ft)) * b_woody_carb
b_tot_carb = sapw_c + struct_c + leaf_c + fnrt_c
b_canopy_carb = leaf_c
b_bg_carb = (1.0_r8-EDPftvarcon_inst%allom_agb_frac(ft)) * b_tot_carb
b_canopy_biom = b_canopy_carb * C2B
! NOTE: SLATOP currently does not use any vertical scaling functions
! but that may not be so forever. ie sla = slatop (RGK-082017)
! m2/gC * cm2/m2 -> cm2/gC
sla = EDPftvarcon_inst%slatop(ft) * cm2_per_m2
! empirical regression data from leaves at Caxiuana (~ 8 spp)
denleaf = -2.3231_r8*sla/C2B + 781.899_r8
v_canopy = b_canopy_biom / denleaf
ccohort_hydr%v_ag(1:n_hypool_leaf) = v_canopy / real(n_hypool_leaf,r8)
b_stem_carb = b_tot_carb - b_bg_carb - b_canopy_carb
b_stem_biom = b_stem_carb * C2B ! kg DM
!BOC...may be needed for testing/comparison w/ v_sapwood
! kg / ( g cm-3 * cm3/m3 * kg/g ) -> m3
v_stem = b_stem_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8 )
! calculate the sapwood cross-sectional area
call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%canopy_trim,a_sapwood_target,bsw_target)
a_sapwood = a_sapwood_target
! Alternative ways to calculate sapwood cross section
! or ....
! a_sapwood = a_sapwood_target * ccohort%bsw / bsw_target
! a_sapwood = a_leaf_tot / EDPftvarcon_inst%allom_latosa_int(ft)*1.e-4_r8
! m2 sapwood = m2 leaf * cm2 sapwood/m2 leaf *1.0e-4m2
! or ...
!a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * ccohort%hite ) * 1.e-4_r8
call CrownDepth(ccohort%hite,crown_depth)
z_stem = ccohort%hite - crown_depth
v_sapwood = a_sapwood * z_stem
ccohort_hydr%v_ag(n_hypool_leaf+1:n_hypool_ag) = v_sapwood / n_hypool_stem
! Determine belowground biomass as a function of total (sapwood, heartwood,
! leaf, fine root) biomass then subtract out the fine root biomass to get
! coarse (transporting) root biomass
b_troot_carb = b_woody_bg_carb
b_troot_biom = b_troot_carb * C2B
v_troot = b_troot_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8)
!! BOC not sure if/how we should multiply this by the sapwood fraction
ccohort_hydr%v_troot(:) = v_troot / n_hypool_troot
! Estimate absorbing root total length (all layers)
! ------------------------------------------------------------------------------
ccohort_hydr%l_aroot_tot = fnrt_c*C2B*EDPftvarcon_inst%hydr_srl(ft)
! Estimate absorbing root volume (all layers)
! ------------------------------------------------------------------------------
ccohort_hydr%v_aroot_tot = pi_const * (EDPftvarcon_inst%hydr_rs2(ft)**2._r8) * &
ccohort_hydr%l_aroot_tot
! Partition the total absorbing root lengths and volumes into the active soil layers
! ------------------------------------------------------------------------------
if(nlevsoi_hyd == 1) then
ccohort_hydr%l_aroot_layer(nlevsoi_hyd) = ccohort_hydr%l_aroot_tot
ccohort_hydr%v_aroot_layer(nlevsoi_hyd) = ccohort_hydr%v_aroot_tot
else
do j=1,nlevsoi_hyd
if(j == 1) then
rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j))
else
rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - &
zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j-1))
end if
ccohort_hydr%l_aroot_layer(j) = rootfr*ccohort_hydr%l_aroot_tot
ccohort_hydr%v_aroot_layer(j) = rootfr*ccohort_hydr%v_aroot_tot
end do
end if
! ------------------------------------------------------------------------------
! Part II. Set maximum (size-dependent) hydraulic conductances
! ------------------------------------------------------------------------------
! first estimate cumulative (petiole to node k) conductances
! without taper as well as the chi taper function
do k=n_hypool_leaf,n_hypool_ag
dz_node1_lowerk = ccohort_hydr%z_node_ag(n_hypool_leaf) &
- ccohort_hydr%z_lower_ag(k)
if(k < n_hypool_ag) then
dz_node1_nodekplus1 = ccohort_hydr%z_node_ag(n_hypool_leaf) &
- ccohort_hydr%z_node_ag(k+1)
else
dz_node1_nodekplus1 = ccohort_hydr%z_node_ag(n_hypool_leaf) &
- ccohort_hydr%z_node_troot(1)
end if
kmax_node1_nodekplus1(k) = EDPftvarcon_inst%hydr_kmax_node(ft,2) * a_sapwood / dz_node1_nodekplus1
kmax_node1_lowerk(k) = EDPftvarcon_inst%hydr_kmax_node(ft,2) * a_sapwood / dz_node1_lowerk
chi_node1_nodekplus1(k) = xylemtaper(taper_exponent, dz_node1_nodekplus1)
chi_node1_lowerk(k) = xylemtaper(taper_exponent, dz_node1_lowerk)
if(.not.do_kbound_upstream) then
if(crown_depth == 0._r8) then
write(fates_log(),*) 'do_kbound_upstream requires a nonzero canopy depth '
call endrun(msg=errMsg(sourcefile, __LINE__))
end if
end if
enddo
! then calculate the conductances at node boundaries as the difference of cumulative conductances
do k=n_hypool_leaf,n_hypool_ag
if(k == n_hypool_leaf) then
ccohort_hydr%kmax_bound(k) = kmax_node1_nodekplus1(k) * chi_node1_nodekplus1(k)
ccohort_hydr%kmax_lower(k) = kmax_node1_lowerk(k) * chi_node1_lowerk(k)
else
ccohort_hydr%kmax_bound(k) = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - &
1._r8/(kmax_node1_nodekplus1(k-1)*chi_node1_nodekplus1(k-1)) ) ** (-1._r8)
ccohort_hydr%kmax_lower(k) = ( 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) - &
1._r8/(kmax_node1_nodekplus1(k-1)*chi_node1_nodekplus1(k-1)) ) ** (-1._r8)
end if
if(k < n_hypool_ag) then
ccohort_hydr%kmax_upper(k+1) = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - &
1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) ) ** (-1._r8)
else if(k == n_hypool_ag) then
ccohort_hydr%kmax_upper_troot = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - &
1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) ) ** (-1._r8)
end if
!!!!!!!!!! FOR TESTING ONLY
!ccohort_hydr%kmax_bound(:) = 0.02_r8 ! Diurnal lwp variation in coldstart: -0.1 MPa
! Diurnal lwp variation in large-tree (50cmDBH) coldstart: less than -0.01 MPa
!ccohort_hydr%kmax_bound(:) = 0.0016_r8 ! Diurnal lwp variation in coldstart: -0.8 - 1.0 MPa
! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -1.5 - 2.0 MPa [seemingly unstable]
!ccohort_hydr%kmax_bound(:) = 0.0008_r8 ! Diurnal lwp variation in coldstart: -1.5 - 2.0 MPa
! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -2.0 - 3.0 MPa [seemingly unstable]
!ccohort_hydr%kmax_bound(:) = 0.0005_r8 ! Diurnal lwp variation in coldstart: -2.0 - 3.0 MPa and one -5 MPa outlier
! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -3.0 - 4.0 MPa and one -10 MPa outlier [Unstable]
!!!!!!!!!!
enddo
! finally, estimate the remaining tree conductance belowground as a residual
kmax_treeag_tot = sum(1._r8/ccohort_hydr%kmax_bound(n_hypool_leaf:n_hypool_ag))**(-1._r8)
kmax_tot = EDPftvarcon_inst%hydr_rfrac_stem(ft) * kmax_treeag_tot
ccohort_hydr%kmax_treebg_tot = ( 1._r8/kmax_tot - 1._r8/kmax_treeag_tot ) ** (-1._r8)
if(nlevsoi_hyd == 1) then
allocate(rootfrs(bc_in%nlevsoil))
call set_root_fraction(rootfrs(:), ft, bc_in%zi_sisl, &
icontext = i_hydro_rootprof_context)
ccohort_hydr%kmax_treebg_layer(:) = ccohort_hydr%kmax_treebg_tot * rootfrs(:)
deallocate(rootfrs)
else
do j=1,nlevsoi_hyd
if(j == 1) then
rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j))
else
rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - &
zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j-1))
end if
ccohort_hydr%kmax_treebg_layer(j) = rootfr*ccohort_hydr%kmax_treebg_tot
end do
end if
end if !check for bleaf
end subroutine UpdateTreeHydrLenVolCond
!=====================================================================================
subroutine UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in)
! -----------------------------------------------------------------------------------
! This subroutine calculates update the conductivity for the soil-root interface,
! depending on the plant water uptake/loss.
! we assume that the conductivitity for water uptake is larger than
! water loss due to composite regulation of resistance the roots
! hydraulic vs osmostic with and without transpiration
! Steudle, E. Water uptake by roots: effects of water deficit.
! J Exp Bot 51, 1531-1542, doi:DOI 10.1093/jexbot/51.350.1531 (2000).
! -----------------------------------------------------------------------------------
! Arguments
type(ed_site_type) , intent(in) :: currentSite ! Site target
type(ed_cohort_type),intent(inout) :: ccohort ! cohort target
integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions
type(ed_cohort_hydr_type),pointer :: ccohort_hydr ! Plant hydraulics structure
type(ed_site_hydr_type),pointer :: csite_hydr
integer :: j,k
real(r8) :: hksat_s ! hksat converted to units of 10^6sec
real(r8) :: kmax_root_surf_total ! maximum conducitivity for total root surface(kg water/Mpa/s)
real(r8) :: kmax_soil_total ! maximum conducitivity for from root surface to soil shell(kg water/Mpa/s)
! which is equiv to [kg m-1 s-1 MPa-1]
real(r8) :: kmax_root_surf ! maximum conducitivity for unit root surface (kg water/m2 root area/Mpa/s)
ccohort_hydr => ccohort%co_hydr
csite_hydr => currentSite%si_hydr
k = 1 !only for the first soil shell
do j=1, nlevsoi_hyd
hksat_s = bc_in%hksat_sisl(j) * 1.e-3_r8 * 1/grav * 1.e6_r8
if(ccohort_hydr%psi_aroot(j)<csite_hydr%psisoi_liq_innershell(j))then
kmax_root_surf = hydr_kmax_rsurf1
else
kmax_root_surf = hydr_kmax_rsurf2
endif
kmax_root_surf_total = kmax_root_surf*2._r8*pi_const *csite_hydr%rs1(j)* &
csite_hydr%l_aroot_layer(j)
if(csite_hydr%r_node_shell(j,1) <= csite_hydr%rs1(j)) then
!csite_hydr%kmax_upper_shell(j,k) = large_kmax_bound
!csite_hydr%kmax_bound_shell(j,k) = large_kmax_bound
!csite_hydr%kmax_lower_shell(j,k) = large_kmax_bound
ccohort_hydr%kmax_innershell(j) = kmax_root_surf_total
else
kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_bound_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_lower_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
ccohort_hydr%kmax_innershell(j) = (1._r8/kmax_root_surf_total + &
1._r8/kmax_soil_total)**(-1._r8)
end if
end do
end subroutine UpdateWaterDepTreeHydrCond
! =====================================================================================
subroutine updateSizeDepTreeHydStates(currentSite,ccohort)
!
! !DESCRIPTION:
!
! !USES:
use EDTypesMod , only : AREA
! !ARGUMENTS:
type(ed_site_type) , intent(in) :: currentSite ! Site stuff
type(ed_cohort_type) , intent(inout) :: ccohort
!
! !LOCAL VARIABLES:
type(ed_cohort_hydr_type), pointer :: ccohort_hydr
type(ed_site_hydr_type),pointer :: csite_hydr
integer :: j,k,FT ! indices
integer :: err_code = 0
real(r8) :: th_ag_uncorr( n_hypool_ag) ! uncorrected aboveground water content[m3 m-3]
real(r8) :: th_troot_uncorr(n_hypool_troot) ! uncorrected transporting root water content[m3 m-3]
real(r8) :: th_aroot_uncorr(currentSite%si_hydr%nlevsoi_hyd) ! uncorrected absorbing root water content[m3 m-3]
real(r8), parameter :: small_theta_num = 1.e-7_r8 ! avoids theta values equalling thr or ths [m3 m-3]
integer :: nstep !number of time steps
!-----------------------------------------------------------------------
ccohort_hydr => ccohort%co_hydr
FT = cCohort%pft
! MAYBE ADD A NAN CATCH? If updateSizeDepTreeHydProps() was not called twice prior to the first
! time this routine is called for a new cohort, then v_ag_init(k) will be a nan.
! It should be ok, but may be vulnerable if code is changed (RGK 02-2017)
! UPDATE WATER CONTENTS (assume water for growth comes from within tissue itself -- apply water mass conservation)
do k=1,n_hypool_ag
th_ag_uncorr(k) = ccohort_hydr%th_ag(k) * &
ccohort_hydr%v_ag_init(k) /ccohort_hydr%v_ag(k)
ccohort_hydr%th_ag(k) = constrain_water_contents(th_ag_uncorr(k), small_theta_num, ft, k)
enddo
do k=1,n_hypool_troot
th_troot_uncorr(k) = ccohort_hydr%th_troot(k) * &
ccohort_hydr%v_troot_init(k) /ccohort_hydr%v_troot(k)
ccohort_hydr%th_troot(k) = constrain_water_contents(th_troot_uncorr(k), small_theta_num, ft, 3)
enddo
do j=1,currentSite%si_hydr%nlevsoi_hyd
th_aroot_uncorr(j) = ccohort_hydr%th_aroot(j) * &
ccohort_hydr%v_aroot_layer_init(j)/ccohort_hydr%v_aroot_layer(j)
ccohort_hydr%th_aroot(j) = constrain_water_contents(th_aroot_uncorr(j), small_theta_num, ft, 4)
ccohort_hydr%errh2o_growturn_aroot(j) = ccohort_hydr%th_aroot(j) - th_aroot_uncorr(j)
enddo
! Storing mass balance error
! + means water created; - means water destroyed
ccohort_hydr%errh2o_growturn_ag(:) = ccohort_hydr%th_ag(:) - th_ag_uncorr(:)
ccohort_hydr%errh2o_growturn_troot(:) = ccohort_hydr%th_troot(:) - th_troot_uncorr(:)
csite_hydr =>currentSite%si_hydr
csite_hydr%h2oveg_growturn_err = csite_hydr%h2oveg_growturn_err + &
(sum(ccohort_hydr%errh2o_growturn_ag(:)*ccohort_hydr%v_ag(:)) + &
sum(ccohort_hydr%errh2o_growturn_troot(:)*ccohort_hydr%v_troot(:)) + &
sum(ccohort_hydr%errh2o_growturn_aroot(:)*ccohort_hydr%v_aroot_layer(:)))* &
denh2o*cCohort%n/AREA
! UPDATES OF WATER POTENTIALS ARE DONE PRIOR TO RICHARDS' SOLUTION WITHIN FATESPLANTHYDRAULICSMOD.F90
end subroutine updateSizeDepTreeHydStates
! =====================================================================================