! Original code by P.J. Ross 2001-2007 See: Ross, P.J. (2003) ! Modeling soil water and solute transport - fast, simplified numerical solutions. Agron. J. 95:1352-1361. ! ! Modified by V. Haverd 2008: matrix expanded by factor of 2 to allow solution of coupled heat/moisture fluxes ! Explicit calculation of heat/moisture fluxes at surface ! Switchable option for litter ! Isotope subroutine (V. Haverd and M. Cuntz, September 2008) ! ! Frozen Soil (V. Haverd, M. Cuntz Sep 2010 - Jan 2011) ! Pond lumped with top soil layer (V. Haverd Jan 2011) ! Snow included in layers -2:0 (V. Haverd Jan 2011) ! Convergence test for dthetaldT (V.Haverd Feb 2011) ! Include heat advection by liquid water flux (V.Haverd Feb 2011) ! Code cleaned and isotope routine updated (V. Haverd, M. Cuntz Aug-Sep 2012) ! lots of changes, especially relating to ponding and freezing together ! Outstanding problems with surface runoff and advection (set h0max large and turn off advection for now) ! ! Sep 17 2012 (V.Haverd) ! remove separate solution for maxpond (h0 > h0max) ! Instead, convert excess pond to runoff at t=tfin, and correct energy stores for loss of pond ! Advection and surface runoff now functioning ! ! Dec 31 2012 (V.Haverd) ! bug fixes to improve energy conservation associated with change in ice, ! especially correction at time of pond disappearance (removal of negative pond) ! ! Jan 5 2013 (V.Haverd) ! Remove iteration over frozen soil within "(while iok==0)" loop ! For frozen soil, at the time of updating variables, evaluate new T and ice content, ! consistent with J0 + deltaJ, where deltaJ is the change in energy obtained from the matrix solution (=LHS_h*dt) ! ! Jan 7-8 2013 (V.Haverd) ! Calls to soil_snow albedo and density ! Improve iteration convergence at time of updating new T and ice content in frozen soil ! Bug fix relating to moisture conservation in snow ! ! Jan 14-15 2013 (V. Haverd) ! Modified definition of Sliq for frozen soil (in hyofS): now defined relative to (theta_sat - theta_r - theta_ice) ! Allow for snow-pack initalisation in the absence of pond ! ! Feb 7 2013 (V.Haverd) ! Revised formulation of dphidT for frozen soil (see hyofS): fixes a few remaining -ve T spikes in frozen soil ! ! Feb 9 2013 (V.Haverd) ! Moved snow layer adjustments to subroutine ! ! Feb 24 2013 (V. Haverd) ! Snow structure expanded to accommodate 3 layers ! ! March 30 (V. Haverd) ! ! Jan 7 2014 (V. Haverd) ! Adjust snow surface bcs such that snow surface temperature can't be more than 0 ! ! Jan 8 2014 (V Haverd) ! enable 2nd snow layer so that if snowpack exceeds 3 cm , the excess goes into the layer below. ! ! March 2014 (V. Haverd) ! Extract SEB calcs to subroutine in utils, and trial Force-Restore method ! ! October 2016 (S Wales) ! Restructured code of solve into several subroutines MODULE sli_solve USE cable_def_types_mod, ONLY: r_2, i_d USE sli_numbers, ONLY: & experiment, & zero, one, two, half, thousand, e3, e5, & ! numbers Tzero, rlambda, lambdaf, lambdas, Dva, rhocp, rhow, gf, hmin, & ! parameters csice , cswat, & snmin, nsnow_max, fsnowliq_max, & params, vars_aquifer, vars_met, vars, vars_snow, solve_type, & ! types dSfac, h0min, Smax, h0max, dSmax, dSmaxr, dtmax, dSmax, dSmaxr, & ! numerical limits dtmax, dtmin, dTsoilmax, dTLmax, nsteps_ice_max, nsteps_max,tol_dthetaldT, & hbot, botbc, & Mw, Rgas, & ! boundary condition freezefac USE sli_utils, ONLY: & x, Sofh, hyofh, hyofS, litter_props, massman_sparse, tri, & aquifer_props, Tfrz, thetalmax, Tthetalmax, dthetalmaxdTh, & getfluxes_vp, getheatfluxes, flux, sol, & csat, slope_csat, potential_evap, tri, setsol, zerovars, & esat_ice, slope_esat_ice, Tfrozen, rtbis_Tfrozen, GTFrozen, & JSoilLayer, esat, forcerestore, SEB USE cable_IO_vars_module, ONLY: logn IMPLICIT NONE PRIVATE PUBLIC :: solve ! solution routine INTEGER(i_d), DIMENSION(:), ALLOCATABLE :: nless, n_noconverge ! global counters ! Definitions of public entities and private parameters (see above for default ! values): ! botbc - bottom boundary condn for water; "constant head", "free drainage", ! "seepage", or "zero flux". Constant head means that matric head h ! is specified. Free drainage means zero gradient of matric head, ! i.e. unit hydraulic gradient. Seepage means zero flux when the ! matric head is below zero and an upper limit of zero for the head. ! h0max - max pond depth allowed before runoff. ! hbot - matric head at bottom of profile when botbc set to "constant head". ! dSmax - max change in S (the "effective saturation") of any unsaturated ! layer to aim for each time step; controls time step size. ! dSmaxr - maximum negative relative change in S each time step. This ! parameter helps avoid very small or negative S. ! dtmax - max time step allowed. ! dsmmax - max solute change per time step (see dSmax); user should set this ! according to solute units used. Units for different solutes can be ! scaled by the user (e.g. to an expected max of around 1.0). ! dSfac - a change in S of up to dSfac*dSmax is accepted. ! dpmaxr - relative change in matric flux potential (MFP) phi that is ! accepted for convergence when finding head h at soil interfaces. ! h0min - min (negative) value for surface pond when it empties. ! Smax - max value for layer saturation to allow some overshoot. ! solve - sub to call to solve RE ! CONTAINS SUBROUTINE get_fluxes_and_derivs( & irec, mp, qprec, qprec_snow, n, dx, h0, & Tsoil, S, & qh, vmet, vlit, vsnow, var, T0, Tsurface, & dxL, SL, Tl, & plit, par, & qali, & qvh, & qevap, & again, getq0,getqn,init, ns, nsat, & nsatlast, & phip, qpme, hint, phimin, & qexd, & q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, qadv, qadvya, qadvyb, & qadvTa, qadvTb, vtmp, qliq, qv, qvT, qlya, qlyb, & qvya, qvyb, qlTb, qvTa, qvTb, & vtop, vbot, & lE0, G0, & Epot, surface_case, & iflux, litter, j, kk, & advection, dTqwdTa, dTqwdTb, Tqw, keff, & hice & ) IMPLICIT NONE INTEGER(i_d) :: irec, mp REAL(r_2), DIMENSION(1:mp) :: qprec REAL(r_2), DIMENSION(1:mp) :: qprec_snow INTEGER(i_d) :: n REAL(r_2), DIMENSION(1:n) :: dx REAL(r_2), DIMENSION(1:mp) :: h0 REAL(r_2), DIMENSION(1:n) :: Tsoil, S REAL(r_2), DIMENSION(-nsnow_max:n) :: qh TYPE(vars_met), DIMENSION(1:mp) :: vmet TYPE(vars), DIMENSION(1:mp) :: vlit TYPE(vars_snow), DIMENSION(1:mp) :: vsnow TYPE(vars), DIMENSION(1:n) :: var REAL(r_2), DIMENSION(1:mp) :: T0, Tsurface REAL(r_2), DIMENSION(1:mp) :: dxL REAL(r_2), DIMENSION(1:mp) :: SL, Tl TYPE(params), DIMENSION(1:mp) :: plit TYPE(params), DIMENSION(1:n) :: par REAL(r_2), DIMENSION(1:mp), OPTIONAL :: qali REAL(r_2), DIMENSION(-nsnow_max:n) :: qvh REAL(r_2), DIMENSION(1:mp) :: qevap LOGICAL, DIMENSION(1:mp) :: again, getq0,getqn,init INTEGER(i_d), DIMENSION(1:mp) :: ns, nsat, nsatlast REAL(r_2), DIMENSION(1:mp) :: phip REAL(r_2), DIMENSION(1:mp) :: qpme REAL(r_2), DIMENSION(1:n) :: hint, phimin, qexd REAL(r_2), DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb REAL(r_2), DIMENSION(-nsnow_max:n) :: qadv, qadvya, qadvyb, qadvTa, qadvTb TYPE(vars) :: vtmp REAL(r_2), DIMENSION(-nsnow_max:n) :: qliq, qv, qvT, qlya, qlyb, qvya, qvyb, qlTb, qvTa, qvTb TYPE(vars), DIMENSION(1:mp) :: vtop, vbot REAL(r_2), DIMENSION(1:mp) :: lE0, G0, Epot INTEGER(i_d), DIMENSION(1:mp) :: surface_case INTEGER(i_d), DIMENSION(1:mp) :: iflux LOGICAL :: litter INTEGER(i_d) :: j, kk INTEGER(i_d) :: advection ! switches REAL(r_2) :: dTqwdTa, dTqwdTb, Tqw, keff REAL(r_2), DIMENSION(1:mp) :: hice !----- get fluxes and derivs ! get surface condition ! ns==1 if no pond or full pond, i.e. do not solve for change in pond height ! ns==0 then change for pond height IF ((var(1)%phi <= phip(kk) .AND. h0(kk) <= zero .AND. nsat(kk) < n) .OR. & (var(1)%isat==0)) THEN ! no ponding ns(kk) = 1 ! start index for eqns ELSE ! ponding ns(kk) = 0 var(1)%phi = (one+e5)*var(1)%phie+(h0(kk)-hice(kk))*var(1)%Ksat TL(kk) = vmet(kk)%Ta ! initialise pond T vtop(kk)%isat = 1 vtop(kk)%h = h0(kk) vtop(kk)%phi = MAX((var(1)%phie -var(1)%he*var(1)%Ksat), & (one+e5)*var(1)%phie)+(h0(kk)-hice(kk))*var(1)%Ksat vtop(kk)%K = var(1)%Ksat vtop(kk)%rh = one vtop(kk)%lambdav = rlambda vtop(kk)%lambdaf = lambdaf ! var(1)%phi = var(1)%phie ! calculates phi1,eff (pond + top soil layer) !vh! does this get used??? CALL flux(par(1), vtop(kk), var(1), half*dx(1), & q(0), qya(0), qyb(0), qTa(0), qTb(0)) ENDIF IF (Sl(kk) >= one) THEN vlit(kk)%isat = 1 vlit(kk)%h = plit(kk)%he ENDIF ! get bottom boundary condn IF (botbc=="seepage") THEN vtmp = zerovars() vbot = SPREAD(vtmp,1,mp) IF (var(n)%h > -half*gf*dx(n)) THEN getqn(kk) = .TRUE. vbot(kk)%isat = 1 vbot(kk)%phi = (zero-var(n)%he)*var(n)%Ksat+var(n)%phie ! (special for frozen soil) vbot(kk)%K = var(n)%Ksat vbot(kk)%rh = one vbot(kk)%lambdav = rlambda vbot(kk)%lambdaf = lambdaf ELSE getqn(kk) = .FALSE. ENDIF END IF ! Fluxes and derivatives at the air/surface interface IF (vsnow(kk)%nsnow > 0) THEN surface_case(kk) = 2 ELSE surface_case(kk) = 1 ENDIF SELECT CASE (surface_case(kk)) CASE (1) ! no snow CALL SEB(n, par(:), vmet(kk), vsnow(kk), var(:), qprec(kk), qprec_snow(kk), dx(:), & h0(kk), Tsoil(:), & Tsurface(kk), G0(kk), lE0(kk),Epot(kk), & q(0), qevap(kk), qliq(0), qv(0), & qyb(0), qTb(0), qlyb(0), qvyb(0), qlTb(0), qvTb(0), qh(0), & qadv(0), qhyb(0), qhTb(0), qadvyb(0), qadvTb(0), irec) qya(0) = zero qTa(0) = zero qlya(0) = zero qvya(0) = zero qvTa(0) = zero qhya(0) = zero qhTa(0) = zero qadvya(0) = zero qadvTa(0) = zero CASE (2) ! snow CALL SEB(n, par(:), vmet(kk), vsnow(kk), var(:), qprec(kk), qprec_snow(kk), dx(:), & h0(kk), Tsoil(:), & Tsurface(kk), G0(kk), lE0(kk), Epot(kk), & q(-vsnow(kk)%nsnow), qevap(kk), qliq(-vsnow(kk)%nsnow), qv(-vsnow(kk)%nsnow), & qyb(-vsnow(kk)%nsnow), qTb(-vsnow(kk)%nsnow), qlyb(-vsnow(kk)%nsnow), & qvyb(-vsnow(kk)%nsnow), qlTb(-vsnow(kk)%nsnow), qvTb(-vsnow(kk)%nsnow), & qh(-vsnow(kk)%nsnow), qadv(-vsnow(kk)%nsnow), qhyb(-vsnow(kk)%nsnow), & qhTb(-vsnow(kk)%nsnow), qadvyb(-vsnow(kk)%nsnow), qadvTb(-vsnow(kk)%nsnow), irec) qya(-vsnow(kk)%nsnow) = zero qTa(-vsnow(kk)%nsnow) = zero qlya(-vsnow(kk)%nsnow) = zero qvya(-vsnow(kk)%nsnow) = zero qvTa(-vsnow(kk)%nsnow) = zero qhya(-vsnow(kk)%nsnow) = zero qhTa(-vsnow(kk)%nsnow) = zero qadvya(-vsnow(kk)%nsnow) = zero qadvTa(-vsnow(kk)%nsnow) = zero CASE default WRITE(*,*) "solve: illegal surface case." STOP 2 END SELECT ! surface_case qpme(kk) = q(0) ! water flux into top of soil column ! finished all the surfaces ! get moisture fluxes and derivatives (at time t=0, i.e. q0 etc.) CALL getfluxes_vp(n, dx(1:n), vtop(kk), vbot(kk), par(1:n), var(1:n), & ! moisture fluxes hint(1:n), phimin(1:n), q(0:n), qya(0:n), qyb(0:n), qTa(0:n), qTb(0:n), & qliq(0:n), qlya(0:n), qlyb(0:n), qv(0:n), qvT(0:n), qvh(0:n), qvya(0:n), & qvyb(0:n), & iflux(kk), init(kk), getq0(kk), getqn(kk), Tsoil(1:n), T0(kk), nsat(kk), nsatlast(kk)) qTa(n) = zero qTb(n) = zero qvTa(1:n) = qTa(1:n) qvTb(1:n) = qTb(1:n) qlTb(1:n) = zero ! get fluxes heat and derivatives (at time t=0, i.e. q0 etc.) CALL getheatfluxes(n, dx(1:n), dxL(kk), & qh(0:n), qhya(0:n), qhyb(0:n), qhTa(0:n), qhTb(0:n), & var(1:n), vlit(kk), Tsoil(1:n), TL(kk), litter, & q(0:n), qya(0:n), qyb(0:n), qTa(0:n), qTb(0:n), & qadv(0:n),qadvya(0:n), qadvyb(0:n), qadvTa(0:n), qadvTb(0:n), & advection) ! heat fluxes ! get heat and vapour fluxes and derivatives in snow-pack IF (vsnow(kk)%nsnow>1) THEN ! vapour flux at interfaces between snow layers DO j=1, vsnow(kk)%nsnow-1 keff = 2_r_2*((vsnow(kk)%kE(j)/(thousand*lambdaf))*(vsnow(kk)%kE(j+1)/(thousand*lambdaf))/ & ((vsnow(kk)%kE(j)/(thousand*lambdaf))*vsnow(kk)%depth(j+1)+(vsnow(kk)%kE(j+1)/ & (thousand*lambdaf))*vsnow(kk)%depth(j)) ) q(j-vsnow(kk)%nsnow) = keff*(vsnow(kk)%tsn(j)-vsnow(kk)%tsn(j+1)) qTa(j-vsnow(kk)%nsnow) = MERGE(keff,zero,vsnow(kk)%hliq(j)<=zero) qTb(j-vsnow(kk)%nsnow) = MERGE(-keff,zero,vsnow(kk)%hliq(j+1)<=zero) qya(j-vsnow(kk)%nsnow) = zero qyb(j-vsnow(kk)%nsnow) = zero ! conductive heat flux at interface between snow layers keff = 2_r_2*(vsnow(kk)%kth(j+1)*vsnow(kk)%kth(j))/ & (vsnow(kk)%kth(j+1)*vsnow(kk)%depth(j)+vsnow(kk)%kth(j)*vsnow(kk)%depth(j+1)) ! check this! qh(j-vsnow(kk)%nsnow) = keff*(vsnow(kk)%tsn(j)-vsnow(kk)%tsn(j+1)) qhTa(j-vsnow(kk)%nsnow) = MERGE(zero,keff,vsnow(kk)%hliq(j)>zero) qhTb(j-vsnow(kk)%nsnow) = MERGE(zero,-keff,vsnow(kk)%hliq(j+1)>zero) ! advective heat flux at interface between snow layers Tqw = MERGE(vsnow(kk)%tsn(j), vsnow(kk)%tsn(j+1), q(j-vsnow(kk)%nsnow)>zero) dTqwdTb = MERGE(zero,one, q(j-vsnow(kk)%nsnow)>zero) dTqwdTa = MERGE(one,zero, q(j-vsnow(kk)%nsnow)>zero) IF (vsnow(kk)%hliq(j)>zero) THEN qadv(j-vsnow(kk)%nsnow) = rhow*q(j-vsnow(kk)%nsnow)*cswat*Tqw qadvTa(j-vsnow(kk)%nsnow) = zero ELSE qadv(j-vsnow(kk)%nsnow) = rhow*q(j-vsnow(kk)%nsnow)*cswat*Tqw qadvTa(j-vsnow(kk)%nsnow) = rhow*cswat*q(j-vsnow(kk)%nsnow)*dTqwdTa + & rhow*cswat*Tqw*qTa(j-vsnow(kk)%nsnow) ENDIF qadvTb(0) = rhow*cswat*q(j-vsnow(kk)%nsnow)*dTqwdTb + rhow*cswat*Tqw*qTb(j-vsnow(kk)%nsnow) qh(j-vsnow(kk)%nsnow) = qh(j-vsnow(kk)%nsnow) + qadv(j-vsnow(kk)%nsnow) qhTa(j-vsnow(kk)%nsnow) = qhTa(j-vsnow(kk)%nsnow) + qadvTa(j-vsnow(kk)%nsnow) qhTb(j-vsnow(kk)%nsnow) = qhTb(j-vsnow(kk)%nsnow) + qadvTb(j-vsnow(kk)%nsnow) ENDDO ENDIF ! end fluxes at snow/snow interfaces IF (vsnow(kk)%nsnow.GE.1) THEN ! vapour flux at soil/snow interface IF (var(1)%isat==1) THEN keff = vsnow(kk)%kE(vsnow(kk)%nsnow)/(thousand*lambdaf)/vsnow(kk)%depth(vsnow(kk)%nsnow)/2_r_2 ENDIF IF (var(1)%isat==0) THEN keff = 2_r_2*((vsnow(kk)%kE(vsnow(kk)%nsnow)/(thousand*lambdaf))*(var(1)%kE/thousand/var(1)%lambdav))/ & ((vsnow(kk)%kE(vsnow(kk)%nsnow)/(thousand*lambdaf))*dx(1)+(var(1)%kE/thousand/var(1)%lambdav)* & vsnow(kk)%depth(vsnow(kk)%nsnow)) ENDIF q(0) = keff*(vsnow(kk)%tsn(vsnow(kk)%nsnow)-Tsoil(1)) qTa(0) = keff qTb(0) =-keff qya(0) = zero qyb(0) = zero qv(0) = q(0) qvyb(0) = qyb(0) qvTb(0) = qTb(0) qvya(0) = qya(0) qvTa(0) = qTa(0) qliq(0) = zero qlyb(0) = zero qlTb(0) = zero qlya(0) = zero qya(0) = zero; IF (vsnow(kk)%hliq(vsnow(kk)%nsnow)>zero) THEN qhTa(0) = zero qTa(0) = zero ENDIF ! conductive heat flux at snow/soil interface ! check this! keff = 2._r_2*(vsnow(kk)%kth(vsnow(kk)%nsnow)*var(1)%kth)/ & (vsnow(kk)%kth(vsnow(kk)%nsnow)*dx(1)+var(1)%kth*vsnow(kk)%depth(vsnow(kk)%nsnow)) qh(0) = keff*(vsnow(kk)%tsn(vsnow(kk)%nsnow)-Tsoil(1)) IF (vsnow(kk)%hliq(1)>zero) THEN qhTa(0) = zero ELSE qhTa(0) = keff ENDIF qhTb(0) = -keff ! advective heat flux at snow/soil interface Tqw = MERGE(vsnow(kk)%tsn(vsnow(kk)%nsnow), Tsoil(1), q(0)>zero) dTqwdTb = MERGE(zero,one, q(0)>zero) dTqwdTa = MERGE(one,zero, q(0)>zero) IF (vsnow(kk)%hliq(vsnow(kk)%nsnow)>zero) THEN qadv(0) = rhow*q(0)*cswat*Tqw qadvTa(0) = zero ELSE qadv(0) = rhow*q(0)*cswat*Tqw qadvTa(0) = rhow*cswat*q(0)*dTqwdTa + rhow*cswat*Tqw*qTa(0) ENDIF qadvTb(0) = rhow*cswat*q(0)*dTqwdTb + rhow*cswat*Tqw*qTb(0) ENDIF ! end of heat and vapour fluxes at soil/snow interface IF (ns(kk)==0) THEN ! pond included in top soil layer ! change qya(1) from dq/dphi (returned by getfluxes) to dq/dh qya(1) = var(1)%Ksat*qya(1) IF (advection==1) THEN qhya(1) = qhya(1) - qadvya(1) Tqw = MERGE(Tsoil(1), Tsoil(2), q(1)>zero) qadvya(1) = rhow*cswat*qya(1)*Tqw ! apply corrected qya(1) to qadvya(1) qhya(1) = qhya(1) + qadvya(1) ENDIF ENDIF ! adjust for bottom boundary condition IF (botbc=="zero flux") THEN qliq(n) = zero qv(n) = zero q(n) = zero qya(n) = zero qlya(n) = zero qvya(n) = zero ENDIF ! specify mositure flux at bottom of soil column (heat flux set to zero) IF (botbc /= "constant head") THEN SELECT CASE (botbc) CASE ("zero flux") q(n) = zero qya(n) = zero CASE ("free drainage") q(n) = gf*var(n)%K IF (var(n)%isat == 0) THEN qya(n) = gf*var(n)%KS ELSE qya(n) = zero END IF CASE ("seepage") IF (var(n)%h <= -half*gf*dx(n)) THEN q(n) = zero qya(n) = zero END IF CASE default WRITE(*,*) "solve: illegal bottom boundary condition." STOP 2 END SELECT END IF IF (PRESENT(qali)) THEN IF (qali(kk)>zero) THEN q(n) = -qali(kk) qya(n) = zero END IF ENDIF IF (experiment==7 .OR. experiment==8) THEN q(n) = q(0) qya(n) = zero ENDIF IF (experiment==8) qh(n) = G0(kk) ! adjust lower heat flux for advection IF (advection==1) THEN qadv(n) = rhow*cswat*(Tsoil(n))*q(n) qadvya(n) = rhow*cswat*(Tsoil(n))*qya(n) qadvTa(n)= rhow*cswat*q(n) qh(n) = qh(n) + qadv(n) qhya(n) = qhya(n) + qadvya(n) qhTa(n) = qhTa(n) + qadvTa(n) ELSE qadv(:) = zero qadvya(:) = zero qadvyb(:) = zero qadvTa(:) = zero qadvTb(:) = zero ENDIF qexd(1:n) = zero ! time derivative for root extraction (assumed fixed at input value) again(kk) = .FALSE. ! flag for recalcn of fluxes (default=false) !----- end get fluxes and derivs END SUBROUTINE get_fluxes_and_derivs SUBROUTINE estimate_timestep( & irec, tfin, mp, n, dx, h0, & qh, qhTa, qadvTa, qhTb, qadvTb, & nsteps, vlit, vsnow, var, & dxL, & plit, par, & deltaTa, & qL, qhL, & again, ns, nsat, & nsatlast, nsteps0, dmax, dt, & qpme, t, & q, qadv, & tmp2d1, tmp2d2, & tmp1d1, tmp1d3, & iflux, litter, kk, & advection, & iqex & ) IMPLICIT NONE REAL(r_2) :: tfin INTEGER(i_d) :: mp, irec INTEGER(i_d) :: n REAL(r_2), DIMENSION(1:n) :: dx REAL(r_2), DIMENSION(1:mp) :: h0 REAL(r_2), DIMENSION(-nsnow_max:n) :: qh INTEGER(i_d), DIMENSION(1:mp) :: nsteps TYPE(vars), DIMENSION(1:mp) :: vlit TYPE(vars_snow), DIMENSION(1:mp) :: vsnow TYPE(vars), DIMENSION(1:n) :: var REAL(r_2), DIMENSION(1:mp) :: dxL TYPE(params), DIMENSION(1:mp) :: plit TYPE(params), DIMENSION(1:n) :: par REAL(r_2), DIMENSION(1:mp) :: deltaTa REAL(r_2), DIMENSION(1:mp) :: qL, qhL LOGICAL, DIMENSION(1:mp) :: again INTEGER(i_d), DIMENSION(1:mp) :: ns, nsat, nsatlast, nsteps0 REAL(r_2), DIMENSION(1:mp) :: dmax, dt REAL(r_2), DIMENSION(1:mp) :: qpme, t REAL(r_2), DIMENSION(-nsnow_max:n) :: q, qhTa, qhTb REAL(r_2), DIMENSION(-nsnow_max:n) :: qadv, qadvTa, qadvTb REAL(r_2), DIMENSION(0:n) :: tmp2d1, tmp2d2, deltaTmax REAL(r_2), DIMENSION(1:mp) :: tmp1d1, tmp1d3 INTEGER(i_d), DIMENSION(1:mp) :: iflux LOGICAL :: litter INTEGER(i_d) :: kk INTEGER(i_d) :: advection ! switches REAL(r_2), DIMENSION(1:n) :: iqex !----- first estimate of time step dt before the calculation ! gets revised after the calculation dmax(kk) = zero tmp2d1(:) = zero ! 1st order estim of rate of change of moisture storage tmp2d2(:) = zero ! 1st order estim of rate of change of temperature ! estimate rate of change of moisture storage [m/s] WHERE (var(1:n)%isat==0.AND.var(1:n)%iice==0.) tmp2d1(1:n) = & ABS(q(1:n)-q(0:n-1)-iqex(1:n))/(par(1:n)%thre*dx(1:n)) WHERE (var(1:n)%iice==1) tmp2d1(1:n) = tmp2d1(1:n)/2._r_2 ! estimate rate of change of temperature [K/s] tmp2d2(1:n) = ABS(qh(1:n)-qh(0:n-1))/(var(1:n)%csoileff*dx(1:n)) IF (advection==1) THEN ! first order estimate of rate of temperature change tmp2d2(1:n) = ABS((qh(1:n)-qadv(1:n))-(qh(0:n-1)-qadv(0:n-1)))/(var(1:n)%csoileff*dx(1:n)) deltaTmax(1:n) = tmp2d2(1:n)*(tfin-t(kk)) ! maximum T change ENDIF IF (litter .AND. ns(kk)==1 ) THEN ! litter , no pond IF (vlit(kk)%isat==0) THEN ! estimate rate of change of moisture storage [m/s] WRITE(*,*) 'Should not be here - QL 01 ', qL(kk) tmp2d1(0) = ABS(q(0) - qL(kk))/(plit(kk)%thre*dxL(kk)) ENDIF WRITE(*,*) 'Should not be here - QHL 01 ', qhL(kk) tmp2d2(0) = ABS(qh(0) - qhL(kk))/(vlit(kk)%csoileff*dxL(kk)) ! estimate rate of change of heat storage [K/s] tmp1d3(kk) = (dTLmax-ABS(deltaTa(kk))) / tmp2d2(0) ELSE tmp2d1(0) = zero tmp2d2(0) = zero tmp1d3(kk) = dtmax ENDIF dmax(kk) = MAXVAL(tmp2d1(1:n),1) ! max derivative |dS/dt| tmp1d1(kk) = MAXVAL(tmp2d2(1:n),1) ! max derivative |dTsoil/dt| IF (dmax(kk) > zero) THEN dt(kk) = MIN(dSmax/dmax(kk), dTLmax/tmp1d1(kk), tmp1d3(kk)) ! constrained either by moisture or temp ELSE ! steady state flow IF (qpme(kk)>=q(n)) THEN ! if saturated soil columnn and more precip then drainige -> finish dt(kk) = tfin-t(kk) ! step to finish ELSE ! otherwise adjust dt because change of pond height dt(kk) = -(h0(kk)-half*h0min)/(qpme(kk)-q(n)) END IF dt(kk) = MIN(dt(kk), dTLmax/tmp1d1(kk), tmp1d3(kk)) ! constrained by temp END IF ! check that time step is short enough to prevent melting a top snow layer IF (vsnow(kk)%nsnow.GT.0 ) THEN ! energy required to melt ice in snow-pack tmp1d1(kk) = -rhow*(vsnow(kk)%hsnow(1)-vsnow(kk)%hliq(1))*(csice*vsnow(kk)%tsn(1) - lambdaf) IF (qh(-vsnow(kk)%nsnow)*dt(kk).GT.tmp1d1(kk)) THEN dt(kk) = 0.9_r_2*tmp1d1(kk)/qh(-vsnow(kk)%nsnow) ENDIF ENDIF IF (dt(kk)>dtmax) dt(kk) = dtmax ! user's limit ! if initial step, improve phi where S>=1 ! might be that you get better derivatives, especially at sat/non-sat interfaces IF (nsteps(kk)==nsteps0(kk) .AND. nsat(kk)>0 .AND. iflux(kk)==1) THEN again(kk) = .TRUE. dt(kk) = dtmin END IF ! if fully saturated but was not fully saturated before, adjust even within each time step iteration IF (nsat(kk)==n .AND. nsatlast(kk)<n .AND. iflux(kk)==1) THEN again(kk) = .TRUE. ! profile has just become saturated so adjust phi values dt(kk) = dtmin END IF ! sprint to the end IF (t(kk)+1.1_r_2*dt(kk)>tfin) THEN ! step to finish dt(kk) = tfin-t(kk) t(kk) = tfin ELSE t(kk) = t(kk)+dt(kk) ! tentative update IF (again(kk)) t(kk) = t(kk)-dt(kk) END IF !----- end estimate time step dt END SUBROUTINE estimate_timestep SUBROUTINE iflux_loop( & tfin, irec, mp, qprec, qprec_snow, n, dx, h0, S, thetai, & Jsensible, Tsoil, evap, infil, drainage, discharge, & qh, nsteps, vmet, vlit, vsnow, var, T0, Tsurface, Hcum, lEcum, & Gcum, Qadvcum, Jcol_sensible, & Jcol_latent_S, Jcol_latent_T, csoil, kth, phi, dxL, zdelta, SL, Tl, & plit, par, wex, ciso_snow, cisoice_snow, & qali, & qvsig, qlsig, qvTsig, qvh, deltaTa, & precip, qevap, qL, qhL, qybL, qTbL, qhTbL, qhybL, rexcol, wcol, & again, getq0,getqn,init, again_ice, ih0, iok, itmp, ns, nsat, & nsatlast, nsteps0, accel, dmax, dt, dwinfil, dwoff, fac, & phip, qpme, rsig, rsigdt, sig, t, hint, phimin, & qexd, aa, bb, cc, dd, ee, ff, gg, dy, aah, bbh, cch, ddh, eeh, ffh, ggh, & de, q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, qadv, qadvya, qadvyb, & qadvTa, qadvTb, vtmp, qsig, qhsig, qadvsig, qliq, qv, qvT, qlya, qlyb, & qvya, qvyb, qlTb, qvTa, qvTb, deltaS, dTsoil, tmp2d1, tmp2d2, & cv0, & cisoliqice_snow, & dthetaldT, thetal, isave, nsteps_ice, imelt, vtop, vbot, v_aquifer, & dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge, & dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S, & deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T, qevapsig, & qrunoff, tmp1d1, tmp1d2, tmp1d3, tmp1d4, deltah0, deltaSL, & lE0, G0, & Epot, Tfreezing, dtdT, LHS, RHS, LHS_h, RHS_h, surface_case, nns, & iflux, litter, i, j, k, kk, condition, littercase, isotopologue, & advection, c2, theta, dTqwdTa, dTqwdTb, Tqw, keff, cp, & cpeff, hice, h0_0, hice_0, h0_tmp, hice_tmp, qmelt, & qtransfer, delta_snowcol, delta_snowT, & delta_snowliq, thetai_0, J0, tmp1, tmp2, iqex, & nfac1, nfac2, nfac3, nfac4, nfac5, nfac6, nfac7, nfac8, nfac9, & nfac10, nfac11, nfac12, J0snow, wcol0snow, h_ex, wpi, err & ) IMPLICIT NONE REAL(r_2) :: tfin INTEGER(i_d) :: irec, mp REAL(r_2), DIMENSION(1:mp) :: qprec REAL(r_2), DIMENSION(1:mp) :: qprec_snow INTEGER(i_d) :: n REAL(r_2), DIMENSION(1:n) :: dx REAL(r_2), DIMENSION(1:mp) :: h0 REAL(r_2), DIMENSION(1:n) :: S REAL(r_2), DIMENSION(1:n) :: thetai REAL(r_2), DIMENSION(1:n) :: Jsensible REAL(r_2), DIMENSION(1:n) :: Tsoil REAL(r_2), DIMENSION(1:mp) :: evap REAL(r_2), DIMENSION(1:mp) :: infil REAL(r_2), DIMENSION(1:mp) :: drainage, discharge REAL(r_2), DIMENSION(1:mp,-nsnow_max:n) :: qh INTEGER(i_d), DIMENSION(1:mp) :: nsteps TYPE(vars_met), DIMENSION(1:mp) :: vmet TYPE(vars), DIMENSION(1:mp) :: vlit TYPE(vars_snow), DIMENSION(1:mp) :: vsnow TYPE(vars), DIMENSION(1:n) :: var REAL(r_2), DIMENSION(1:mp) :: T0, Tsurface REAL(r_2), DIMENSION(1:mp) :: Hcum, lEcum REAL(r_2), DIMENSION(1:mp) :: Gcum, Qadvcum REAL(r_2), DIMENSION(1:mp) :: Jcol_sensible, Jcol_latent_S, Jcol_latent_T REAL(r_2), DIMENSION(1:n) :: csoil, kth REAL(r_2), DIMENSION(1:n) :: phi REAL(r_2), DIMENSION(1:mp) :: dxL REAL(r_2), DIMENSION(1:mp) :: zdelta, SL, Tl TYPE(params), DIMENSION(1:mp) :: plit TYPE(params), DIMENSION(1:n) :: par REAL(r_2), DIMENSION(1:mp,1:n), OPTIONAL :: wex REAL(r_2), DIMENSION(1:mp,1:nsnow_max) :: ciso_snow REAL(r_2), DIMENSION(1:mp,1:nsnow_max) :: cisoice_snow REAL(r_2), DIMENSION(1:mp), OPTIONAL :: qali REAL(r_2), DIMENSION(1:mp,-nsnow_max:n) :: qvsig, qlsig, qvTsig, qvh REAL(r_2), DIMENSION(1:mp) :: deltaTa REAL(r_2), DIMENSION(1:mp) :: precip, qevap REAL(r_2), DIMENSION(1:mp) :: qL, qhL, qybL, qTbL, qhTbL, qhybL, rexcol, wcol LOGICAL, DIMENSION(1:mp) :: again, getq0,getqn,init LOGICAL, DIMENSION(1:mp,1:n) :: again_ice INTEGER(i_d), DIMENSION(1:mp) :: ih0, iok, itmp, ns, nsat, nsatlast, nsteps0 REAL(r_2), DIMENSION(1:mp) :: accel, dmax, dt, dwinfil, dwoff, fac, phip REAL(r_2), DIMENSION(1:mp) :: qpme, rsig, rsigdt, sig, t REAL(r_2), DIMENSION(1:n) :: hint, phimin, qexd REAL(r_2), DIMENSION(-nsnow_max+1:n) :: aa, bb, cc, dd, ee, ff, gg, dy REAL(r_2), DIMENSION(-nsnow_max+1:n) :: aah, bbh, cch, ddh, eeh, ffh, ggh, de REAL(r_2), DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb REAL(r_2), DIMENSION(-nsnow_max:n) :: qadv, qadvya, qadvyb, qadvTa, qadvTb TYPE(vars) :: vtmp REAL(r_2), DIMENSION(-nsnow_max:n) :: qsig, qhsig, qadvsig REAL(r_2), DIMENSION(-nsnow_max:n) :: qliq, qv, qvT, qlya, qlyb, qvya, qvyb, qlTb, qvTa, qvTb REAL(r_2), DIMENSION(1:n) :: deltaS, dTsoil REAL(r_2), DIMENSION(0:n) :: tmp2d1, tmp2d2 REAL(r_2), DIMENSION(1:n) :: cv0 REAL(r_2), DIMENSION(1:nsnow_max) :: cisoliqice_snow REAL(r_2), DIMENSION(1:n) :: dthetaldT, thetal INTEGER(i_d), DIMENSION(1:n) :: isave, nsteps_ice, imelt TYPE(vars), DIMENSION(1:mp) :: vtop, vbot TYPE(vars_aquifer), DIMENSION(1:mp) :: v_aquifer REAL(r_2), DIMENSION(1:mp) :: dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge REAL(r_2), DIMENSION(1:mp) :: dJcol_latent_S, dJcol_latent_T, dJcol_sensible REAL(r_2), DIMENSION(1:n) :: deltaJ_latent_S, deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T REAL(r_2), DIMENSION(1:mp) :: qevapsig REAL(r_2), DIMENSION(1:mp) :: qrunoff REAL(r_2), DIMENSION(1:mp) :: tmp1d1, tmp1d2, tmp1d3, tmp1d4 REAL(r_2), DIMENSION(1:mp) :: deltah0 REAL(r_2), DIMENSION(1:mp) :: deltaSL REAL(r_2), DIMENSION(1:mp) :: lE0, G0, Epot REAL(r_2), DIMENSION(1:mp) :: Tfreezing REAL(r_2), DIMENSION(1:mp) :: dtdT REAL(r_2), DIMENSION(-nsnow_max+1:n) :: LHS, RHS, LHS_h, RHS_h INTEGER(i_d), DIMENSION(1:mp) :: surface_case INTEGER(i_d), DIMENSION(1:mp) :: nns, iflux LOGICAL :: litter INTEGER(i_d) :: i, j, k, kk, condition INTEGER(i_d) :: littercase, isotopologue, advection ! switches REAL(r_2) :: c2, theta REAL(r_2) :: dTqwdTa, dTqwdTb, Tqw, keff REAL(r_2), DIMENSION(1:mp) :: cp, cpeff, hice, h0_0, hice_0, h0_tmp, hice_tmp REAL(r_2), DIMENSION(nsnow_max) :: qmelt, hsnow REAL(r_2), DIMENSION(1:mp) :: qtransfer REAL(r_2), DIMENSION(nsnow_max) :: delta_snowcol, delta_snowT, delta_snowliq REAL(r_2), DIMENSION(1:n) :: thetai_0, J0 REAL(r_2) :: tmp1, tmp2 REAL(r_2), DIMENSION(1:n) :: iqex INTEGER(i_d), DIMENSION(1:mp) :: nfac1, nfac2, nfac3, nfac4, nfac5, & nfac6, nfac7, nfac8, nfac9, nfac10, nfac11, nfac12 REAL(r_2), DIMENSION(1:mp) :: J0snow, wcol0snow REAL(r_2), DIMENSION(1:n) :: h_ex REAL(r_2) :: wpi ! Error flag if nstep of SLI > nsteps_max: err=0 -> no error; err/=0 -> error INTEGER(i_d), INTENT(INOUT), OPTIONAL :: err hsnow(:) = zero DO WHILE (again(kk)) ! sometimes need twice to adjust phi at satn nsatlast(kk) = nsat(kk) ! for detecting onset of profile saturation nsat(kk) = SUM(var(:)%isat,1) ! no. of sat layers sig(kk) = half IF (nsat(kk) /= 0) sig(kk) = one ! time weighting sigma rsig(kk) = one/sig(kk) ! update variables IF (iflux(kk)==1) THEN ! Calc flux matric potentials (and derivatives) from S ! this set var-structure isave(:) = var(:)%isat ! Debug for mp=1: remove elemental from hyofS and do loop instead of next line CALL hyofS(S(1:n), Tsoil(1:n), par(1:n), var(1:n)) !do i=1, n ! call hyofS(S(i), Tsoil(i), par(i), var(i)) !end do ! End debug hyofS cp(kk) = REAL(1-var(1)%iice,r_2)*cswat*rhow & ! heat capacity of pond + REAL(var(1)%iice,r_2)*rhow* & ((one-var(1)%thetai/par(1)%thre)*cswat + (var(1)%thetai/par(1)%thre)*csice) cpeff(kk) = cp(kk) + rhow*lambdaf*var(1)%dthetaldT/par(1)%thre ! phip(kk) = max(var(1)%phie-var(1)%he*var(1)%Ksat, (one+e5)*var(1)%phie) !at onset of ponding phip(kk) = (one+e5)*var(1)%phie !at onset of ponding var(:)%isat = isave(:) thetai(:) = var(:)%thetai ! ice profile thetal(:) = var(:)%thetal ! liq water profile thetai_0(:) = thetai(:) ! initial ice profile dthetaldT(:) = var(:)%dthetaldT hice(kk) = h0(kk)*var(1)%thetai/par(1)%thre hice_0(kk) = hice(kk) h0_0(kk) = h0(kk) ! sensible heat stored in top layer + pond Jsensible(1) = (var(1)%csoil* dx(1)+h0(kk)*cp(kk))*(Tsoil(1)) ! sensible heat stored in soil column (2:n) Jsensible(2:n) = var(2:n)%csoil*(Tsoil(2:n))* dx(2:n) ! calculate this here in case snow pack disappears: ciso of transferred water needs to be retained IF ((isotopologue /= 0).AND.(vsnow(kk)%hsnow(1).GT.zero)) THEN cisoliqice_snow(1:nsnow_max) = (ciso_snow(kk,1:nsnow_max)*vsnow(kk)%hliq(1:nsnow_max) + & cisoice_snow(kk,1:nsnow_max)*(vsnow(kk)%hsnow(1:nsnow_max)-vsnow(kk)%hliq(1:nsnow_max))) / & vsnow(kk)%hsnow(1:nsnow_max) ENDIF CALL snow_adjust(irec, mp, n, kk, ns, h0, hice, thetai(:), dx(:), vsnow, var, par, S(:), Tsoil(:), & Jcol_latent_S, Jcol_latent_T, Jcol_sensible, deltaJ_sensible_S(:), qmelt(:), qtransfer, j0snow) thetai(1) = var(1)%thetai ! this is the value of thetaice prior to matrix call CALL hyofS(S(1), Tsoil(1), par(1), var(1)) ENDIF ! iflux==1 ! initialise litter vars IF (iflux(kk)==1 .AND. (littercase==1 .OR. littercase == 2)) THEN ! call litter_props(Sl(kk), Tl(kk), vlit(kk), plit(kk), h0(kk)) CALL litter_props(Sl(kk), Tsoil(1), vlit(kk), plit(kk), h0(kk)) ENDIF ! ! phi is solution var at satn, so h calc from phi where S>=1 - done in hyofS above for S<1 ! where (S(:) >= one) & ! (special for frozen soil) ! var(:)%h = var(:)%he + (var(:)%phi-var(:)%phie)/var(:)%Ksat CALL get_fluxes_and_derivs( & irec, mp, qprec, qprec_snow, n, dx(:), h0, & Tsoil(:), S(:), & qh(kk,:), vmet, vlit, vsnow, var, T0, Tsurface, & dxL, SL, Tl, & plit, par, & qali, & qvh(kk,:), & qevap, & again, getq0,getqn,init, ns, nsat, & nsatlast, & phip, qpme, hint(:), phimin(:), & qexd(:), & q(:), qya(:), qyb(:), qTa(:), qTb(:),qhya(:), qhyb(:), qhTa(:), qhTb(:), qadv(:), qadvya(:), qadvyb(:), & qadvTa(:), qadvTb(:), vtmp, qliq(:), qv(:), qvT(:), qlya(:), qlyb(:), & qvya(:), qvyb(:), qlTb(:), qvTa(:), qvTb(:), & vtop, vbot, & lE0, G0, & Epot, surface_case, & iflux, litter, j, kk, & advection, dTqwdTa, dTqwdTb, Tqw, keff, & hice & ) CALL estimate_timestep( & irec, tfin, mp, n, dx(:), h0, & qh(kk,:), qhTa(:), qadvTa(:), qhTb(:), qadvTb(:), & nsteps, vlit, vsnow, var, & dxL, & plit, par, & deltaTa, & qL, qhL, & again, ns, nsat, & nsatlast, nsteps0, dmax, dt, & qpme, t, & q(:), qadv(:), & tmp2d1(:), tmp2d2(:), & tmp1d1, tmp1d3, & iflux, litter, kk, & advection, & iqex(:) & ) !----- get and solve eqns rsigdt(kk) = one/(sig(kk)*dt(kk)) IF (.NOT. again(kk)) THEN dwoff(kk) = MAX(h0(kk)*(one-var(1)%thetai/par(1)%thre)-h0max,zero) qrunoff(kk) = dwoff(kk)/dt(kk) ELSE qrunoff(kk) = zero ENDIF ! aa, bb, cc and dd hold coeffs and rhs of linear equation eqn set CALL get_and_solve_eqn( & irec, mp, qprec, n, dx(:), h0, S(:), thetai(:), & Tsoil(:), & qh(kk,:), nsteps, vlit, vsnow, var, & dxL, & plit, par, & deltaTa, & qevap, qL, qhL, qybL, qTbL, qhTbL, qhybL, & again, again_ice, iok, itmp, ns, & accel, dt, fac, & phip, rsig, rsigdt, sig, t, & qexd(:), aa(:), bb(:), cc(:), dd(:), ee(:), ff(:), gg(:), dy(:), & aah(:), bbh(:), cch(:), ddh(:), eeh(:), ffh(:), ggh(:), & de(:), q(:), qya(:), qyb(:), qTa(:), qTb(:),qhya(:), qhyb(:), qhTa(:), qhTb(:), qadv(:), qadvya(:), qadvyb(:), & qadvTa(:), qadvTb(:), qsig(:), qhsig(:), qadvsig(:), & dTsoil(:), tmp2d1(:), & dthetaldT(:), nsteps_ice, imelt, & deltaJ_latent_T(:), deltaJ_sensible_S(:), & qrunoff, tmp1d1, tmp1d2, tmp1d3, tmp1d4, & G0, & Tfreezing, LHS(:), RHS(:), LHS_h(:), RHS_h(:), nns, & iflux, litter, i, kk, condition, & advection, c2, theta, cp, & cpeff, hice, h0_tmp, hsnow(:), & delta_snowT(:), & delta_snowliq(:), J0(:), iqex(:), & nfac1, nfac2, nfac3, nfac4, nfac5, nfac6, nfac7, nfac8, nfac9, & nfac10, nfac11, nfac12, err & ) IF (err /= 0) RETURN CALL update_unknowns( & irec, mp, qprec, qprec_snow, n, dx(:), h0, S(:), & Tsoil(:), evap, infil, drainage, discharge, & qh(kk,:), nsteps, vmet, vlit, vsnow, var, Hcum, lEcum, & Gcum, Qadvcum, & csoil(:), kth(:), phi(:), zdelta, SL, Tl, & par, wex, & qvsig(kk,:), qlsig(kk,:), qvTsig(kk,:), & precip, qevap, rexcol, wcol, & again, ih0, ns, & dt, dwinfil, dwoff, & sig, & qexd(:), dy(:), & de(:), q(:), qya(:), qyb(:), qTa(:), qTb(:),qhya(:), qhyb(:), qhTa(:), qhTb(:), & qsig(:), qhsig(:), qadvsig(:), qliq(:), qv(:), qvT(:), qlyb(:), & qvya(:), qvyb(:), qlTb(:), qvTa(:), qvTb(:), dTsoil(:), & v_aquifer, & dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge, & dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S(:), & deltaJ_latent_T(:), deltaJ_sensible_S(:), deltaJ_sensible_T(:), qevapsig, & qrunoff, tmp1d1, tmp1d2, tmp1d3, deltah0, deltaSL, & LHS_h(:), surface_case, & litter, i, j, kk, & advection, theta, Tqw, cp, & hsnow(:), & delta_snowcol(:), delta_snowT(:), & delta_snowliq(:), J0(:), iqex(:), & J0snow, wcol0snow & ) CALL update_s_t( & mp, n, dx(:), h0, S(:), thetai(:), & Tsoil(:), infil, & nsteps, var, & par, & qlsig(kk,:), & again, ih0, ns, & dt, & sig, & dy(:), & qhsig(:), & deltaS(:), dTsoil(:), & dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S(:), & deltaJ_latent_T(:), deltaJ_sensible_S(:), deltaJ_sensible_T(:), & qrunoff, tmp1d1, tmp1d2, tmp1d3, tmp1d4, deltah0, & Tfreezing, dtdT, LHS_h(:), & i, j, k, kk, & theta, cp, & hice, hice_tmp, & J0(:), tmp1, tmp2, & h_ex, wpi & ) IF (.NOT. again(kk)) THEN cv0(1:n) = var(1:n)%cv ! save cv before updating isave(1:n) = var(1:n)%iice ! save isat before updating ! update variables (particularly thetaice) ! Debug for mp=1: remove elemental from hyofS and do loop instead of next line CALL hyofS(S(1:n), Tsoil(1:n), par(1:n), var(1:n)) !do i=1, n ! call hyofS(S(i), Tsoil(i), par(i), var(i)) !end do ! End debug hyofS var(1:n)%iice = isave(1:n) ENDIF ! if .not.again iflux(kk) = iflux(kk) + 1 END DO ! do while (again(kk)) ! iflux loop END SUBROUTINE iflux_loop SUBROUTINE update_s_t( & mp, n, dx, h0, S, thetai, & Tsoil, infil, & nsteps, var, & par, & qlsig, & again, ih0, ns, & dt, & sig, & dy, & qhsig, & deltaS, dTsoil, & dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S, & deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T, & qrunoff, tmp1d1, tmp1d2, tmp1d3, tmp1d4, deltah0, & Tfreezing, dtdT, LHS_h, & i, j, k, kk, & theta, cp, & hice, hice_tmp, & J0, tmp1, tmp2, & h_ex, wpi & ) IMPLICIT NONE INTEGER(i_d) :: mp INTEGER(i_d) :: n REAL(r_2), DIMENSION(1:n) :: dx REAL(r_2), DIMENSION(1:mp) :: h0 REAL(r_2), DIMENSION(1:n) :: S REAL(r_2), DIMENSION(1:n) :: thetai REAL(r_2), DIMENSION(1:n) :: Tsoil REAL(r_2), DIMENSION(1:mp) :: infil INTEGER(i_d), DIMENSION(1:mp) :: nsteps TYPE(vars), DIMENSION(1:n) :: var TYPE(params), DIMENSION(1:n) :: par REAL(r_2), DIMENSION(-nsnow_max:n) :: qlsig LOGICAL, DIMENSION(1:mp) :: again INTEGER(i_d), DIMENSION(1:mp) :: ih0, ns REAL(r_2), DIMENSION(1:mp) :: dt REAL(r_2), DIMENSION(1:mp) :: sig REAL(r_2), DIMENSION(-nsnow_max+1:n) :: dy REAL(r_2), DIMENSION(-nsnow_max:n) :: qhsig REAL(r_2), DIMENSION(1:n) :: deltaS, dTsoil REAL(r_2), DIMENSION(1:mp) :: dJcol_latent_S, dJcol_latent_T, dJcol_sensible REAL(r_2), DIMENSION(1:n) :: deltaJ_latent_S, deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T REAL(r_2), DIMENSION(1:mp) :: qrunoff REAL(r_2), DIMENSION(1:mp) :: tmp1d1, tmp1d2, tmp1d3, tmp1d4 REAL(r_2), DIMENSION(1:mp) :: deltah0 REAL(r_2), DIMENSION(1:mp) :: Tfreezing REAL(r_2), DIMENSION(1:mp) :: dtdT REAL(r_2), DIMENSION(-nsnow_max+1:n) :: LHS_h INTEGER(i_d) :: i, j, k, kk REAL(r_2) :: theta REAL(r_2), DIMENSION(1:mp) :: cp, hice, hice_tmp REAL(r_2), DIMENSION(1:n) :: J0 REAL(r_2) :: tmp1, tmp2 REAL(r_2), DIMENSION(1:n) :: h_ex REAL(r_2) :: wpi ! update variables (S,T) to end of time step ! update variables to sig for use in isotope routine DO i=1, n IF (.NOT.again(kk)) Tsoil(i) = Tsoil(i) + dTsoil(i) IF (var(i)%isat==0) THEN IF (.NOT.again(kk)) THEN deltaS(i) = dy(i) !required for isotope subroutine S(i) = S(i)+dy(i) IF (S(i)>one .AND. dy(i)>zero) THEN ! onset of saturation of layer var(i)%isat = 1 var(i)%K = var(i)%Ksat var(i)%phi = var(i)%phie ENDIF ENDIF ELSE deltaS(i) = zero !required for isotope subroutine IF (i==1) THEN var(i)%phi = var(i)%phi + dy(i)*REAL(ns(kk),r_2) ! pond included in top soil layer ELSE var(i)%phi = var(i)%phi + dy(i) ENDIF ! var(i)%phi = zero from Ross ! VH thinks it is o.k. because hyofS is called later again ! MC think that we should probably get rid of -he*Ksat in phip etc. ! because we merged the pond with the first layer. IF (i==1 .AND. ih0(kk)/=0 .AND. var(i)%phi>=var(i)%phie) var(i)%phi = zero ! pond gone IF (var(i)%phi<var(i)%phie .AND. var(i)%iice==0) THEN ! desaturation of layer var(i)%isat = 0 var(i)%K = var(i)%Ksat var(i)%phi = var(i)%phie var(i)%KS = par(i)%KSe var(i)%phiS = par(i)%phiSe ELSEIF (var(i)%phi<var(i)%phie .AND. var(i)%iice==1) THEN var(i)%isat = 0 var(i)%K = var(i)%Ksat var(i)%phi = var(i)%phie var(i)%KS = var(i)%KS var(i)%phiS = var(i)%phiS ENDIF END IF IF (.NOT. again(kk)) THEN !$ if (h0(kk)<zero .and. var(1)%isat==0 .and. (i==1)) then ! start negative pond correction !$ infil(kk) = infil(kk)+h0(kk) !$ ! negative pond converted to loss of soil moisture in top two layers !$ S(1) = S(1) + h0(kk)/(dx(1)*par(1)%thre)* dx(1)/(dx(1)+dx(2)) !$ deltaS(1) = h0(kk)/(dx(1)*par(1)%thre)* dx(1)/(dx(1)+dx(2)) !$ ! adjust flux between layer 1 and layer 2 !$ qlsig(1) = qlsig(1) + h0(kk)* dx(2)/(dx(1)+dx(2))/dt(kk) !$ deltaS(2) = h0(kk)/(dx(2)*par(2)%thre)* dx(2)/(dx(1)+dx(2)) !$ if (var(2)%isat.eq.0) then !$ dy(2) = dy(2) + deltaS(2) !$ else !$ dy(2) = deltaS(2) !$ var(2)%isat = 0 !$ endif !$ !$ deltah0(kk) = -(h0(kk)-deltah0(kk)) ! whole pond lost (new change = - original pond height) !$ h0(kk) = zero ! zero pond remaining !$ endif ! water in profile initially wpi = SUM((par(:)%thr + (par(:)%the-par(:)%thr)*S(:))*dx(:)) IF (h0(kk)<zero .AND. var(1)%isat==0 .AND. (i==1)) THEN ! start negative pond correction wpi = SUM((par(:)%thr + (par(:)%the-par(:)%thr)*S(:))*dx(:)) infil(kk) = infil(kk)+h0(kk) DO j=1,n h_ex(j) = -h0(kk)* (par(j)%the-par(j)%thr)*S(j)* dx(j)/wpi ENDDO DO j=1,n-1 qlsig(j) = qlsig(j) +SUM(h_ex(j:n)) ENDDO DO j = 1,n IF (j<n .AND. ( S(j)*(dx(j)*par(j)%thre) + h_ex(j)).LT.zero) THEN h_ex(j+1) = h_ex(j+1) +(h_ex(j)- S(j)*(dx(j)*par(j)%thre)) h_ex(j) = S(j)*(dx(j)*par(j)%thre) ENDIF deltaS(j)= -MIN(h_ex(j)/(dx(j)*par(j)%thre),S(j)) IF (var(j)%isat.EQ.0) THEN dy(j) = dy(j) + deltaS(j) ELSE dy(j) = deltaS(j) var(j)%isat = 0 ENDIF IF (j==1) S(j)=S(j)+deltaS(j) ENDDO deltah0(kk) = -(h0(kk)-deltah0(kk)) ! whole pond lost (new change = - original pond height) h0(kk) = zero ! zero pond remaining ENDIF IF (S(i)<zero .AND. i<n) THEN ! start correction for negative soil moisture in any layer wpi = SUM((par(i+1:n)%thr + (par(i+1:n)%the-par(i+1:n)%thr)*S(i+1:n))*dx(i+1:n)) DO j=i+1,n h_ex(j) = -deltaS(i)*dx(i)*par(i)%thre * & (par(j)%thr + (par(j)%the-par(j)%thr)*S(j))*dx(j) /wpi ENDDO DO j=i+1,n qlsig(j) = qlsig(j) +SUM(h_ex(j:n)) ENDDO DO j = i+1,n IF (j<n .AND. ( S(j)*(dx(j)*par(j)%thre) - h_ex(j)).LT.zero) THEN h_ex(j+1) = h_ex(j+1) +(h_ex(j)- S(j)*(dx(j)*par(j)%thre)) h_ex(j) = S(j)*(dx(j)*par(j)%thre) ENDIF deltaS(j)= -MIN(h_ex(j)/(dx(j)*par(j)%thre),S(j)) IF (var(j)%isat.EQ.0) THEN dy(j) = dy(j) + deltaS(j) ELSE dy(j) = deltaS(j) var(j)%isat = 0 ENDIF IF (j==i) S(j)=S(j)+deltaS(j) ENDDO S(i)=S(i)-deltaS(i) deltaS(i) = zero ENDIF ! corrections for onset of freezing and melting ! calculate freezing point temperature at new S Tfreezing(kk) = Tfrz(S(i), par(i)%he, one/(par(i)%lambc*freezefac)) ! check for onset of freezing and adjust (increase) temperature to account for latent heat ! release by ice formation IF (experiment /= 184) THEN IF (Tsoil(i) < Tfreezing(kk) .AND. var(i)%iice==0) THEN ! start correction for onset of freezing dtdT(kk) = dthetalmaxdTh(Tfreezing(kk), S(i), par(i)%he, one/(par(i)%lambc*freezefac), & par(i)%thre,par(i)%the) tmp1d3(kk) = 1000._r_2 tmp1d2(kk) = Tsoil(i) - dTsoil(i) k = 1 IF ((i==1).AND.(h0(kk)- deltah0(kk))>zero) THEN h0(kk) = h0(kk) - deltah0(kk) ! calculate stored heat using h0 at t-dt DO WHILE ((k<nsteps_ice_max).AND.(tmp1d3(kk)>tol_dthetaldT)) tmp1d1(kk) = LHS_h(i)*dt(kk)-deltaJ_sensible_S(i) + tmp1d2(kk)* & (var(i)%csoil*dx(i)+cp(kk)*h0(kk)) + & Tfreezing(kk)*rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * & dtdT(kk)*(dx(i)+h0(kk)/par(i)%thre) tmp1d1(kk) = tmp1d1(kk) /( (var(i)%csoil*dx(i)+cp(kk)*h0(kk)) + & rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * & dtdT(kk)*(dx(i)+h0(kk)/par(i)%thre)) dtdT(kk) = (thetalmax(tmp1d1(kk), S(i), par(i)%he, one/(par(i)%lambc*freezefac), & par(i)%thre, par(i)%the) & - (S(i)*par(i)%thre + (par(i)%the-par(i)%thre)))/ & (tmp1d1(kk) - Tfreezing(kk)) tmp1d3(kk) = ABS((tmp1d1(kk)- Tfreezing(kk))*dtdT(kk)) k = k + 1 ENDDO ! end of do while dTsoil(i) = tmp1d1(kk) - tmp1d2(kk) Tsoil(i) = tmp1d1(kk) deltaJ_sensible_T(i) = (var(i)%csoil*dx(i)+ cp(kk)*h0(kk))*dTsoil(i) + & (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*(dx(1)+h0(kk)/par(1)%thre)* & ((cswat-csice)*(tmp1d2(kk))) deltaJ_latent_T(i) = (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*lambdaf* & (dx(1)+h0(kk)/par(1)%thre) h0(kk) = h0(kk) + deltah0(kk) ! change in heat stored in soil column dJcol_latent_S(kk) = SUM(deltaJ_latent_S(1:n)) dJcol_latent_T(kk) = SUM(deltaJ_latent_T(1:n)) dJcol_sensible(kk) = SUM(deltaJ_sensible_T(1:n)) + SUM(deltaJ_sensible_S(1:n)) ! tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + & ! dJcol_sensible(kk) - (qhsig(0)-qhsig(6))*dt(kk) + & ! qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk) tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + & dJcol_sensible(kk) - (qhsig(0)-qhsig(n))*dt(kk) + & qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk) ! if (abs(tmp1d1(kk)).gt.10.) then ! write(*,*) 'E imbalance check 2' ! endif ELSE DO WHILE ((k<nsteps_ice_max).AND.(tmp1d3(kk)>tol_dthetaldT)) ! zeroth estimate of corrected temperature tmp1d1(kk) = LHS_h(i)*dt(kk)-deltaJ_sensible_S(i) + tmp1d2(kk)* & (var(i)%csoil*dx(i)) + Tfreezing(kk)*rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * & dtdT(kk)*dx(i) tmp1d1(kk) = tmp1d1(kk) /( (var(i)%csoil*dx(i)) + & rhow*(lambdaf+(tmp1d2(kk))*(cswat-csice)) * dtdT(kk)*dx(i)) dtdT(kk) = (thetalmax(tmp1d1(kk), S(i), par(i)%he, & one/(par(i)%lambc*freezefac), par(i)%thre, par(i)%the) & - (S(i)*par(i)%thre + (par(i)%the-par(i)%thre)))/ & (tmp1d1(kk) - Tfreezing(kk)) tmp1d3(kk) = ABS((tmp1d1(kk)- Tfreezing(kk))*dtdT(kk)) k = k + 1 ENDDO ! end of do while dTsoil(i) = tmp1d1(kk)- (Tsoil(i)-dTsoil(i)) Tsoil(i) = tmp1d1(kk) deltaJ_sensible_T(i) = (var(i)%csoil*dx(i))*dTsoil(i) + & (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*dx(i)* & ((cswat-csice)*(tmp1d2(kk))) deltaJ_latent_T(i) = (Tsoil(i)-Tfreezing(kk))*rhow*dtdT(kk)*lambdaf*dx(i) ENDIF ENDIF ! end correction for onset of freezing ENDIF IF (var(i)%iice==1) THEN theta = S(i)*(par(i)%thre) + (par(i)%the - par(i)%thre) tmp1d1(kk) = Tsoil(i) - dTsoil(i) ! energy for complete melting IF ((i==1).AND.(h0(kk)>zero)) THEN tmp1d3(kk) = thetai(i)*dx(i)*rhow*lambdaf + & thetai(i)*(h0(kk)- deltah0(kk))/par(i)%thre*rhow*lambdaf tmp1d2(kk) = (Tfreezing(kk))*(rhow*cswat*(theta*dx(1)+h0(kk))+ & par(1)%rho*par(1)%css*dx(1)) - & (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(1)+(h0(kk)- deltah0(kk))* & (one-thetai(1)/par(1)%thre))+par(1)%rho*par(1)%css*dx(1)) - & (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(1)+(h0(kk)- deltah0(kk))* & thetai(1)/par(1)%thre)) ELSE tmp1d3(kk) = thetai(i)*dx(i)*rhow*lambdaf tmp1d2(kk) = (Tfreezing(kk))*(rhow*cswat*(theta*dx(i))+par(i)%rho*par(i)%css*dx(i)) - & (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(i))+par(i)%rho*par(i)%css*dx(i)) - & (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(i))) ENDIF ENDIF ! check for onset of thawing and decrease temperature to account for latent heat required to melt ice IF ((var(i)%iice==1).AND.((J0(i) + LHS_h(i)*dt(kk)).GE.JSoilLayer(Tfreezing(kk), & dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, & par(i)%he, one/(par(i)%lambc*freezefac)))) THEN ! start correction for onset of thawing ! Correct for "overthawing". This extra energy is used to heat the soil even more ! The correction comes from equating the energy balances before and after correction. tmp1d1(kk) = Tsoil(i) - dTsoil(i) IF ((i==1) .AND. (h0(kk)>zero)) THEN deltaJ_latent_T(i)= tmp1d3(kk) deltaJ_latent_S(i) = zero h0(kk) = h0(kk) - deltah0(kk) deltaJ_sensible_S(i) = (tmp1d1(kk))*(rhow*cswat*(theta*dx(1)+(h0(kk)+deltah0(kk)))+ & par(1)%rho*par(1)%css*dx(1)) - & (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(1)+(h0(kk))* & (one-thetai(1)/par(1)%thre))+par(1)%rho*par(1)%css*dx(1)) - & (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(1)+(h0(kk))* & thetai(1)/par(1)%thre)) dTsoil(i) = (LHS_h(i)*dt(kk) - (deltaJ_latent_T(i) + deltaJ_sensible_S(i)))/ & (rhow*cswat*(theta*dx(1)+(h0(kk)+deltah0(kk)))+par(1)%rho*par(1)%css*dx(1)) deltaJ_sensible_T(i) = dTsoil(i)* & (rhow*cswat*(theta*dx(1)+(h0(kk)+deltah0(kk)))+par(1)%rho*par(1)%css*dx(1)) h0(kk) = h0(kk) + deltah0(kk) ELSE deltaJ_latent_T(i) = thetai(i)*dx(i)*rhow*lambdaf deltaJ_latent_S(i) = zero theta = S(i)*(par(i)%thre) + (par(i)%the - par(i)%thre) var(i)%csoil = cswat*theta*rhow + par(i)%css*par(i)%rho deltaJ_sensible_S(i) = (tmp1d1(kk))*(rhow*cswat*(theta*dx(i))+ & par(i)%rho*par(i)%css*dx(i)) - & (tmp1d1(kk))*(rhow*cswat*(var(i)%thetal*dx(i))+par(i)%rho*par(i)%css*dx(i)) - & (tmp1d1(kk))*(rhow*csice*(thetai(i)*dx(i))) dTsoil(i) = (LHS_h(i)*dt(kk) - (deltaJ_latent_T(i) + deltaJ_sensible_S(i)))/ & (dx(i)*var(i)%csoil) deltaJ_sensible_T(i) = var(i)%csoil*dTsoil(i)*dx(i) ENDIF Tsoil(i) = tmp1d1(kk) + dTsoil(i) ! correct thetai and T in frozen soil ELSEIF ((var(i)%iice==1).AND.((J0(i) + LHS_h(i)*dt(kk))<JSoilLayer(Tfreezing(kk), & dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, & par(i)%he, one/(par(i)%lambc*freezefac)))) THEN tmp1d2(kk) = J0(i) + LHS_h(i)*dt(kk) ! total energy in soil layer IF (Tsoil(i) .LT. -40.0_r_2) THEN ! assume no lquid below min temp threshhold var(i)%thetal = 0.0_r_2 var(i)%thetai = theta IF (i.EQ.1) hice(kk) = h0(kk) tmp1d3(kk) = (tmp1d2(kk) + rhow*lambdaf*(theta*dx(i) + MERGE(h0(kk),zero,i==1)))/ & (dx(i)*par(i)%css*par(i)%rho + rhow*csice*(theta*dx(i) + & MERGE(h0(kk),zero,i==1))) ELSE !check there is a zero tmp1 = GTfrozen(Tsoil(i)-dTsoil(i)-50._r_2, tmp1d2(kk), dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, par(i)%he, one/(par(i)%lambc*freezefac)) tmp2 = GTFrozen(Tfreezing(kk), tmp1d2(kk), dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, par(i)%he, one/(par(i)%lambc*freezefac)) ! there is a zero in between IF ((tmp1*tmp2) < zero) THEN tmp1d3(kk) = rtbis_Tfrozen(tmp1d2(kk), dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, & par(i)%he, one/(par(i)%lambc*freezefac), & Tsoil(i)-dTsoil(i)-50._r_2, Tfreezing(kk), 0.0001_r_2) tmp1d4(kk) = thetalmax(tmp1d3(kk), S(i), par(i)%he, one/(par(i)%lambc*freezefac), & par(i)%thre, par(i)%the) ! liquid content at solution for Tsoil ELSE WRITE(logn,*) "Found no solution for Tfrozen 1. ", kk, i WRITE(logn,*) "Assume soil is totally frozen" var(i)%thetal = 0.0_r_2 var(i)%thetai = theta IF (i.EQ.1) hice(kk) = h0(kk) tmp1d3(kk) = (tmp1d2(kk) + rhow*lambdaf*(theta*dx(i) + MERGE(h0(kk),zero,i==1)))/ & (dx(i)*par(i)%css*par(i)%rho + rhow*csice*(theta*dx(i) + & MERGE(h0(kk),zero,i==1))) WRITE(logn,*) "frozen soil temperature: ", tmp1d3(kk) WRITE(logn,*) nsteps(kk), S(i), Tsoil(i), dTsoil(i), h0(kk), tmp1, tmp2, tmp1d2(kk), theta, & JSoilLayer(Tfreezing(kk), & dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, & par(i)%he, one/(par(i)%lambc*freezefac)), J0(i) + LHS_h(i)*dt(kk), Tfreezing(kk) ENDIF var(i)%thetal = tmp1d4(kk) var(i)%thetai = theta - tmp1d4(kk) ENDIF ! Tsoil <-40 IF (i==1) THEN hice_tmp(kk) = hice(kk) hice(kk) = h0(kk)*var(1)%thetai/par(1)%thre ! correct total energy stored in pond + soil deltaJ_latent_S(i) = - rhow*lambdaf*((hice(kk)-hice_tmp(kk)) + & dx(1)*(var(1)%thetai-thetai(1))) ELSE ! correct total energy stored in soil deltaJ_latent_S(i) = - rhow*lambdaf*( + & dx(i)*(var(i)%thetai-thetai(i))) ENDIF deltaJ_latent_T(i) = zero deltaJ_sensible_T(i) = JSoilLayer(tmp1d3(kk), & dx(i), theta,par(i)%css, par(i)%rho, & MERGE(h0(kk),zero,i==1), par(i)%thre, par(i)%the, & par(i)%he, one/(par(i)%lambc*freezefac)) - & J0(i)- & deltaJ_latent_S(i) deltaJ_sensible_S(i) = 0._r_2 thetai(i) = var(i)%thetai Tsoil(i) = tmp1d3(kk) ENDIF ! soil remains frozen ENDIF ! if .not.again END DO ! i=1, n => update variables S,T END SUBROUTINE update_s_t SUBROUTINE update_unknowns( & irec, mp, qprec, qprec_snow, n, dx, h0, S, & Tsoil, evap, infil, drainage, discharge, & qh, nsteps, vmet, vlit, vsnow, var, Hcum, lEcum, & Gcum, Qadvcum, & csoil, kth, phi, zdelta, SL, Tl, & par, wex, & qvsig, qlsig, qvTsig, & precip, qevap, rexcol, wcol, & again, ih0, ns, & dt, dwinfil, dwoff, & sig, & qexd, dy, & de, q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, & qsig, qhsig, qadvsig, qliq, qv, qvT, qlyb, & qvya, qvyb, qlTb, qvTa, qvTb, dTsoil, & v_aquifer, & dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge, & dJcol_latent_S, dJcol_latent_T, dJcol_sensible, deltaJ_latent_S, & deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T, qevapsig, & qrunoff, tmp1d1, tmp1d2, tmp1d3, deltah0, deltaSL, & LHS_h, surface_case, & litter, i, j, kk, & advection, theta, Tqw, cp, & hsnow, & delta_snowcol, delta_snowT, & delta_snowliq, J0, iqex, & J0snow, wcol0snow & ) IMPLICIT NONE INTEGER(i_d) :: irec, mp REAL(r_2), DIMENSION(1:mp) :: qprec REAL(r_2), DIMENSION(1:mp) :: qprec_snow INTEGER(i_d) :: n REAL(r_2), DIMENSION(1:n) :: dx REAL(r_2), DIMENSION(1:mp) :: h0 REAL(r_2), DIMENSION(1:n) :: S REAL(r_2), DIMENSION(1:n) :: Tsoil REAL(r_2), DIMENSION(1:mp) :: evap REAL(r_2), DIMENSION(1:mp) :: infil REAL(r_2), DIMENSION(1:mp) :: drainage, discharge REAL(r_2), DIMENSION(-nsnow_max:n) :: qh INTEGER(i_d), DIMENSION(1:mp) :: nsteps TYPE(vars_met), DIMENSION(1:mp) :: vmet TYPE(vars), DIMENSION(1:mp) :: vlit TYPE(vars_snow), DIMENSION(1:mp) :: vsnow TYPE(vars), DIMENSION(1:n) :: var REAL(r_2), DIMENSION(1:mp) :: Hcum, lEcum REAL(r_2), DIMENSION(1:mp) :: Gcum, Qadvcum REAL(r_2), DIMENSION(1:n) :: csoil, kth REAL(r_2), DIMENSION(1:n) :: phi REAL(r_2), DIMENSION(1:mp) :: zdelta, SL, Tl TYPE(params), DIMENSION(1:n) :: par REAL(r_2), DIMENSION(1:mp,1:n), OPTIONAL :: wex REAL(r_2), DIMENSION(-nsnow_max:n) :: qvsig, qlsig, qvTsig REAL(r_2), DIMENSION(1:mp) :: precip, qevap REAL(r_2), DIMENSION(1:mp) :: rexcol, wcol LOGICAL, DIMENSION(1:mp) :: again INTEGER(i_d), DIMENSION(1:mp) :: ih0, ns REAL(r_2), DIMENSION(1:mp) :: dt, dwinfil, dwoff REAL(r_2), DIMENSION(1:mp) :: sig REAL(r_2), DIMENSION(1:n) :: qexd REAL(r_2), DIMENSION(-nsnow_max+1:n) :: dy REAL(r_2), DIMENSION(-nsnow_max+1:n) :: de REAL(r_2), DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb REAL(r_2), DIMENSION(-nsnow_max:n) :: qsig, qhsig, qadvsig REAL(r_2), DIMENSION(-nsnow_max:n) :: qliq, qv, qvT, qlyb, qvya, qvyb, qlTb, qvTa, qvTb REAL(r_2), DIMENSION(1:n) :: dTsoil TYPE(vars_aquifer), DIMENSION(1:mp) :: v_aquifer REAL(r_2), DIMENSION(1:mp) :: dwcol, dwdrainage, drn,inlit, dwinlit, drexcol, dwdischarge REAL(r_2), DIMENSION(1:mp) :: dJcol_latent_S, dJcol_latent_T, dJcol_sensible REAL(r_2), DIMENSION(1:n) :: deltaJ_latent_S, deltaJ_latent_T, deltaJ_sensible_S, deltaJ_sensible_T REAL(r_2), DIMENSION(1:mp) :: qevapsig REAL(r_2), DIMENSION(1:mp) :: qrunoff REAL(r_2), DIMENSION(1:mp) :: tmp1d1, tmp1d2, tmp1d3 REAL(r_2), DIMENSION(1:mp) :: deltah0 REAL(r_2), DIMENSION(1:mp) :: deltaSL REAL(r_2), DIMENSION(-nsnow_max+1:n) :: LHS_h INTEGER(i_d), DIMENSION(1:mp) :: surface_case LOGICAL :: litter INTEGER(i_d) :: i, j, kk INTEGER(i_d) :: advection ! switches REAL(r_2) :: theta REAL(r_2) :: Tqw REAL(r_2), DIMENSION(1:mp) :: cp REAL(r_2), DIMENSION(nsnow_max) :: hsnow REAL(r_2), DIMENSION(nsnow_max) :: delta_snowcol, delta_snowT, delta_snowliq REAL(r_2), DIMENSION(1:n) :: J0 REAL(r_2), DIMENSION(1:n) :: iqex REAL(r_2), DIMENSION(1:mp) :: J0snow, wcol0snow !----- update unknowns ! i.e. update state variables to end of time step ! cumulate some surface fluxes ih0(kk) = 0 IF (.NOT. again(kk)) THEN ! get surface fluxes at sigma of the time step dwoff(kk) = zero deltah0(kk) = zero precip(kk) = precip(kk) + (qprec(kk)+qprec_snow(kk))*dt(kk) dwcol(kk) = SUM(par(1:n)%thre*dx(1:n)*dy(1:n)*(ABS(var(1:n)%isat-1)),1) + & REAL(1-ns(kk),r_2)*dy(1) IF (vsnow(kk)%nsnow==1) dwcol(kk) = dwcol(kk) + dy(0) ! absolute heat stored in soil column before update DO j=1, n theta = S(j)*(par(j)%thre) + (par(j)%the - par(j)%thre) J0(j) = JSoilLayer(Tsoil(j), & dx(j), theta,par(j)%css, par(j)%rho, & MERGE(h0(kk),zero,j.EQ.1), par(j)%thre, par(j)%the, & par(j)%he, one/(par(j)%lambc*freezefac)) END DO ! change in heat stored in soil column DO j=1, n IF (j==1) THEN ! pond included in top soil layer deltaJ_latent_S(j) = -dx(j)*dy(j)*par(j)%thre*REAL(var(j)%iice,r_2)* & REAL(1-var(j)%isat,r_2)*rhow*lambdaf - & REAL(1-ns(kk),r_2)*dy(1)*REAL(var(j)%iice,r_2)*rhow*lambdaf*(var(1)%thetai/par(1)%thre) deltaJ_latent_T(j) = var(j)%dthetaldT*dx(j)*dTsoil(j)*rhow*lambdaf* & REAL(var(j)%iice,r_2) + & h0(kk)/par(1)%thre*var(j)%dthetaldT*dTsoil(j)*rhow*lambdaf*REAL(var(j)%iice,r_2) deltaJ_sensible_T(j) = var(j)%csoil*dx(j)*dTsoil(j) + & cp(kk)*h0(kk)*dTsoil(j) + & var(1)%iice*var(j)%dthetaldT*(Tsoil(j))*rhow*(cswat-csice)*dx(j)*dTsoil(j) + & var(j)%iice*var(j)%dthetaldT*(Tsoil(j))*rhow*(cswat-csice)* & h0(kk)/par(j)%thre*dTsoil(j) IF (advection==1) THEN deltaJ_sensible_S(j) = rhow*cswat*(Tsoil(j))*dx(j)*par(j)%thre*dy(j)* & REAL(1-var(j)%isat,r_2)*REAL(1-var(j)%iice,r_2) + & rhow*csice*(Tsoil(j))*dx(j)*par(j)%thre*dy(j)* & REAL(1-var(j)%isat,r_2)*REAL(var(j)%iice,r_2) + & REAL(1-ns(kk),r_2)*REAL(1-var(1)%iice,r_2)* & dy(1)*rhow*cswat*(Tsoil(1)) + & REAL(1-ns(kk),r_2)*REAL(var(1)%iice,r_2)* & dy(1)*rhow*csice*(Tsoil(1))*(var(1)%thetai/par(1)%thre) + & REAL(1-ns(kk),r_2)*REAL(var(1)%iice,r_2)* & dy(1)*rhow*cswat*(Tsoil(1))*(one-(var(1)%thetai/par(1)%thre)) ENDIF ELSE ! (j>1) deltaJ_latent_S(j) = -dx(j)*dy(j)*par(j)%thre*REAL(var(j)%iice,r_2)* & REAL(1-var(j)%isat,r_2)*rhow*lambdaf deltaJ_sensible_T(j) = var(j)%csoil*dx(j)*dTsoil(j) + & var(j)%iice*var(j)%dthetaldT*(Tsoil(j))*rhow*(cswat-csice)*dx(j)*dTsoil(j) IF (advection==1) THEN deltaJ_sensible_S(j) = rhow*cswat*(Tsoil(j))*dx(j)*dy(j)*par(j)%thre* & REAL(1-var(j)%isat,r_2)*REAL(1-var(j)%iice,r_2) & + rhow*csice*(Tsoil(j))*dx(j)*dy(j)*par(j)%thre & *REAL(1-var(j)%isat,r_2)*REAL(var(j)%iice,r_2) ENDIF deltaJ_latent_T(j) = var(j)%dthetaldT*dx(j)*dTsoil(j)*rhow*lambdaf*REAL(var(j)%iice,r_2) ENDIF END DO ! change in heat stored in soil column dJcol_latent_S(kk) = SUM(deltaJ_latent_S(1:n)) dJcol_latent_T(kk) = SUM(deltaJ_latent_T(1:n)) dJcol_sensible(kk) = SUM(deltaJ_sensible_T(1:n)) + SUM(deltaJ_sensible_S(1:n)) ! tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + dJcol_sensible(kk) - (qhsig(0)-qhsig(6))*dt(kk) & ! + qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk) tmp1d1(kk) = dJcol_latent_S(kk) + dJcol_latent_T(kk) + dJcol_sensible(kk) - (qhsig(0)-qhsig(n))*dt(kk) & + qrunoff(kk)*cswat*rhow*(Tsoil(1) + sig(kk)*dTsoil(1))*dt(kk) IF (ns(kk)<1) THEN ! change in pond height h0(kk) = h0(kk) + dy(1) deltah0(kk) = dy(1) IF (h0(kk)<zero .AND. dy(1)<zero) THEN ih0(kk)=1 ! pond gone ENDIF ENDIF ! cumulate evaporation from top of soil column or top of litter/pond or top of snow pack SELECT CASE (surface_case(kk)) CASE(1) ! no snow evap(kk) = evap(kk) +qevap(kk)*dt(kk) - sig(kk)*(qyb(0)*dy(1)+qTb(0)*dTsoil(1))*dt(kk) qevapsig(kk) = qevap(kk) - sig(kk)*(qyb(0)*dy(1)+qTb(0)*dTsoil(1)) Gcum(kk) = Gcum(kk)+(qhsig(0)-qadvsig(0))*dt(kk) Qadvcum(kk) = Qadvcum(kk) + qadvsig(0)*dt(kk) - qadvsig(n)*dt(kk) lEcum(kk) = lEcum(kk) + qevapsig(kk)*thousand*var(1)%lambdav*dt(kk) Hcum(kk) = Hcum(kk) + (vmet(kk)%Rn*dt(kk)-(qhsig(0)-qadvsig(0))*dt(kk)- & qevapsig(kk)*thousand*var(1)%lambdav*dt(kk)) dwinfil(kk) = (qprec(kk)+qprec_snow(kk)-qevapsig(kk))*dt(kk) CASE(2) ! dedicated snow layer evap(kk) = evap(kk) +qevap(kk)*dt(kk) - sig(kk)*(qyb(-vsnow(kk)%nsnow)*dy(1-vsnow(kk)%nsnow)+ & qTb(-vsnow(kk)%nsnow)*de(1-vsnow(kk)%nsnow))*dt(kk) qevapsig(kk) = qevap(kk) - sig(kk)*(qyb(-vsnow(kk)%nsnow)*dy(1-vsnow(kk)%nsnow)+ & qTb(-vsnow(kk)%nsnow)*de(1-vsnow(kk)%nsnow)) Gcum(kk) = Gcum(kk)+(qhsig(-vsnow(kk)%nsnow)-qadvsig(-vsnow(kk)%nsnow))*dt(kk) Qadvcum(kk) = Qadvcum(kk) + qadvsig(-vsnow(kk)%nsnow)*dt(kk) - qadvsig(n)*dt(kk) IF (vsnow(kk)%hliq(1).GT.zero) THEN lEcum(kk) = lEcum(kk) + qevapsig(kk)*thousand*lambdaf*dt(kk) Hcum(kk) = Hcum(kk) + (vmet(kk)%Rn*dt(kk)-(qhsig(-vsnow(kk)%nsnow)-qadvsig(-vsnow(kk)%nsnow))*dt(kk)- & qevapsig(kk)*thousand*lambdaf*dt(kk)) ELSE lEcum(kk) = lEcum(kk) + qevapsig(kk)*thousand*lambdas*dt(kk) Hcum(kk) = Hcum(kk) + (vmet(kk)%Rn*dt(kk)-(qhsig(-vsnow(kk)%nsnow)-qadvsig(-vsnow(kk)%nsnow))*dt(kk)- & qevapsig(kk)*thousand*lambdas*dt(kk)) ENDIF dwinfil(kk) = (qprec(kk)+qprec_snow(kk)-qevapsig(kk))*dt(kk) END SELECT ! surface_case dwdrainage(kk) = q(n)*dt(kk) +sig(kk)*dt(kk)*(qya(n)*dy(n)+qTa(n)*dTsoil(n)) IF (botbc=="aquifer" .AND. v_aquifer(kk)%isat==0) THEN dwdischarge(kk) = v_aquifer(kk)%discharge*dt(kk) ELSE dwdischarge(kk) = dwdrainage(kk) ENDIF drexcol(kk) = SUM(iqex(:),1)*dt(kk) ! cumulative fluxes infil(kk) = infil(kk) + dwinfil(kk) inlit(kk) = inlit(kk) + dwinlit(kk) wcol(kk) = wcol(kk) + dwcol(kk) drainage(kk) = drainage(kk) + dwdrainage(kk) discharge(kk) = discharge(kk) + qrunoff(kk)*dt(kk) rexcol(kk) = rexcol(kk) + drexcol(kk) Qadvcum(kk) = Qadvcum(kk) - SUM(iqex(:)*(Tsoil(:)),1)*dt(kk)*rhow*cswat - & qrunoff(kk)*(Tsoil(1))*dt(kk)*rhow*cswat ! if (irec.eq.109) then ! write(*,"(1i8,16e16.6)") irec, infil(kk)-(wcol(kk)+discharge(kk)+drainage(kk))-rexcol(kk), & ! wcol(kk), deltah0(kk), h0(kk), discharge(kk), evap(kk), & ! drainage(kk), rexcol(kk) ! write(*,"(1i8,16e16.6)") irec, dwinfil(kk)-dwcol(kk)-dwdrainage(kk)-drexcol(kk)-qrunoff(kk)*dt(kk) ! endif ! evaluate soil fluxes at sigma of time step for use in isotope routine SELECT CASE (surface_case(kk)) CASE(1) ! no snow qsig(0) = q(0) + sig(kk)*(qyb(0)*dy(1) + qya(0)*dy(0) + qTb(0)*de(1) & + qTa(0)*de(0)) qhsig(0) = qh(0) + sig(kk)*(qhyb(0)*dy(1) + qhya(0)*dy(0) + qhTb(0)*de(1) & + qhTa(0)*de(0)) qvsig(0) = qv(0)+sig(kk)*qvyb(0)*dy(1)+sig(kk)*qvTb(0)*de(1) qlsig(0) = qliq(0)+sig(kk)*qlyb(0)*dy(1)+sig(kk)*qlTb(0)*de(1) qsig(1:n-1) = q(1:n-1) + sig(kk)*(qya(1:n-1)*dy(1:n-1) + qyb(1:n-1)*dy(2:n) & + qTa(1:n-1)*de(1:n-1) + qTb(1:n-1)*de(2:n)) qsig(n) = q(n) + sig(kk)*(qya(n)*dy(n) + qTa(n)*de(n)) qhsig(1:n-1) = qh(1:n-1) + sig(kk)*(qhya(1:n-1)*dy(1:n-1) + qhyb(1:n-1)*dy(2:n) & + qhTa(1:n-1)*de(1:n-1) + qhTb(1:n-1)*de(2:n)) qhsig(n) = qh(n) + sig(kk)*(qhya(n)*dy(n) + qhTa(n)*de(n)) qvsig(1:n-1) = qv(1:n-1) + sig(kk)*(qvya(1:n-1)*dy(1:n-1) + qvyb(1:n-1)*dy(2:n) & + qvTa(1:n-1)*de(1:n-1) + qvTb(1:n-1)*de(2:n)) qvsig(n) = zero qvTsig(0) = qvT(0) + sig(kk)*qvTb(0)*de(1) qvTsig(1:n-1) = qvT(1:n-1) + sig(kk)*(qvTa(1:n-1)*de(1:n-1) + qvTb(1:n-1)*de(2:n)) qvTsig(n) = zero qlsig(1:n-1) = qsig(1:n-1) - qvsig(1:n-1) qlsig(n) = qsig(n) CASE(2) ! dedicated snow layer qsig(-vsnow(kk)%nsnow) = q(-vsnow(kk)%nsnow) + & sig(kk)*(qyb(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow+1) + & ! qya(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow) + & qTb(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow+1)) ! & ! + qTa(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow)) qhsig(-vsnow(kk)%nsnow) = qh(-vsnow(kk)%nsnow) + & sig(kk)*(qhyb(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow+1) + & ! qhya(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow) + & qhTb(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow+1)) !& ! + qhTa(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow)) qvsig(-vsnow(kk)%nsnow) = qv(-vsnow(kk)%nsnow) + & sig(kk)*qvyb(-vsnow(kk)%nsnow)*dy(-vsnow(kk)%nsnow+1) + & sig(kk)*qvTb(-vsnow(kk)%nsnow)*de(-vsnow(kk)%nsnow+1) qlsig(-vsnow(kk)%nsnow:-1) = zero qlsig(0) = qliq(0)+sig(kk)*qlyb(0)*dy(1)+sig(kk)*qlTb(0)*de(1) qsig(-vsnow(kk)%nsnow+1:n-1) = q(-vsnow(kk)%nsnow+1:n-1) + & sig(kk)*(qya(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+1:n-1) + & qyb(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+2:n) & + qTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + & qTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n)) qsig(n) = q(n) + sig(kk)*(qya(n)*dy(n) + qTa(n)*de(n)) qhsig(-vsnow(kk)%nsnow+1:n-1) = qh(-vsnow(kk)%nsnow+1:n-1) + & sig(kk)*(qhya(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+1:n-1) + & qhyb(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+2:n) & + qhTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + & qhTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n)) qhsig(n) = qh(n) + sig(kk)*(qhya(n)*dy(n) + qhTa(n)*de(n)) qvsig(-vsnow(kk)%nsnow+1:n-1) = qv(-vsnow(kk)%nsnow+1:n-1) + & sig(kk)*(qvya(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+1:n-1) + & qvyb(-vsnow(kk)%nsnow+1:n-1)*dy(-vsnow(kk)%nsnow+2:n) & + qvTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + & qvTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n)) qvsig(n) = zero qvTsig(0) = qvT(0) + sig(kk)*qvTb(0)*de(1) qvTsig(-vsnow(kk)%nsnow+1:n-1) = qvT(-vsnow(kk)%nsnow+1:n-1) + & sig(kk)*(qvTa(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+1:n-1) + & qvTb(-vsnow(kk)%nsnow+1:n-1)*de(-vsnow(kk)%nsnow+2:n)) qvTsig(n) = zero qlsig(-vsnow(kk)%nsnow+1:n-1) = qsig(-vsnow(kk)%nsnow+1:n-1) - qvsig(-vsnow(kk)%nsnow+1:n-1) qlsig(n) = qsig(n) END SELECT ! surface_case IF (botbc=="aquifer") THEN ! update aquifer props IF (v_aquifer(kk)%isat==0) THEN IF (v_aquifer(kk)%WA+(dwdrainage(kk)-dwdischarge(kk))*dt(kk) > & (v_aquifer(kk)%zzero-v_aquifer(kk)%zsoil)*v_aquifer(kk)%Sy) THEN v_aquifer(kk)%isat = 1 ! aquifer saturated S(n) = S(n) + (-(v_aquifer(kk)%WA+(dwdrainage(kk)-dwdischarge(kk))*dt(kk)) & +(v_aquifer(kk)%zzero-v_aquifer(kk)%zsoil)*v_aquifer(kk)%Sy)/(dx(n)*par(n)%thre) ENDIF v_aquifer(kk)%WA = MIN(v_aquifer(kk)%WA+(dwdrainage(kk)-dwdischarge(kk))*dt(kk), & (v_aquifer(kk)%zzero-v_aquifer(kk)%zsoil)*v_aquifer(kk)%Sy) v_aquifer(kk)%zdelta = v_aquifer(kk)%zzero - v_aquifer(kk)%Wa/v_aquifer(kk)%Sy ! new discharge rate v_aquifer(kk)%discharge = v_aquifer(kk)%Rsmax*EXP(-v_aquifer(kk)%f*v_aquifer(kk)%zdelta) ELSEIF (v_aquifer(kk)%isat==1) THEN ! check for desat of aquifer IF (dwdrainage(kk) < v_aquifer(kk)%Rsmax*EXP(-v_aquifer(kk)%f*v_aquifer(kk)%zsoil)) THEN v_aquifer(kk)%isat = 0 v_aquifer(kk)%zdelta = v_aquifer(kk)%zsoil ! new discharge rate v_aquifer(kk)%discharge = v_aquifer(kk)%Rsmax*EXP(-v_aquifer(kk)%f*v_aquifer(kk)%zdelta) ENDIF ENDIF zdelta(kk) = v_aquifer(kk)%zdelta ENDIF ! end update aquifer props IF (botbc=="constant head") THEN drn(kk) = drn(kk)+(q(n)+sig(kk)*qya(n)*dy(n))*dt(kk) ELSE drn(kk) = drn(kk)+(q(n)+sig(kk)*qya(n)*dy(n))*dt(kk) END IF IF (PRESENT(wex)) THEN wex(kk,1:n) = wex(kk,1:n)+(iqex(1:n)+SPREAD(sig(kk),1,n)*qexd(1:n)*dy(1:n))*SPREAD(dt(kk),1,n) END IF IF (litter .AND. ns(kk)==1 ) THEN ! adjust litter moisture content if no ponding SL(kk) = SL(kk) + dy(0) deltaSL(kk) = dy(0) ELSE deltaSL(kk) = zero ENDIF TL(kk) = TL(kk) + de(0) ! Tlitter assumed the same as pond T IF (SL(kk) >= one) vlit(kk)%isat = 1 ! check for litter saturation csoil(:) = var(1:n)%csoil kth(:) = var(1:n)%kth phi(:) = var(1:n)%phi ! update snow pack water content and temperature (or liquid water) IF (vsnow(kk)%nsnow>0) THEN vsnow(kk)%Qadv_snow = vsnow(kk)%Qadv_snow + rhow*(qprec_snow(kk))* & (csice*(MIN(vmet(kk)%Ta,zero))-lambdaf)*dt(kk) vsnow(kk)%Qadv_rain = vsnow(kk)%Qadv_rain + rhow*(qprec(kk))*cswat*(MAX(vmet(kk)%Ta,zero))*dt(kk) ! update heat flux components Tqw = MERGE(vmet(kk)%Ta, vsnow(kk)%tsn(vsnow(kk)%nsnow), -qevap(kk)>zero) !vsnow(kk)%Qadv_vap = vsnow(kk)%Qadv_vap + & ! rhow*(-qevapsig(kk))*cswat*Tqw*dt(kk) - & ! qadvsig(0)*dt(kk) vsnow(kk)%Qadv_vap =vsnow(kk)%Qadv_vap + (qadvsig(-vsnow(kk)%nsnow)-qadvsig(0))*dt(kk) - & rhow*(qprec_snow(kk))*(csice*(MIN(vmet(kk)%Ta,zero))-lambdaf)*dt(kk) - & rhow*(qprec(kk))*cswat*(vmet(kk)%Ta)*dt(kk) vsnow(kk)%Qcond_net = vsnow(kk)%Qcond_net + (qhsig(-vsnow(kk)%nsnow) - qhsig(0))*dt(kk) - & (qadvsig(-vsnow(kk)%nsnow) - qadvsig(0))*dt(kk) DO i=1, vsnow(kk)%nsnow IF (vsnow(kk)%hliq(i)>zero) THEN IF ((vsnow(kk)%hliq(i) + de(i-vsnow(kk)%nsnow))>(dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i))) THEN ! account here for new hliq exceeding new hsnow tmp1d1(kk) = vsnow(kk)%hliq(i)*cswat*rhow*zero + & (vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow*(zero*csice-lambdaf) tmp1d2(kk) = ((tmp1d1(kk) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk)))/ & ((dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i))*cswat*rhow) vsnow(kk)%tsn(i) = tmp1d2(kk) IF (vsnow(kk)%tsn(i)>zero) THEN WRITE(*,*) 'T>zero ', irec, nsteps(kk), i, vsnow(kk)%tsn(i) ENDIF vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) + & lambdaf*(vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk) - & lambdaf*(vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow vsnow(kk)%hliq(i) = vsnow(kk)%hsnow(i) + dy(i-vsnow(kk)%nsnow) ELSEIF ((vsnow(kk)%hliq(i)+de(i-vsnow(kk)%nsnow))<zero) THEN ! account here for new hliq less than zero tmp1d1(kk) = (csice*(vsnow(kk)%tsn(i))-lambdaf)*rhow*(vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i)) + & cswat*(vsnow(kk)%tsn(i))*rhow*vsnow(kk)%hliq(i) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk) delta_snowT(i) = (tmp1d1(kk)/rhow/(dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i))+lambdaf)/csice vsnow(kk)%tsn(i) = delta_snowT(i) vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) - lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + & lambdaf*vsnow(kk)%hliq(i)*rhow vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + LHS_h(i-vsnow(kk)%nsnow)*dt(kk) - & (- lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + lambdaf*vsnow(kk)%hliq(i)*rhow) vsnow(kk)%hliq(i) = zero ELSE vsnow(kk)%hliq(i) = vsnow(kk)%hliq(i) + de(i-vsnow(kk)%nsnow) delta_snowliq(i) = de(i-vsnow(kk)%nsnow) delta_snowT(i) = zero vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) - & lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + lambdaf*de(i-vsnow(kk)%nsnow)*rhow vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + & dy(i-vsnow(kk)%nsnow)*rhow*csice*(vsnow(kk)%tsn(i)) + & de(i-vsnow(kk)%nsnow)*rhow*(cswat-csice)*(vsnow(kk)%tsn(i)) ENDIF ELSE !hliq = zero IF ((vsnow(kk)%tsn(i) + de(i-vsnow(kk)%nsnow))<zero) THEN vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) - lambdaf*dy(i-vsnow(kk)%nsnow)*rhow vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) + & dy(i-vsnow(kk)%nsnow)*rhow*csice*(vsnow(kk)%tsn(i)) + & de(i-vsnow(kk)%nsnow)*csice*rhow*hsnow(i) vsnow(kk)%tsn(i) = vsnow(kk)%tsn(i) + de(i-vsnow(kk)%nsnow) delta_snowT(i) = de(i-vsnow(kk)%nsnow) delta_snowliq(i) = zero vsnow(kk)%hliq(i) = zero ELSE ! vsnow(kk)%tsn + de(0))>zero delta_snowT(i) = zero - vsnow(kk)%tsn(i) delta_snowliq(i) = (LHS_h(i-vsnow(kk)%nsnow)*dt(kk) - & delta_snowT(i)*csice*hsnow(i) *rhow - & dy(i-vsnow(kk)%nsnow)*rhow*(csice*(vsnow(kk)%tsn(i))-lambdaf))/ & (rhow*(zero*(cswat-csice)+lambdaf)) vsnow(kk)%hliq(i) = vsnow(kk)%hliq(i) + delta_snowliq(i) vsnow(kk)%deltaJlatent(i) = vsnow(kk)%deltaJlatent(i) & - lambdaf*dy(i-vsnow(kk)%nsnow)*rhow + & lambdaf*rhow*delta_snowliq(i) vsnow(kk)%deltaJsensible(i) = vsnow(kk)%deltaJsensible(i) & + dy(i-vsnow(kk)%nsnow)*rhow*csice*(vsnow(kk)%tsn(i)) + & delta_snowT(i)*csice*rhow*hsnow(i) vsnow(kk)%tsn(i) = vsnow(kk)%tsn(i) + delta_snowT(i) ENDIF ! vsnow(kk)%tsn + de(0))>zero ENDIF ! hliq = zero vsnow(kk)%hsnow(i) = dy(i-vsnow(kk)%nsnow) + vsnow(kk)%hsnow(i) vsnow(kk)%Jsensible(i) = (vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow*csice*(vsnow(kk)%Tsn(i)) vsnow(kk)%Jlatent(i) = (vsnow(kk)%hsnow(i)-vsnow(kk)%hliq(i))*rhow*(-lambdaf) delta_snowcol(i) = dy(i-vsnow(kk)%nsnow) ENDDO ! loop over snow layers ! total internal energy and water content of snowpack vsnow(kk)%wcol = SUM(vsnow(kk)%hsnow(1:vsnow(kk)%nsnow)) vsnow(kk)%deltawcol = vsnow(kk)%wcol - wcol0snow(kk) vsnow(kk)%Qprec = vsnow(kk)%Qprec + qprec_snow(kk)*dt(kk) + qprec(kk)*dt(kk) vsnow(kk)%Qevap = vsnow(kk)%Qevap + qevapsig(kk)*dt(kk) vsnow(kk)%Qvap = vsnow(kk)%Qvap + qsig(0)*dt(kk) vsnow(kk)%J = SUM(vsnow(kk)%Jsensible(1:vsnow(kk)%nsnow)+vsnow(kk)%Jlatent(1:vsnow(kk)%nsnow)) vsnow(kk)%deltaJ = vsnow(kk)%J - J0snow(kk) IF (vsnow(kk)%nsnow.LT.nsnow_max) THEN vsnow(kk)%hsnow(vsnow(kk)%nsnow+1:nsnow_max) = zero ENDIF ! update densities and depths of snow layers ! effect of new snowfall ! snowfall density (tmp1d1), LaChapelle 1969 IF (vmet(kk)%Ta > 2.0_r_2) THEN tmp1d1(kk) = 189.0_r_2 ELSEIF ((vmet(kk)%Ta > -15.0_r_2) .AND. (vmet(kk)%Ta <= 2.0_r_2)) THEN tmp1d1(kk) = 50.0_r_2 + 1.7_r_2*(vmet(kk)%Ta+15.0_r_2)**1.5_r_2 ELSE tmp1d1(kk) = 50.0_r_2 ENDIF IF ((vsnow(kk)%hsnow(1)-qprec_snow(kk)*dt(kk)) > 1.e-3_r_2) THEN vsnow(kk)%dens(1) = (vsnow(kk)%dens(1)*(vsnow(kk)%hsnow(1)-qprec_snow(kk)*dt(kk)) & + tmp1d1(kk)*qprec_snow(kk)*dt(kk))/vsnow(kk)%hsnow(1) ENDIF DO i=1, vsnow(kk)%nsnow vsnow(kk)%depth(i) = vsnow(kk)%hsnow(i)/(vsnow(kk)%dens(i)/rhow) ! adjust depth for compaction IF (vsnow(kk)%hliq(i) .GT. zero) THEN tmp1d1(kk) = MERGE(one, EXP(-0.06_r_2*(vsnow(kk)%dens(i)-150._r_2)), vsnow(kk)%dens(i)<150._r_2) tmp1d2(kk) = 2 ELSE tmp1d1(kk) = MERGE(one, EXP(-0.046_r_2*(vsnow(kk)%dens(i)-150._r_2)), vsnow(kk)%dens(i)<150._r_2) tmp1d2(kk) = 1 ENDIF tmp1d3(kk) = vsnow(kk)%hsnow(i)*rhow/2._r_2 ! overburden kg/m2 IF (i.GT.1) THEN tmp1d3(kk) = tmp1d3(kk) + SUM(vsnow(kk)%hsnow(1:i-1)*rhow) ENDIF IF (vsnow(kk)%depth(i)*dt(kk)* & (2.8e-6_r_2*tmp1d1(kk)*tmp1d2(kk)*EXP(-vsnow(kk)%tsn(i)/25.0_r_2) + & tmp1d3(kk)/3.6e6_r_2*EXP(-0.08_r_2*vsnow(kk)%tsn(i))*EXP(-0.023_r_2*vsnow(kk)%dens(i))) & .LT.0.5_r_2*vsnow(kk)%depth(i)) THEN vsnow(kk)%depth(i) = vsnow(kk)%depth(i) - vsnow(kk)%depth(i)*dt(kk)* & (2.8e-6_r_2*tmp1d1(kk)*tmp1d2(kk)*EXP(-vsnow(kk)%tsn(i)/25.0_r_2) + & ! metamorphism tmp1d3(kk)/3.6e6_r_2*EXP(-0.08_r_2*vsnow(kk)%tsn(i))*EXP(-0.023_r_2*vsnow(kk)%dens(i))) ! overburden ENDIF vsnow(kk)%dens(i) = vsnow(kk)%hsnow(i)*rhow/vsnow(kk)%depth(i) ENDDO ENDIF ! end update snowlayers vsnow(kk)%MoistureFluxDivergence = vsnow(kk)%Qprec - vsnow(kk)%Qevap - vsnow(kk)%Qvap & - vsnow(kk)%Qmelt + vsnow(kk)%Qtransfer vsnow(kk)%FluxDivergence = vsnow(kk)%Qadv_rain + vsnow(kk)%Qadv_snow + vsnow(kk)%Qadv_vap + & vsnow(kk)%Qadv_melt + vsnow(kk)%Qcond_net + vsnow(kk)%Qadv_transfer ENDIF ! .not. again END SUBROUTINE update_unknowns SUBROUTINE get_and_solve_eqn( & irec, mp, qprec, n, dx, h0, S, thetai, & Tsoil, & qh, nsteps, vlit, vsnow, var, & dxL, & plit, par, & deltaTa, & qevap, qL, qhL, qybL, qTbL, qhTbL, qhybL, & again, again_ice, iok, itmp, ns, & accel, dt, fac, & phip, rsig, rsigdt, sig, t, & qexd, aa, bb, cc, dd, ee, ff, gg, dy, aah, bbh, cch, ddh, eeh, ffh, ggh, & de, q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb, qadv, qadvya, qadvyb, & qadvTa, qadvTb, qsig, qhsig, qadvsig, & dTsoil, tmp2d1, & dthetaldT, nsteps_ice, imelt, & deltaJ_latent_T, deltaJ_sensible_S, & qrunoff, tmp1d1, tmp1d2, tmp1d3, tmp1d4, & G0, & Tfreezing, LHS, RHS, LHS_h, RHS_h, nns, & iflux, litter, i, kk, condition, & advection, c2, theta, cp, & cpeff, hice, h0_tmp, hsnow, & delta_snowT, & delta_snowliq, J0, iqex, & nfac1, nfac2, nfac3, nfac4, nfac5, nfac6, nfac7, nfac8, nfac9, & nfac10, nfac11, nfac12, err & ) IMPLICIT NONE INTEGER(i_d) :: irec, mp REAL(r_2), DIMENSION(1:mp) :: qprec INTEGER(i_d) :: n REAL(r_2), DIMENSION(1:n) :: dx REAL(r_2), DIMENSION(1:mp) :: h0 REAL(r_2), DIMENSION(1:n) :: S REAL(r_2), DIMENSION(1:n) :: thetai REAL(r_2), DIMENSION(1:n) :: Tsoil REAL(r_2), DIMENSION(-nsnow_max:n) :: qh INTEGER(i_d), DIMENSION(1:mp) :: nsteps TYPE(vars), DIMENSION(1:mp) :: vlit TYPE(vars_snow), DIMENSION(1:mp) :: vsnow TYPE(vars), DIMENSION(1:n) :: var REAL(r_2), DIMENSION(1:mp) :: dxL TYPE(params), DIMENSION(1:mp) :: plit TYPE(params), DIMENSION(1:n) :: par REAL(r_2), DIMENSION(1:mp) :: deltaTa REAL(r_2), DIMENSION(1:mp) :: qevap REAL(r_2), DIMENSION(1:mp) :: qL, qhL, qybL, qTbL, qhTbL, qhybL LOGICAL, DIMENSION(1:mp) :: again LOGICAL, DIMENSION(1:mp,1:n) :: again_ice INTEGER(i_d), DIMENSION(1:mp) :: iok, itmp, ns REAL(r_2), DIMENSION(1:mp) :: accel, dt, fac, phip REAL(r_2), DIMENSION(1:mp) :: rsig, rsigdt, sig, t REAL(r_2), DIMENSION(1:n) :: qexd REAL(r_2), DIMENSION(-nsnow_max+1:n) :: aa, bb, cc, dd, ee, ff, gg, dy REAL(r_2), DIMENSION(-nsnow_max+1:n) :: aah, bbh, cch, ddh, eeh, ffh, ggh, de REAL(r_2), DIMENSION(-nsnow_max:n) :: q, qya, qyb, qTa, qTb,qhya, qhyb, qhTa, qhTb REAL(r_2), DIMENSION(-nsnow_max:n) :: qadv, qadvya, qadvyb, qadvTa, qadvTb REAL(r_2), DIMENSION(-nsnow_max:n) :: qsig, qhsig, qadvsig REAL(r_2), DIMENSION(1:n) :: dTsoil REAL(r_2), DIMENSION(0:n) :: tmp2d1 REAL(r_2), DIMENSION(1:n) :: dthetaldT INTEGER(i_d), DIMENSION(1:n) :: nsteps_ice, imelt REAL(r_2), DIMENSION(1:n) :: deltaJ_latent_T, deltaJ_sensible_S REAL(r_2), DIMENSION(1:mp) :: qrunoff REAL(r_2), DIMENSION(1:mp) :: tmp1d1, tmp1d2, tmp1d3, tmp1d4 REAL(r_2), DIMENSION(1:mp) :: G0 REAL(r_2), DIMENSION(1:mp) :: Tfreezing REAL(r_2), DIMENSION(-nsnow_max+1:n) :: LHS, RHS, LHS_h, RHS_h INTEGER(i_d), DIMENSION(1:mp) :: nns, iflux LOGICAL :: litter INTEGER(i_d) :: i, kk, condition INTEGER(i_d) :: advection ! switches REAL(r_2) :: c2, theta REAL(r_2), DIMENSION(1:mp) :: cp, cpeff, hice, h0_tmp REAL(r_2), DIMENSION(nsnow_max) :: hsnow REAL(r_2), DIMENSION(nsnow_max) :: delta_snowT, delta_snowliq REAL(r_2), DIMENSION(1:n) :: J0 REAL(r_2), DIMENSION(1:n) :: iqex INTEGER(i_d), DIMENSION(1:mp) :: nfac1, nfac2, nfac3, nfac4, nfac5, & nfac6, nfac7, nfac8, nfac9, nfac10, nfac11, nfac12 ! Error flag if nstep of SLI > nsteps_max: err=0 -> no error; err/=0 -> error INTEGER(i_d), INTENT(INOUT), OPTIONAL :: err iok(kk) = 0 ! flag for time step test itmp(kk) = 0 ! counter to abort if not getting solution again_ice(kk,:) = .FALSE. nsteps_ice(1:n) = 0 h0_tmp(kk) = h0(kk) DO WHILE (iok(kk)==0) ! keep reducing time step until all ok itmp(kk) = itmp(kk) + 1 accel(kk) = one - 0.05_r_2*REAL(MIN(10,MAX(0,itmp(kk)-4)),r_2) ! acceleration [0.5,1], start with 1 IF (itmp(kk) > 1000) THEN WRITE(logn,*) "Solve: too many iterations of equation solution" WRITE(logn,*) " irec, kk, S" WRITE(logn,*) irec, kk, S(:) WRITE(logn,*) " irec, kk, Tsoil" WRITE(logn,*) irec, kk, Tsoil(:) WRITE(logn,*) " irec, kk, qex" WRITE(logn,*) irec, kk, iqex(:) WRITE(logn,*) " irec, kk, h0, hsnow, hsnowliq" WRITE(logn,*) irec, kk, h0(kk), vsnow(kk)%hsnow, vsnow(kk)%hliq WRITE(logn,*) nfac1(kk), nfac2(kk), nfac3(kk), nfac4(kk), nfac5(kk), & nfac6(kk), nfac7(kk), nfac8(kk), nfac9(kk), nfac10(kk), nfac11(kk), nfac12(kk) err = 1 RETURN END IF ! prelim estimate of new top snow layer depth for use in energy cons eq'n ! IF ((vsnow(kk)%nsnow>0)) THEN hsnow(1:vsnow(kk)%nsnow) = vsnow(kk)%hsnow(1:vsnow(kk)%nsnow) + & (-q(1-vsnow(kk)%nsnow:0) + q(-vsnow(kk)%nsnow:-1))*dt(kk) ENDIF aa(1:n) = qya(0:n-1) ee(0:n-1) = -qyb(0:n-1) bb(1:n) = qTa(0:n-1) ff(0:n-1) = -qTb(0:n-1) gg(1:n) = -(q(0:n-1)-q(1:n)-iqex(1:n))*rsig(kk) gg(1) = gg(1)+qrunoff(kk)*rsig(kk) aah(1:n) = qhya(0:n-1) eeh(0:n-1) = -qhyb(0:n-1) bbh(1:n) = qhTa(0:n-1) ffh(0:n-1) = -qhTb(0:n-1) ggh(1:n) = -(qh(0:n-1)-qh(1:n))*rsig(kk) IF (advection==1) THEN ggh(1:n) = ggh(1:n) + iqex(1:n)*rsig(kk)*(Tsoil(1:n))*cswat*rhow ggh(1) = ggh(1) + qrunoff(kk)*rsig(kk)*(Tsoil(1))*cswat*rhow ENDIF IF (litter) THEN ! full litter model: litter in zeroth layer ! only use zeroth layer for litter (pond included in layer 1) WRITE(*,*) 'Should not be here - QYBL 01 ', qybl(kk) cc(0) = -qya(0) - rsigdt(kk)*plit(kk)%thre*dxL(kk) + qybL(kk) gg(0) = -(qprec(kk)-qevap(kk)-q(0))*rsig(kk) ggh(0) = -(G0(kk)-qh(0))*rsig(kk) dd(0) = -qTa(0) ENDIF aa(0) = qya(0) aah(0) = zero bbh(0) = zero WHERE (var(1:n)%isat==0) ! unsaturated layers cc(1:n) = qyb(0:n-1) - qya(1:n) - par(1:n)%thre*dx(1:n)*rsigdt(kk) - qexd(1:n) ELSEWHERE ! saturated layers cc(1:n) = qyb(0:n-1) - qya(1:n) - qexd(1:n) endwhere IF (ns(kk)<1) THEN ! pond included in top soil layer, solving for change in pond height cc(1) = -qya(1)-rsigdt(kk) -qexd(1) ENDIF cch(1:n) = qhyb(0:n-1)-qhya(1:n) + & REAL(var(1:n)%iice,r_2)*REAL(1-var(1:n)%isat,r_2)*rhow*lambdaf*par(1:n)%thre*dx(1:n)*rsigdt(kk) IF (ns(kk)==0) THEN ! change in pond height (top layer saturated) cch(1) = cch(1) + & REAL(var(1)%iice,r_2)*rhow*lambdaf*rsigdt(kk)*(var(1)%thetai/par(1)%thre) ENDIF ! modification to cch for advection IF (advection==1) THEN IF (ns(kk)==0) THEN ! changing pond included in top soil layer cch(1) = cch(1) -rhow*cswat*(Tsoil(1))*rsigdt(kk)*REAL(1-var(1)%iice,r_2) & -rhow*csice*(Tsoil(1))*rsigdt(kk)*REAL(var(1)%iice,r_2) & *(var(1)%thetai/par(1)%thre) & -rhow*cswat*(Tsoil(1))*rsigdt(kk)*REAL(var(1)%iice,r_2)* & (one-(var(1)%thetai/par(1)%thre)) ELSE ! no pond cch(1) = cch(1) -rhow*(Tsoil(1))*par(1)%thre*dx(1)*rsigdt(kk)* & REAL(1-var(1)%isat,r_2)*(cswat*REAL(1-var(1)%iice,r_2) & +csice*REAL(var(1)%iice,r_2)) ENDIF cch(2:n) = cch(2:n) -rhow*(Tsoil(2:n))*par(2:n)%thre*dx(2:n)*rsigdt(kk)* & REAL(1-var(2:n)%isat,r_2)*(cswat*REAL(1-var(2:n)%iice,r_2) & +csice*REAL(var(2:n)%iice,r_2)) ENDIF dd(1:n) = qTb(0:n-1)-qTa(1:n) ddh(1:n) = qhTb(0:n-1)-qhTa(1:n) - & ! Only apply latent heat component of heat capacity to total deltaT if soil remains frozen var(1:n)%csoileff*dx(1:n)*rsigdt(kk)- & (cswat-csice)*dx(1:n)*var(1:n)%dthetaldt*rhow*(Tsoil(1:n))* & rsigdt(kk)*REAL(var(1:n)%iice,r_2) ! modification to ddh(1) for pond ddh(1) = ddh(1) - cpeff(kk)*h0_tmp(kk)*rsigdt(kk) - & (cswat-csice)*h0(kk)/par(1)%thre*var(1)%dthetaldt*rhow*(Tsoil(1))* & rsigdt(kk)*REAL(var(1)%iice,r_2) IF (advection==1) THEN ddh(1:n) = ddh(1:n) - iqex(1:n)*cswat*rhow ddh(1) = ddh(1) - qrunoff(kk)*cswat*rhow ENDIF ! modification of matrix to incorporate single snow layer IF (vsnow(kk)%nsnow==1) THEN cc(0) = qyb(-1)-qya(0)-rsigdt(kk) dd(0) = qTb(-1)-qTa(0) ee(0) = -qyb(0) ff(0) = -qTb(0) gg(0) = -(q(-1)-q(0))*rsig(kk) IF (vsnow(kk)%hliq(1)>zero) THEN ! liquid phase present, solve for change in liq content ddh(0) = - rhow*((vsnow(kk)%tsn(1))*(cswat-csice)+lambdaf)*rsigdt(kk) ELSE ! solid phase only; solve for change in snow t ddh(0) = qhTb(-1) - qhTa(0) - rhow*csice*hsnow(1)*rsigdt(kk) ENDIF cch(0) = qhyb(-1)-qhya(0) - rhow*(csice*(vsnow(kk)%tsn(1))-lambdaf)*rsigdt(kk) eeh(0) = -qhyb(0) ffh(0) = -qhTb(0) ggh(0) = -(qh(-1)-qh(0))*rsig(kk) ENDIF ! modification of matrix to incorporate more than one snow layer IF (vsnow(kk)%nsnow>1) THEN aa(1-vsnow(kk)%nsnow:0) = qya(-vsnow(kk)%nsnow:-1) bb(1-vsnow(kk)%nsnow:0) = qTa(-vsnow(kk)%nsnow:-1) cc(1-vsnow(kk)%nsnow:0) = qyb(-vsnow(kk)%nsnow:-1)-qya(1-vsnow(kk)%nsnow:0)-rsigdt(kk) dd(1-vsnow(kk)%nsnow:0) = qTb(-vsnow(kk)%nsnow:-1)-qTa(1-vsnow(kk)%nsnow:0) ee(1-vsnow(kk)%nsnow:0) = -qyb(1-vsnow(kk)%nsnow:0) ff(1-vsnow(kk)%nsnow:0) = -qTb(1-vsnow(kk)%nsnow:0) gg(1-vsnow(kk)%nsnow:0) = -(q(-vsnow(kk)%nsnow:-1)-q(1-vsnow(kk)%nsnow:0))*rsig(kk) ddh(1-vsnow(kk)%nsnow:0) = MERGE( & ! liquid phase present, solve for change in liq content -rhow*((vsnow(kk)%tsn(1:vsnow(kk)%nsnow))*(cswat-csice)+lambdaf)*rsigdt(kk), & qhTb(-vsnow(kk)%nsnow:-1) - qhTa(1-vsnow(kk)%nsnow:0) - & rhow*csice*hsnow(1:vsnow(kk)%nsnow)*rsigdt(kk), & ! solid phase only; solve for change in snow t vsnow(kk)%hliq(1:vsnow(kk)%nsnow) > zero) aah(1-vsnow(kk)%nsnow:0) = qhya(-vsnow(kk)%nsnow:-1) bbh(1-vsnow(kk)%nsnow:0) = qhTa(-vsnow(kk)%nsnow:-1) cch(1-vsnow(kk)%nsnow:0) = qhyb(-vsnow(kk)%nsnow:-1)-qhya(1-vsnow(kk)%nsnow:0) - & rhow*(csice*(vsnow(kk)%tsn(1:vsnow(kk)%nsnow))-lambdaf)*rsigdt(kk) eeh(1-vsnow(kk)%nsnow:0) = -qhyb(1-vsnow(kk)%nsnow:0) ffh(1-vsnow(kk)%nsnow:0) = -qhTb(1-vsnow(kk)%nsnow:0) ggh(1-vsnow(kk)%nsnow:0) = -(qh(-vsnow(kk)%nsnow:-1)-qh(1-vsnow(kk)%nsnow:0))*rsig(kk) ENDIF IF (litter .AND. ns(kk)==1) THEN ! litter and no pond ! watch for deltaTa WRITE(*,*) 'Should not be here - QYBL 02 ', qybl(kk) cc(0) = qybL(kk) -qya(0) -rsigdt(kk)*plit(kk)%thre*dxL(kk) dd(0) = qTbL(kk) - qTa(0) ee(0) = -qyb(0) ff(0) = -qTb(0) gg(0) = (q(0) - qL(kk))*rsig(kk) + deltaTa(kk)*(qTbL(kk)-qTa(0)) gg(1) = -(q(0)-q(1)-iqex(1))*rsig(kk) + deltaTa(kk)*qTa(0) cch(0) = qhybL(kk) - qhya(0) ddh(0) = -qhTa(0)-vlit(kk)%csoil*dxL(kk)*rsigdt(kk)+qhTbL(kk) eeh(0) = -qhyb(0) ffh(0) = -qhTb(0) ggh(0) = (qh(0)-qhL(kk))*rsig(kk) + deltaTa(kk)*(qhTbL(kk) - qhTa(0)) ggh(1) = -(qh(0)-qh(1))*rsig(kk) + deltaTa(kk)*qhTa(0) ENDIF ! litter and pond !vh! need to check this now that pond is lumped with top soil layer IF (litter .AND. ns(kk)==0) THEN ddh(0) = -qhTa(0)-cswat*rhow*h0(kk)*rsigdt(kk) -vlit(kk)%csoil*dxL(kk)*rsigdt(kk) ENDIF nns(kk) = 1 ! pond included in top soil layer IF (vsnow(kk)%nsnow>0) THEN nns(kk) = 1-vsnow(kk)%nsnow ENDIF CALL massman_sparse(aa(nns(kk)+1:n), aah(nns(kk)+1:n), bb(nns(kk)+1:n), bbh(nns(kk)+1:n), & cc(nns(kk):n), cch(nns(kk):n), dd(nns(kk):n), ddh(nns(kk):n), ee(nns(kk):n-1), & eeh(nns(kk):n-1), & ff(nns(kk):n-1), ffh(nns(kk):n-1), gg(nns(kk):n), ggh(nns(kk):n), & dy(nns(kk):n), de(nns(kk):n), condition=condition, err=err) IF (err /= 0) THEN WRITE(logn,*) "Sparse matrix solution failed ", irec, kk WRITE(logn,*) Tsoil(1), S(1) RETURN ENDIF dTsoil(1:n) = de(1:n) WHERE (vsnow(kk)%hliq(:)>zero) delta_snowT(:) = zero ELSEWHERE delta_snowT(:) = de(LBOUND(de,1):0) endwhere ! evaluate soil fluxes at sigma of time step qsig(0) = q(0)+sig(kk)*qyb(0)*dy(1) + sig(kk)*qya(0)*dy(0) + & sig(kk)*qTb(0)*dTsoil(1) + sig(kk)*qTa(0)*de(0) qhsig(0) = qh(0) + sig(kk)*qhyb(0)*dy(1) + sig(kk)*qhya(0)*dy(0) + & sig(kk)*qhTb(0)*dTsoil(1) + sig(kk)*qhTa(0)*de(0) qadvsig(0) = qadv(0) + sig(kk)*qadvyb(0)*dy(1) & + sig(kk)*qadvya(0)*dy(0) + sig(kk)*qadvTb(0)*dTsoil(1) + sig(kk)*qadvTa(0)*de(0) qsig(1:n-1) = q(1:n-1) + sig(kk)*(qya(1:n-1)*dy(1:n-1)+qyb(1:n-1)*dy(2:n) & +qTa(1:n-1)*dTsoil(1:n-1)+qTb(1:n-1)*dTsoil(2:n)) qsig(n) = q(n) + sig(kk)*(qya(n)*dy(n)+qTa(n)*dTsoil(n)) qhsig(1:n-1) = qh(1:n-1) + sig(kk)*(qhya(1:n-1)*dy(1:n-1)+qhyb(1:n-1)*dy(2:n) & +qhTa(1:n-1)*dTsoil(1:n-1)+qhTb(1:n-1)*dTsoil(2:n)) qhsig(n) = qh(n) + sig(kk)*(qhya(n)*dy(n)+qhTa(n)*dTsoil(n)) qadvsig(1:n-1) = qadv(1:n-1) + sig(kk)*(qadvya(1:n-1)*dy(1:n-1)+qadvyb(1:n-1)*dy(2:n) & +qadvTa(1:n-1)*dTsoil(1:n-1)+qadvTb(1:n-1)*dTsoil(2:n)) qadvsig(n) = qadv(n) + sig(kk)*(qadvya(n)*dy(n)+qadvTa(n)*dTsoil(n)) LHS_h(1) = (dx(1)*var(1)%csoileff+h0_tmp(kk)*cpeff(kk))*dTsoil(1)/dt(kk) - & dy(1)*par(1)%thre*dx(1)*rhow*lambdaf/dt(kk)*var(1)%iice* & REAL(1-var(1)%isat,r_2) - & dy(1)*(var(1)%thetai/par(1)%thre)*rhow*lambdaf/dt(kk)*REAL(1-ns(kk),r_2)* & var(1)%iice + & var(1)%dthetaldT*(Tsoil(1))*rhow*(cswat-csice)*dx(1)*dTsoil(1)/dt(kk)* & var(1)%iice + & var(1)%dthetaldT*(Tsoil(1))*rhow*(cswat-csice)*h0(kk)/par(1)%thre* & dTsoil(1)/dt(kk)*var(1)%iice LHS_h(2:n) = (dx(2:n)*var(2:n)%csoileff)*dTsoil(2:n)/dt(kk) - & dy(2:n)*par(2:n)%thre*dx(2:n)*rhow*lambdaf/dt(kk)*var(2:n)%iice* & REAL(1-var(2:n)%isat,r_2) + & var(2:n)%dthetaldT*(Tsoil(2:n))*rhow*(cswat-csice)*dx(2:n)* & dTsoil(2:n)/dt(kk)*var(2:n)%iice RHS_h(1:n) = qhsig(0:n-1) -qhsig(1:n) IF (advection==1) THEN LHS_h(1:n) = LHS_h(1:n) & + REAL(1-var(1:n)%iice,r_2) * REAL(1-var(1:n)%isat,r_2) * & dx(1:n) * dy(1:n)/dt(kk) * par(1:n)%thre * (Tsoil(1:n)) * rhow * cswat & + REAL(var(1:n)%iice,r_2) * REAL(1-var(1:n)%isat,r_2) * & dx(1:n) * dy(1:n)/dt(kk) * par(1:n)%thre * rhow * csice * (Tsoil(1:n)) LHS_h(1) = LHS_h(1) & + REAL(1-ns(kk),r_2) * REAL(1-var(1)%iice,r_2) * & dy(1)/dt(kk) * rhow * cswat * (Tsoil(1)) & + REAL(1-ns(kk),r_2) * REAL(var(1)%iice,r_2)* & dy(1)/dt(kk) * rhow * csice * (Tsoil(1)) * (var(1)%thetai/par(1)%thre) & + REAL(1-ns(kk),r_2) * REAL(var(1)%iice,r_2) * & dy(1)/dt(kk) * rhow * cswat * (Tsoil(1)) * (one-(var(1)%thetai/par(1)%thre)) RHS_h(1:n) = RHS_h(1:n) - iqex(1:n)*cswat*rhow*(