-
Notifications
You must be signed in to change notification settings - Fork 312
/
SoilWaterMovementMod.F90
2232 lines (1890 loc) · 98.5 KB
/
SoilWaterMovementMod.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 SoilWaterMovementMod
#include "shr_assert.h"
!-----------------------------------------------------------------------
! DESCRIPTION
! module contains different subroutines to couple soil and root water interactions
!
! created by Jinyun Tang, Mar 12, 2014
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_sys_mod , only : shr_sys_flush
!
implicit none
private
!
! !PUBLIC MEMBER FUNCTIONS:
public :: SoilWater ! Calculate soil hydrology
public :: init_soilwater_movement
public :: readParams
private :: soilwater_zengdecker2009
private :: soilwater_moisture_form
! private :: soilwater_mixed_form
! private :: soilwater_head_form
private :: compute_hydraulic_properties
private :: compute_moisture_fluxes_and_derivs
private :: compute_RHS_moisture_form
private :: compute_LHS_moisture_form
private :: compute_qcharge
private :: IceImpedance
private :: TridiagonalCol
type, private :: params_type
real(r8) :: e_ice ! Soil ice impedance factor (unitless)
end type params_type
type(params_type), private :: params_inst
!
! The following is only public for the sake of unit testing; it should not be called
! directly by CLM code outside this module
public :: BaseflowSink
public :: use_aquifer_layer
!
! !PUBLIC DATA MEMBERS:
! logical, public :: use_bedrock = .false. ! true => run with spatially variable soil depth
! !PRIVATE DATA MEMBERS:
! Solution method
integer, parameter :: zengdecker_2009 = 0
integer, parameter :: moisture_form = 1
integer, parameter :: mixed_form = 2
integer, parameter :: head_form = 3
! Boundary conditions
integer, parameter :: bc_head = 0
integer, parameter :: bc_flux = 1
integer, parameter :: bc_zero_flux = 2
integer, parameter :: bc_waterTable = 3
integer, parameter :: bc_aquifer = 4
! Soil hydraulic properties
integer, parameter :: soil_hp_clapphornberg_1978=0
integer, parameter :: soil_hp_vanGenuchten_1980=1
real(r8),parameter :: m_to_mm = 1.e3_r8 !convert meters to mm
integer :: soilwater_movement_method ! method for solving richards equation
integer :: upper_boundary_condition ! named variable for the boundary condition
integer :: lower_boundary_condition ! named variable for the boundary condition
! Adaptive time stepping algorithmic control parameters
real(r8) :: dtmin ! minimum time step length (seconds)
real(r8) :: verySmall ! a very small number: used to check for sub step completion
real(r8) :: xTolerUpper ! tolerance to halve length of substep
real(r8) :: xTolerLower ! tolerance to double length of substep
integer :: expensive
integer :: inexpensive
integer :: flux_calculation
character(len=*), parameter, private :: sourcefile = &
__FILE__
!-----------------------------------------------------------------------
contains
!#1
!-----------------------------------------------------------------------
subroutine readParams( ncid )
!
! !USES:
use ncdio_pio, only: file_desc_t
use paramUtilMod, only: readNcdioScalar
!
! !ARGUMENTS:
implicit none
type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id
!
! !LOCAL VARIABLES:
character(len=*), parameter :: subname = 'readParams_SoilWaterMovement'
!--------------------------------------------------------------------
! Soil ice impedance factor (unitless)
call readNcdioScalar(ncid, 'e_ice', subname, params_inst%e_ice)
end subroutine readParams
!#2
!-----------------------------------------------------------------------
subroutine init_soilwater_movement()
!
!DESCRIPTION
!specify method for doing soil&root water interactions
!
! !USES:
use abortutils , only : endrun
use fileutils , only : getavu, relavu
use spmdMod , only : mpicom, masterproc
use shr_mpi_mod , only : shr_mpi_bcast
use clm_varctl , only : iulog, use_bedrock
use controlMod , only : NLFilename
use clm_nlUtilsMod , only : find_nlgroup_name
! !ARGUMENTS:
!------------------------------------------------------------------------------
implicit none
integer :: nu_nml ! unit for namelist file
integer :: nml_error ! namelist i/o error flag
character(*), parameter :: subName = "('init_soilwater_movement')"
!-----------------------------------------------------------------------
! MUST agree with name in namelist and read statement
namelist /soilwater_movement_inparm/ &
soilwater_movement_method, &
upper_boundary_condition, &
lower_boundary_condition, &
dtmin, &
verySmall, &
xTolerUpper, &
xTolerLower, &
expensive, &
inexpensive, &
flux_calculation
! Default values for namelist
soilwater_movement_method = zengdecker_2009
upper_boundary_condition = bc_flux
lower_boundary_condition = bc_aquifer
dtmin=60._r8
verySmall=1.e-8_r8
xTolerUpper=1.e-1_r8
xTolerLower=1.e-2_r8
expensive=42
inexpensive=1
flux_calculation=inexpensive
! Read soilwater_movement namelist
if (masterproc) then
nu_nml = getavu()
open( nu_nml, file=trim(NLFilename), status='old', iostat=nml_error )
call find_nlgroup_name(nu_nml, 'soilwater_movement_inparm', status=nml_error)
if (nml_error == 0) then
read(nu_nml, nml=soilwater_movement_inparm,iostat=nml_error)
if (nml_error /= 0) then
call endrun(subname // ':: ERROR reading soilwater_movement namelist')
end if
else
call endrun(subname // ':: ERROR reading soilwater_movement namelist')
end if
close(nu_nml)
call relavu( nu_nml )
! test for namelist consistency
if((soilwater_movement_method == zengdecker_2009) .and. &
(lower_boundary_condition /= bc_aquifer)) then
call endrun(subname // ':: ERROR inconsistent soilwater_movement namelist: ZD09 must use bc_aquifer lbc')
endif
if((use_bedrock) .and. (lower_boundary_condition /= bc_zero_flux)) then
call endrun(subname // ':: ERROR inconsistent soilwater_movement namelist: use_bedrock requires bc_zero_flux lbc')
endif
endif
call shr_mpi_bcast(soilwater_movement_method, mpicom)
call shr_mpi_bcast(upper_boundary_condition, mpicom)
call shr_mpi_bcast(lower_boundary_condition, mpicom)
call shr_mpi_bcast(dtmin, mpicom)
call shr_mpi_bcast(verySmall, mpicom)
call shr_mpi_bcast(xTolerUpper, mpicom)
call shr_mpi_bcast(xTolerLower, mpicom)
call shr_mpi_bcast(expensive, mpicom)
call shr_mpi_bcast(inexpensive, mpicom)
call shr_mpi_bcast(flux_calculation, mpicom)
if (masterproc) then
write(iulog,*) ' '
write(iulog,*) 'soilwater_movement settings:'
write(iulog,*) ' soilwater_movement_method = ',soilwater_movement_method
write(iulog,*) ' upper_boundary_condition = ',upper_boundary_condition
write(iulog,*) ' lower_boundary_condition = ',lower_boundary_condition
write(iulog,*) ' use_bedrock = ',use_bedrock
write(iulog,*) ' dtmin = ',dtmin
write(iulog,*) ' verySmall = ',verySmall
write(iulog,*) ' xTolerUpper = ',xTolerUpper
write(iulog,*) ' xTolerLower = ',xTolerLower
write(iulog,*) ' expensive = ',expensive
write(iulog,*) ' inexpensive = ',inexpensive
write(iulog,*) ' flux_calculation = ',flux_calculation
endif
end subroutine init_soilwater_movement
!#3
!------------------------------------------------------------------------------
function use_aquifer_layer() result(lres)
!
!DESCRIPTION
! return true if an aquifer layer is used
! otherwise false
implicit none
logical :: lres
if(lower_boundary_condition == bc_aquifer .or. lower_boundary_condition == bc_watertable)then
lres=.true.
else
lres=.false.
endif
return
end function use_aquifer_layer
!#4
!-----------------------------------------------------------------------
subroutine SoilWater(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, &
waterfluxbulk_inst, waterstatebulk_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
!
! DESCRIPTION
! select one subroutine to do the soil and root water coupling
!
!USES
use shr_kind_mod , only : r8 => shr_kind_r8
use clm_varpar , only : nlevsoi
use decompMod , only : bounds_type
use abortutils , only : endrun
use clm_varpar , only : nlevsoi
use SoilHydrologyType , only : soilhydrology_type
use SoilStateType , only : soilstate_type
use TemperatureType , only : temperature_type
use WaterFluxBulkType , only : waterfluxbulk_type
use EnergyFluxType , only : energyflux_type
use WaterStateBulkType , only : waterstatebulk_type
use CanopyStateType , only : canopystate_type
use ColumnType , only : col
use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type
use clm_varcon , only : denh2o, denice
use clm_varctl, only : use_flexibleCN
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds ! bounds
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(:) ! column filter for urban points
type(soilhydrology_type) , intent(inout) :: soilhydrology_inst
type(soilstate_type) , intent(inout) :: soilstate_inst
type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst
type(energyflux_type) , intent(in) :: energyflux_inst
type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst
type(canopystate_type) , intent(inout) :: canopystate_inst
type(temperature_type) , intent(in) :: temperature_inst
class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve
!
! !LOCAL VARIABLES:
character(len=32) :: subname = 'SoilWater' ! subroutine name
real(r8) :: xs(bounds%begc:bounds%endc) !excess soil water above urban ponding limit
real(r8) :: watmin
integer :: fc, c, j
!------------------------------------------------------------------------------
associate( &
wa => waterstatebulk_inst%wa_col , & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm)
dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m)
h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2)
h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2)
h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col & ! Output: [real(r8) (:,:) ] liquid water (kg/m2)
)
select case(soilwater_movement_method)
case (zengdecker_2009)
call soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, &
waterfluxbulk_inst, waterstatebulk_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
case (moisture_form)
call soilwater_moisture_form(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, &
waterfluxbulk_inst, waterstatebulk_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
case (mixed_form)
!!$ call soilwater_mixed_form(bounds, num_hydrologyc, filter_hydrologyc, &
!!$ num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, &
!!$ waterfluxbulk_inst, waterstate_inst, temperature_inst)
case (head_form)
!!$ call soilwater_head_form(bounds, num_hydrologyc, filter_hydrologyc, &
!!$ num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, &
!!$ waterfluxbulk_inst, waterstate_inst, temperature_inst)
case default
call endrun(subname // ':: a SoilWater implementation must be specified!')
end select
if (use_flexibleCN) then
!a work around of the negative liquid water. Jinyun Tang, Jan 14, 2015
watmin = 0.001_r8
do j = 1, nlevsoi-1
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (h2osoi_liq(c,j) < 0._r8) then
xs(c) = watmin - h2osoi_liq(c,j)
else
xs(c) = 0._r8
end if
h2osoi_liq(c,j ) = h2osoi_liq(c,j ) + xs(c)
h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) - xs(c)
end do
end do
j = nlevsoi
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (h2osoi_liq(c,j) < watmin) then
xs(c) = watmin-h2osoi_liq(c,j)
else
xs(c) = 0._r8
end if
h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c)
wa(c) = wa(c) - xs(c)
end do
!update volumetric soil moisture for bgc calculation
do j = 1, nlevsoi
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) &
+ h2osoi_ice(c,j)/(dz(c,j)*denice)
enddo
enddo
end if
end associate
end subroutine SoilWater
!#5
!-----------------------------------------------------------------------
subroutine BaseflowSink(bounds, num_hydrologyc, &
filter_hydrologyc, baseflow_sink, waterfluxbulk_inst, soilstate_inst)
!
! Generic routine to apply baseflow as a sink condition that
! is vertically distributed over the soil column.
!
!USES:
use decompMod , only : bounds_type
use shr_kind_mod , only : r8 => shr_kind_r8
use clm_varpar , only : nlevsoi, max_patch_per_col
use SoilStateType , only : soilstate_type
use WaterFluxBulkType , only : waterfluxbulk_type
use PatchType , only : patch
use ColumnType , only : col
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds ! bounds
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points
real(r8) , intent(out) :: baseflow_sink(bounds%begc:,1:) ! vertically distributed baseflow sink (mm H2O/s) (+ = to rof)
type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst
type(soilstate_type) , intent(inout) :: soilstate_inst
!
! !LOCAL VARIABLES:
integer :: p,c,fc,j ! do loop indices
integer :: pi ! patch index
real(r8) :: temp(bounds%begc:bounds%endc) ! accumulator for rootr weighting
!-----------------------------------------------------------------------
! Enforce expected array sizes
SHR_ASSERT_ALL_FL((ubound(baseflow_sink) == (/bounds%endc, nlevsoi/)), sourcefile, __LINE__)
!this is just a placeholder for now
baseflow_sink = 0.
end subroutine BaseflowSink
!#6
!-----------------------------------------------------------------------
subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, &
waterfluxbulk_inst, waterstatebulk_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
!
! !DESCRIPTION:
! Soil hydrology
! Soil moisture is predicted from a 10-layer model (as with soil
! temperature), in which the vertical soil moisture transport is governed
! by infiltration, runoff, gradient diffusion, gravity, and root
! extraction through canopy transpiration. The net water applied to the
! surface layer is the snowmelt plus precipitation plus the throughfall
! of canopy dew minus surface runoff and evaporation.
! CLM3.5 uses a zero-flow bottom boundary condition.
!
! The vertical water flow in an unsaturated porous media is described by
! Darcy's law, and the hydraulic conductivity and the soil negative
! potential vary with soil water content and soil texture based on the work
! of Clapp and Hornberger (1978) and Cosby et al. (1984). The equation is
! integrated over the layer thickness, in which the time rate of change in
! water mass must equal the net flow across the bounding interface, plus the
! rate of internal source or sink. The terms of water flow across the layer
! interfaces are linearly expanded by using first-order Taylor expansion.
! The equations result in a tridiagonal system equation.
!
! Note: length units here are all millimeter
! (in temperature subroutine uses same soil layer
! structure required but lengths are m)
!
! Richards equation:
!
! d wat d d wat d psi
! ----- = - -- [ k(----- ----- - 1) ] + S
! dt dz dz d wat
!
! where: wat = volume of water per volume of soil (mm**3/mm**3)
! psi = soil matrix potential (mm)
! dt = time step (s)
! z = depth (mm)
! dz = thickness (mm)
! qin = inflow at top (mm h2o /s)
! qout= outflow at bottom (mm h2o /s)
! s = source/sink flux (mm h2o /s)
! k = hydraulic conductivity (mm h2o /s)
!
! d qin d qin
! qin[n+1] = qin[n] + -------- d wat(j-1) + --------- d wat(j)
! d wat(j-1) d wat(j)
! ==================|=================
! < qin
!
! d wat(j)/dt * dz = qin[n+1] - qout[n+1] + S(j)
!
! > qout
! ==================|=================
! d qout d qout
! qout[n+1] = qout[n] + --------- d wat(j) + --------- d wat(j+1)
! d wat(j) d wat(j+1)
!
!
! Solution: linearize k and psi about d wat and use tridiagonal
! system of equations to solve for d wat,
! where for layer j
!
!
! r_j = a_j [d wat_j-1] + b_j [d wat_j] + c_j [d wat_j+1]
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_const_mod , only : SHR_CONST_TKFRZ, SHR_CONST_LATICE, SHR_CONST_G
use decompMod , only : bounds_type
use clm_varcon , only : grav,hfus,tfrz
use clm_varcon , only : denh2o, denice
use clm_varpar , only : nlevsoi, max_patch_per_col, nlevgrnd
use clm_time_manager , only : get_step_size_real, get_nstep
use column_varcon , only : icol_roof, icol_road_imperv
use clm_varctl , only : use_flexibleCN, use_hydrstress
use TridiagonalMod , only : Tridiagonal
use abortutils , only : endrun
use SoilStateType , only : soilstate_type
use SoilHydrologyType , only : soilhydrology_type
use TemperatureType , only : temperature_type
use WaterFluxBulkType , only : waterfluxbulk_type
use EnergyFluxType , only : energyflux_type
use WaterStateBulkType , only : waterstatebulk_type
use CanopyStateType , only : canopystate_type
use SoilWaterRetentionCurveMod , only : soil_water_retention_curve_type
use PatchType , only : patch
use ColumnType , only : col
use clm_varctl , only : iulog
use SoilWaterPlantSinkMod , only : COmpute_EffecRootFrac_And_VertTranSink
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds ! bounds
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(:) ! column filter for urban points
type(soilhydrology_type), intent(inout) :: soilhydrology_inst
type(soilstate_type) , intent(inout) :: soilstate_inst
type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst
type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst
type(canopystate_type) , intent(inout) :: canopystate_inst
type(temperature_type) , intent(in) :: temperature_inst
type(energyflux_type) , intent(in) :: energyflux_inst
class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve
!
! !LOCAL VARIABLES:
character(len=32) :: subname = 'soilwater_zengdecker2009' ! subroutine name
integer :: p,c,fc,j ! do loop indices
integer :: jtop(bounds%begc:bounds%endc) ! top level at each column
real(r8) :: dtime ! land model time step (sec)
real(r8) :: hk(bounds%begc:bounds%endc,1:nlevsoi) ! hydraulic conductivity [mm h2o/s]
real(r8) :: dhkdw(bounds%begc:bounds%endc,1:nlevsoi) ! d(hk)/d(vol_liq)
real(r8) :: amx(bounds%begc:bounds%endc,1:nlevsoi+1) ! "a" left off diagonal of tridiagonal matrix
real(r8) :: bmx(bounds%begc:bounds%endc,1:nlevsoi+1) ! "b" diagonal column for tridiagonal matrix
real(r8) :: cmx(bounds%begc:bounds%endc,1:nlevsoi+1) ! "c" right off diagonal tridiagonal matrix
real(r8) :: rmx(bounds%begc:bounds%endc,1:nlevsoi+1) ! "r" forcing term of tridiagonal matrix
real(r8) :: zmm(bounds%begc:bounds%endc,1:nlevsoi+1) ! layer depth [mm]
real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevsoi+1) ! layer thickness [mm]
real(r8) :: den ! used in calculating qin, qout
real(r8) :: dqidw0(bounds%begc:bounds%endc,1:nlevsoi+1) ! d(qin)/d(vol_liq(i-1))
real(r8) :: dqidw1(bounds%begc:bounds%endc,1:nlevsoi+1) ! d(qin)/d(vol_liq(i))
real(r8) :: dqodw1(bounds%begc:bounds%endc,1:nlevsoi+1) ! d(qout)/d(vol_liq(i))
real(r8) :: dqodw2(bounds%begc:bounds%endc,1:nlevsoi+1) ! d(qout)/d(vol_liq(i+1))
real(r8) :: dsmpdw(bounds%begc:bounds%endc,1:nlevsoi+1) ! d(smp)/d(vol_liq)
real(r8) :: num ! used in calculating qin, qout
real(r8) :: qin(bounds%begc:bounds%endc,1:nlevsoi+1) ! flux of water into soil layer [mm h2o/s]
real(r8) :: qout(bounds%begc:bounds%endc,1:nlevsoi+1) ! flux of water out of soil layer [mm h2o/s]
real(r8) :: s_node ! soil wetness
real(r8) :: s1 ! "s" at interface of layer
real(r8) :: s2 ! k*s**(2b+2)
real(r8) :: smp(bounds%begc:bounds%endc,1:nlevsoi) ! soil matrix potential [mm]
real(r8) :: sdamp ! extrapolates soiwat dependence of evaporation
integer :: pi ! patch index
real(r8) :: temp(bounds%begc:bounds%endc) ! accumulator for rootr weighting
integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-)
real(r8) :: smp1,dsmpdw1,wh,wh_zwt,ka
real(r8) :: dwat2(bounds%begc:bounds%endc,1:nlevsoi+1)
real(r8) :: dzq ! used in calculating qin, qout (difference in equilbirium matric potential)
real(r8) :: zimm(bounds%begc:bounds%endc,0:nlevsoi) ! layer interface depth [mm]
real(r8) :: zq(bounds%begc:bounds%endc,1:nlevsoi+1) ! equilibrium matric potential for each layer [mm]
real(r8) :: vol_eq(bounds%begc:bounds%endc,1:nlevsoi+1) ! equilibrium volumetric water content
real(r8) :: tempi ! temp variable for calculating vol_eq
real(r8) :: temp0 ! temp variable for calculating vol_eq
real(r8) :: voleq1 ! temp variable for calculating vol_eq
real(r8) :: zwtmm(bounds%begc:bounds%endc) ! water table depth [mm]
real(r8) :: imped(bounds%begc:bounds%endc,1:nlevsoi)
real(r8) :: vol_ice(bounds%begc:bounds%endc,1:nlevsoi)
real(r8) :: z_mid
real(r8) :: vwc_zwt(bounds%begc:bounds%endc)
real(r8) :: vwc_liq(bounds%begc:bounds%endc,1:nlevsoi+1) ! liquid volumetric water content
real(r8) :: smp_grad(bounds%begc:bounds%endc,1:nlevsoi+1)
real(r8) :: dsmpds !temporary variable
real(r8) :: dhkds !temporary variable
real(r8) :: hktmp !temporary variable
integer :: nstep
!-----------------------------------------------------------------------
associate(&
z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m)
zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m)
dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m)
origflag => soilhydrology_inst%origflag , & ! Input: constant
qcharge => soilhydrology_inst%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s)
zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m)
fracice => soilhydrology_inst%fracice_col , & ! Input: [real(r8) (:,:) ] fractional impermeability (-)
icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice
hkdepth => soilhydrology_inst%hkdepth_col , & ! Input: [real(r8) (:) ] decay factor (m)
smpmin => soilstate_inst%smpmin_col , & ! Input: [real(r8) (:) ] restriction for min of soil potential (mm)
watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s)
bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm)
eff_porosity => soilstate_inst%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice
smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm]
hk_l => soilstate_inst%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s)
h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice water (kg/m2)
h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2)
h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]
qflx_deficit => waterfluxbulk_inst%qflx_deficit_col , & ! Input: [real(r8) (:) ] water deficit to keep non-negative liquid water content
qflx_infl => waterfluxbulk_inst%qflx_infl_col , & ! Input: [real(r8) (:) ] infiltration (mm H2O /s)
qflx_rootsoi_col => waterfluxbulk_inst%qflx_rootsoi_col , & ! Output: [real(r8) (:,:) ] vegetation/soil water exchange (mm H2O/s) (+ = to atm)
qflx_drain_vr_col => waterfluxbulk_inst%qflx_drain_vr_col , & ! Input: [real(r8) (:,:) ] drainage from soil layers due to plant induced
t_soisno => temperature_inst%t_soisno_col & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin)
)
! Get time step
nstep = get_nstep()
dtime = get_step_size_real()
! Because the depths in this routine are in mm, use local
! variable arrays instead of pointers
do j = 1, nlevsoi
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
zmm(c,j) = z(c,j)*1.e3_r8
dzmm(c,j) = dz(c,j)*1.e3_r8
zimm(c,j) = zi(c,j)*1.e3_r8
! calculate icefrac up here
vol_ice(c,j) = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))
icefrac(c,j) = min(1._r8,vol_ice(c,j)/watsat(c,j))
vwc_liq(c,j) = max(h2osoi_liq(c,j),1.0e-6_r8)/(dz(c,j)*denh2o)
end do
end do
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
zimm(c,0) = 0.0_r8
zwtmm(c) = zwt(c)*1.e3_r8
end do
! compute jwt index
! The layer index of the first unsaturated layer, i.e., the layer right above
! the water table
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
jwt(c) = nlevsoi
! allow jwt to equal zero when zwt is in top layer
do j = 1,nlevsoi
if(zwt(c) <= zi(c,j)) then
jwt(c) = j-1
exit
end if
enddo
! compute vwc at water table depth (mainly for case when t < tfrz)
! this will only be used when zwt is below the soil column
vwc_zwt(c) = watsat(c,nlevsoi)
if(t_soisno(c,jwt(c)+1) < tfrz) then
vwc_zwt(c) = vwc_liq(c,nlevsoi)
do j = nlevsoi,nlevgrnd
if(zwt(c) <= zi(c,j)) then
smp1 = hfus*(tfrz-t_soisno(c,j))/(grav*t_soisno(c,j)) * 1000._r8 !(mm)
!smp1 = max(0._r8,smp1)
smp1 = max(sucsat(c,nlevsoi),smp1)
vwc_zwt(c) = watsat(c,nlevsoi)*(smp1/sucsat(c,nlevsoi))**(-1._r8/bsw(c,nlevsoi))
! for temperatures close to tfrz, limit vwc to total water content
vwc_zwt(c) = min(vwc_zwt(c), 0.5*(watsat(c,nlevsoi) + h2osoi_vol(c,nlevsoi)) )
exit
endif
enddo
endif
end do
! calculate the equilibrium water content based on the water table depth
do j=1,nlevsoi
do fc=1, num_hydrologyc
c = filter_hydrologyc(fc)
if ((zwtmm(c) <= zimm(c,j-1))) then
vol_eq(c,j) = watsat(c,j)
! use the weighted average from the saturated part (depth > wtd) and the equilibrium solution for the
! rest of the layer, the equilibrium solution is based on Clapp-Hornberg parameterization
! and no extension to full range swrc is needed
else if ((zwtmm(c) .lt. zimm(c,j)) .and. (zwtmm(c) .gt. zimm(c,j-1))) then
tempi = 1.0_r8
temp0 = (((sucsat(c,j)+zwtmm(c)-zimm(c,j-1))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
voleq1 = -sucsat(c,j)*watsat(c,j)/(1._r8-1._r8/bsw(c,j))/(zwtmm(c)-zimm(c,j-1))*(tempi-temp0)
vol_eq(c,j) = (voleq1*(zwtmm(c)-zimm(c,j-1)) + watsat(c,j)*(zimm(c,j)-zwtmm(c)))/(zimm(c,j)-zimm(c,j-1))
vol_eq(c,j) = min(watsat(c,j),vol_eq(c,j))
vol_eq(c,j) = max(vol_eq(c,j),0.0_r8)
else
tempi = (((sucsat(c,j)+zwtmm(c)-zimm(c,j))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
temp0 = (((sucsat(c,j)+zwtmm(c)-zimm(c,j-1))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
vol_eq(c,j) = -sucsat(c,j)*watsat(c,j)/(1._r8-1._r8/bsw(c,j))/(zimm(c,j)-zimm(c,j-1))*(tempi-temp0)
vol_eq(c,j) = max(vol_eq(c,j),0.0_r8)
vol_eq(c,j) = min(watsat(c,j),vol_eq(c,j))
endif
zq(c,j) = -sucsat(c,j)*(max(vol_eq(c,j)/watsat(c,j),0.01_r8))**(-bsw(c,j))
zq(c,j) = max(smpmin(c), zq(c,j))
end do
end do
! If water table is below soil column calculate zq for the 11th layer
j = nlevsoi
do fc=1, num_hydrologyc
c = filter_hydrologyc(fc)
if(jwt(c) == nlevsoi) then
tempi = 1._r8
temp0 = (((sucsat(c,j)+zwtmm(c)-zimm(c,j))/sucsat(c,j)))**(1._r8-1._r8/bsw(c,j))
vol_eq(c,j+1) = -sucsat(c,j)*watsat(c,j)/(1._r8-1._r8/bsw(c,j))/(zwtmm(c)-zimm(c,j))*(tempi-temp0)
vol_eq(c,j+1) = max(vol_eq(c,j+1),0.0_r8)
vol_eq(c,j+1) = min(watsat(c,j),vol_eq(c,j+1))
zq(c,j+1) = -sucsat(c,j)*(max(vol_eq(c,j+1)/watsat(c,j),0.01_r8))**(-bsw(c,j))
zq(c,j+1) = max(smpmin(c), zq(c,j+1))
end if
end do
! Hydraulic conductivity and soil matric potential and their derivatives
sdamp = 0._r8
do j = 1, nlevsoi
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
! compute hydraulic conductivity based on liquid water content only
if (origflag == 1) then
s1 = 0.5_r8*(h2osoi_vol(c,j) + h2osoi_vol(c,min(nlevsoi, j+1))) / &
(0.5_r8*(watsat(c,j)+watsat(c,min(nlevsoi, j+1))))
else
s1 = 0.5_r8*(vwc_liq(c,j) + vwc_liq(c,min(nlevsoi, j+1))) / &
(0.5_r8*(watsat(c,j)+watsat(c,min(nlevsoi, j+1))))
endif
s1 = min(1._r8, s1)
s2 = hksat(c,j)*s1**(2._r8*bsw(c,j)+2._r8)
! replace fracice with impedance factor, as in zhao 97,99
if (origflag == 1) then
imped(c,j)=(1._r8-0.5_r8*(fracice(c,j)+fracice(c,min(nlevsoi, j+1))))
else
imped(c,j)=10._r8**(-params_inst%e_ice*(0.5_r8*(icefrac(c,j)+icefrac(c,min(nlevsoi, j+1)))))
endif
hk(c,j) = imped(c,j)*s1*s2
dhkdw(c,j) = imped(c,j)*(2._r8*bsw(c,j)+3._r8)*s2* &
(1._r8/(watsat(c,j)+watsat(c,min(nlevsoi, j+1))))
!compute un-restricted hydraulic conductivity
!call soil_water_retention_curve%soil_hk(hksat(c,j), imped(c,j), s1, bsw(c,j), hktmp, dhkds)
!if(hktmp/=hk(c,j))write(10,*)'diff',hktmp,hk(c,j)
! call endrun('bad in hk')
!endif
!apply ice impedance
!hk(c,j) = imped(c,j)*hk(c,j)
!dhkdw(c,j) = imped(c,j) * dhkds * (1._r8/(watsat(c,j)+watsat(c,min(nlevsoi, j+1))))
! compute matric potential and derivative based on liquid water content only
if (origflag == 1) then
s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
else
s_node = max(vwc_liq(c,j)/watsat(c,j), 0.01_r8)
endif
s_node = min(1.0_r8, s_node)
!call soil_water_retention_curve%soil_suction(sucsat(c,j), s_node, bsw(c,j), smp(c,j), dsmpds)
smp(c,j) = -sucsat(c,j)*s_node**(-bsw(c,j))
smp(c,j) = max(smpmin(c), smp(c,j))
!do not turn on the line below, which will cause bit to bit error, jyt, 2014 Mar 6
!dsmpdw(c,j) = dsmpds/watsat(c,j)
if (origflag == 1) then
dsmpdw(c,j) = -bsw(c,j)*smp(c,j)/(s_node*watsat(c,j))
else
dsmpdw(c,j) = -bsw(c,j)*smp(c,j)/vwc_liq(c,j)
endif
smp_l(c,j) = smp(c,j)
hk_l(c,j) = hk(c,j)
end do
end do
! aquifer (11th) layer
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
zmm(c,nlevsoi+1) = 0.5*(1.e3_r8*zwt(c) + zmm(c,nlevsoi))
if(jwt(c) < nlevsoi) then
dzmm(c,nlevsoi+1) = dzmm(c,nlevsoi)
else
dzmm(c,nlevsoi+1) = (1.e3_r8*zwt(c) - zmm(c,nlevsoi))
end if
end do
! Set up r, a, b, and c vectors for tridiagonal solution
! Node j=1 (top)
j = 1
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
qin(c,j) = qflx_infl(c)
den = (zmm(c,j+1)-zmm(c,j))
dzq = (zq(c,j+1)-zq(c,j))
num = (smp(c,j+1)-smp(c,j)) - dzq
qout(c,j) = -hk(c,j)*num/den
dqodw1(c,j) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqodw2(c,j) = -( hk(c,j)*dsmpdw(c,j+1) + num*dhkdw(c,j))/den
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_rootsoi_col(c,j)
amx(c,j) = 0._r8
bmx(c,j) = dzmm(c,j)*(sdamp+1._r8/dtime) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
end do
! Nodes j=2 to j=nlevsoi-1
do j = 2, nlevsoi - 1
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
den = (zmm(c,j) - zmm(c,j-1))
dzq = (zq(c,j)-zq(c,j-1))
num = (smp(c,j)-smp(c,j-1)) - dzq
qin(c,j) = -hk(c,j-1)*num/den
dqidw0(c,j) = -(-hk(c,j-1)*dsmpdw(c,j-1) + num*dhkdw(c,j-1))/den
dqidw1(c,j) = -( hk(c,j-1)*dsmpdw(c,j) + num*dhkdw(c,j-1))/den
den = (zmm(c,j+1)-zmm(c,j))
dzq = (zq(c,j+1)-zq(c,j))
num = (smp(c,j+1)-smp(c,j)) - dzq
qout(c,j) = -hk(c,j)*num/den
dqodw1(c,j) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqodw2(c,j) = -( hk(c,j)*dsmpdw(c,j+1) + num*dhkdw(c,j))/den
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_rootsoi_col(c,j)
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
end do
end do
! Node j=nlevsoi (bottom)
j = nlevsoi
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if(j > jwt(c)) then !water table is in soil column
den = (zmm(c,j) - zmm(c,j-1))
dzq = (zq(c,j)-zq(c,j-1))
num = (smp(c,j)-smp(c,j-1)) - dzq
qin(c,j) = -hk(c,j-1)*num/den
dqidw0(c,j) = -(-hk(c,j-1)*dsmpdw(c,j-1) + num*dhkdw(c,j-1))/den
dqidw1(c,j) = -( hk(c,j-1)*dsmpdw(c,j) + num*dhkdw(c,j-1))/den
qout(c,j) = 0._r8
dqodw1(c,j) = 0._r8
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_rootsoi_col(c,j)
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = 0._r8
! next set up aquifer layer; hydrologically inactive
rmx(c,j+1) = 0._r8
amx(c,j+1) = 0._r8
bmx(c,j+1) = dzmm(c,j+1)/dtime
cmx(c,j+1) = 0._r8
else ! water table is below soil column
! compute aquifer soil moisture as average of layer 10 and saturation
if(origflag == 1) then
s_node = max(0.5*(1.0_r8+h2osoi_vol(c,j)/watsat(c,j)), 0.01_r8)
else
s_node = max(0.5*((vwc_zwt(c)+vwc_liq(c,j))/watsat(c,j)), 0.01_r8)
endif
s_node = min(1.0_r8, s_node)
! compute smp for aquifer layer
!call soil_water_retention_curve%soil_suction(sucsat(c,j), s_node, bsw(c,j), smp1, dsmpds)
smp1 = -sucsat(c,j)*s_node**(-bsw(c,j))
smp1 = max(smpmin(c), smp1)
! compute dsmpdw for aquifer layer
!dsmpdw1 = dsmpds/watsat(c,j)
dsmpdw1 = -bsw(c,j)*smp1/(s_node*watsat(c,j))
! first set up bottom layer of soil column
den = (zmm(c,j) - zmm(c,j-1))
dzq = (zq(c,j)-zq(c,j-1))
num = (smp(c,j)-smp(c,j-1)) - dzq
qin(c,j) = -hk(c,j-1)*num/den
dqidw0(c,j) = -(-hk(c,j-1)*dsmpdw(c,j-1) + num*dhkdw(c,j-1))/den
dqidw1(c,j) = -( hk(c,j-1)*dsmpdw(c,j) + num*dhkdw(c,j-1))/den
den = (zmm(c,j+1)-zmm(c,j))
dzq = (zq(c,j+1)-zq(c,j))
num = (smp1-smp(c,j)) - dzq
qout(c,j) = -hk(c,j)*num/den
dqodw1(c,j) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqodw2(c,j) = -( hk(c,j)*dsmpdw1 + num*dhkdw(c,j))/den
rmx(c,j) = qin(c,j) - qout(c,j) - qflx_rootsoi_col(c,j)
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
! next set up aquifer layer; den/num unchanged, qin=qout
qin(c,j+1) = qout(c,j)
dqidw0(c,j+1) = -(-hk(c,j)*dsmpdw(c,j) + num*dhkdw(c,j))/den
dqidw1(c,j+1) = -( hk(c,j)*dsmpdw1 + num*dhkdw(c,j))/den
qout(c,j+1) = 0._r8 ! zero-flow bottom boundary condition
dqodw1(c,j+1) = 0._r8 ! zero-flow bottom boundary condition
rmx(c,j+1) = qin(c,j+1) - qout(c,j+1)
amx(c,j+1) = -dqidw0(c,j+1)
bmx(c,j+1) = dzmm(c,j+1)/dtime - dqidw1(c,j+1) + dqodw1(c,j+1)
cmx(c,j+1) = 0._r8
endif
end do
! Solve for dwat
jtop(bounds%begc : bounds%endc) = 1
call Tridiagonal(bounds, 1, nlevsoi+1, &
jtop(bounds%begc:bounds%endc), &
num_hydrologyc, filter_hydrologyc, &
amx(bounds%begc:bounds%endc, :), &
bmx(bounds%begc:bounds%endc, :), &
cmx(bounds%begc:bounds%endc, :), &
rmx(bounds%begc:bounds%endc, :), &
dwat2(bounds%begc:bounds%endc, :) )
! Renew the mass of liquid water
! also compute qcharge from dwat in aquifer layer
! update in drainage for case jwt < nlevsoi
do fc = 1,num_hydrologyc
c = filter_hydrologyc(fc)
do j = 1, nlevsoi
h2osoi_liq(c,j) = h2osoi_liq(c,j) + dwat2(c,j)*dzmm(c,j)
end do
! calculate qcharge for case jwt < nlevsoi
if(jwt(c) < nlevsoi) then
wh_zwt = 0._r8 !since wh_zwt = -sucsat - zq_zwt, where zq_zwt = -sucsat
! Recharge rate qcharge to groundwater (positive to aquifer)
s_node = max(h2osoi_vol(c,jwt(c)+1)/watsat(c,jwt(c)+1), 0.01_r8)
s1 = min(1._r8, s_node)
!scs: this is the expression for unsaturated hk
ka = imped(c,jwt(c)+1)*hksat(c,jwt(c)+1) &
*s1**(2._r8*bsw(c,jwt(c)+1)+3._r8)
!compute unsaturated hk, this shall be tested later, because it
!is not bit for bit
!call soil_water_retention_curve%soil_hk(hksat(c,jwt(c)+1), s1, bsw(c,jwt(c)+1), ka)
!apply ice impedance
!ka = imped(c,jwt(c)+1) * ka
! Recharge rate qcharge to groundwater (positive to aquifer)
smp1 = max(smpmin(c), smp(c,max(1,jwt(c))))
wh = smp1 - zq(c,max(1,jwt(c)))
!scs: original formulation
if(jwt(c) == 0) then
qcharge(c) = -ka * (wh_zwt-wh) /((zwt(c)+1.e-3)*1000._r8)
else
! qcharge(c) = -ka * (wh_zwt-wh)/((zwt(c)-z(c,jwt(c)))*1000._r8)
!scs: 1/2, assuming flux is at zwt interface, saturation deeper than zwt
qcharge(c) = -ka * (wh_zwt-wh)/((zwt(c)-z(c,jwt(c)))*1000._r8*2.0)
endif
! To limit qcharge (for the first several timesteps)
qcharge(c) = max(-10.0_r8/dtime,qcharge(c))
qcharge(c) = min( 10.0_r8/dtime,qcharge(c))
else
! if water table is below soil column, compute qcharge from dwat2(11)
qcharge(c) = dwat2(c,nlevsoi+1)*dzmm(c,nlevsoi+1)/dtime
endif
end do
! compute the water deficit and reset negative liquid water content
! Jinyun Tang
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
qflx_deficit(c) = 0._r8
do j = 1, nlevsoi
if(h2osoi_liq(c,j)<0._r8)then
qflx_deficit(c) = qflx_deficit(c) - h2osoi_liq(c,j)
endif
enddo
enddo
end associate
end subroutine soilwater_zengdecker2009
!#7
!-----------------------------------------------------------------------
subroutine soilwater_moisture_form(bounds, num_hydrologyc, &
filter_hydrologyc, num_urbanc, filter_urbanc, soilhydrology_inst, &
soilstate_inst, waterfluxbulk_inst, waterstatebulk_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
!
! !DESCRIPTION:
! Soil hydrology
! Soil moisture is predicted from a n-layer model (as with soil
! temperature), in which the vertical soil moisture transport is governed