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