Subroutine waterbalance_init
99 :
100 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
103 : !*-----Run initial run--------------------------------------------------
104 :
105 : subroutine waterbalance_init ()
106 : use vom_vegwat_mod
107 : implicit none
108 :
109 : INTEGER :: jj
110 :
111 : dtsu_count = 0
112 : dtmax_count = 0
113 : zw_ = i_cz
114 : wlayer_ = 0
115 :
116 : ! * Adapt the number of unsaturated layers such that ys = zr
117 :
118 : do while (zw_ .gt. i_zr)
119 : wlayer_ = wlayer_ + 1
120 : zw_ = zw_ - s_delz(wlayer_)
121 : enddo
122 : pcap_(wlayer_+1:s_maxlayer) = 0.d0
123 :
124 : ! * Equilibrium pressure head
125 :
126 : do jj = wlayer_, 1, -1
127 : pcap_(jj) = 0.5d0 * s_delz(jj+1) + 0.5d0 * s_delz(jj) + pcap_(jj+1)
128 : enddo
129 :
130 : ! * equilibrium su above a saturated layer (after eq_sueq1 in Watbal3)
131 :
132 : do jj = 1, s_maxlayer - 1
133 : sueq(jj) = (1.d0 / ((0.5d0 * (s_delz(jj+1) + s_delz(jj)) &
134 : & * s_avg(jj)) ** s_nvg(jj) + 1.d0)) ** c_mvg(jj)
135 : enddo
136 : sueq(s_maxlayer) = (1.d0 / ((0.5d0 * s_delz(jj) * s_avg(jj)) &
137 : & ** s_nvg(jj) + 1.d0)) ** c_mvg(jj)
138 : su__(:) = (1.d0 / ((pcap_(:) * s_avg(:)) ** s_nvg(:) + 1.d0)) ** c_mvg(:)
139 :
140 : ! * unsat. hydrol. cond. as a function of su
141 :
142 : kunsat_(:) = ((-su__(:) ** (1.d0 / c_mvg(:)) + 1.d0) ** c_mvg(:) &
143 : & - 1.d0) ** 2.d0 * s_ksat(:) * SQRT(su__(:))
144 : cH2Ol_s(:) = (-su__(:) * s_thetar(:) + su__(:) * s_thetas(:) &
145 : & + s_thetar(:)) * s_delz(:)
146 :
147 : zwnew = zw_
148 : wlayernew = wlayer_
149 : pcapnew = pcap_
150 : sunew = su__
151 : kunsatnew = kunsat_
152 :
153 : return
154 : end subroutine waterbalance_init