MODULE sli_utils USE cable_def_types_mod, ONLY: r_2, i_d USE cable_def_types_mod, ONLY: soil_parameter_type, veg_parameter_type USE sli_numbers, ONLY: & experiment, & zero, half, one, two, four, e3, pi, & Tzero, gravity, Rgas, thousand, Mw, rlambda, Dva, & params, vars_aquifer, vars, rapointer, & rhow, lambdaf, lambdas, csice, cswat, cpa,& dpmaxr, solve_type, & gf, hmin, csol, rhmin, dsmmax, rhocp, vars_snow, vars_met, & freezefac, ithermalcond, rmair, Mw IMPLICIT NONE PRIVATE PUBLIC :: dx, dxL, par, plit, sol, x ! soil water parameters PUBLIC :: bd, dis, isopar, isotype ! soil solute parameters PUBLIC :: aquifer_props, flux, generic_thomas, getfluxes_vp, getheatfluxes, hyofh, hyofS, isosub ! subroutines PUBLIC :: litter_props, massman_sparse, potential_evap, setlitterpar, setpar, setpar_Loetsch, setsol, setx, tri PUBLIC :: csat, csoil, dthetalmaxdT, dthetalmaxdTh, esat, esat_ice, gammln, igamma, phi, rh0_sol, rtbis_rh0 ! functions PUBLIC :: slope_csat, slope_esat,slope_esat_ice, Sofh, Tfrz, thetalmax, weight, zerovars, Tthetalmax, Tfrozen PUBLIC :: rtbis_Tfrozen, GTfrozen, JSoilLayer, forcerestore, SEB PUBLIC :: spline_b, mean, nse REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: dx REAL(r_2), DIMENSION(:), ALLOCATABLE :: dxL TYPE(params), DIMENSION(:,:), ALLOCATABLE :: par TYPE(params), DIMENSION(:), ALLOCATABLE :: plit TYPE(solve_type), DIMENSION(:), ALLOCATABLE :: sol REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: x !MC solute not done yet !REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: bd !REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: dis !TYPE(rapointer), DIMENSION(:,:,:), ALLOCATABLE :: isopar !CHARACTER(LEN=2), DIMENSION(:,:,:), ALLOCATABLE :: isotype REAL(r_2), DIMENSION(:), ALLOCATABLE :: bd ! function solute not done yet REAL(r_2), DIMENSION(:), ALLOCATABLE :: dis TYPE(rapointer), DIMENSION(:,:), ALLOCATABLE :: isopar CHARACTER(LEN=2), DIMENSION(:,:), ALLOCATABLE :: isotype ! Subroutine interfaces INTERFACE generic_thomas MODULE PROCEDURE generic_thomas_1d MODULE PROCEDURE generic_thomas_2d END INTERFACE INTERFACE getfluxes_vp MODULE PROCEDURE getfluxes_vp_1d MODULE PROCEDURE getfluxes_vp_2d END INTERFACE INTERFACE getheatfluxes MODULE PROCEDURE getheatfluxes_1d MODULE PROCEDURE getheatfluxes_2d END INTERFACE INTERFACE massman_sparse MODULE PROCEDURE massman_sparse_1d MODULE PROCEDURE massman_sparse_2d END INTERFACE INTERFACE tri MODULE PROCEDURE tri_1d MODULE PROCEDURE tri_2d END INTERFACE ! Function interfaces INTERFACE mean MODULE PROCEDURE mean_1d, mean_2d END INTERFACE INTERFACE nse MODULE PROCEDURE nse_1d, nse_2d END INTERFACE !P.J. Ross 2005-2007: ! This module implements Brooks-Corey (BC) soil water retention and conductivity ! functions. It illustrates the structure required without the complicating ! detail of more flexible hydraulic property functions. ! Definitions of public entities (see above for default values): ! params - type for water parameters. Params the, thre (=the-thr), he, lam, ! Ke and eta are for the BC functions. Params he, Ke, KSe, phie ! and phiSe are values of variables h, K, KS, phi and phiS at ! saturation (denoted by the "e" for "air entry"), needed by ! module flow (MF). ! vars - type for water variables used by MF and returned by subroutine ! hyofS (except for isat, which is 0 for unsaturated layers and 1 ! for saturated layers). ! gf - gravity factor for flow direction (usually 1 for straight down). ! hmin - minimum matric head h (used by MF). ! par(:) - hydraulic property params for soil types (used by MF). ! allo - subroutine to allocate parameter storage. ! hyofS - subroutine to get water variable from saturation S (where S<1). ! hyofh - subroutine to get some water variables from h. ! Sofh - subroutine to get S from h. ! hypar - subroutine to set soil hydraulic params. ! weight - subroutine to get gravity flow conductivity weight w. !V. Haverd 2008: ! modified by vh to include isothermal vapour conductivity in phi ! change units from cm, h to m, s ! includes thermal properties in soil parameters ! litter_props - subroutine to define litter properties !============================================================================= CONTAINS !============================================================================= ! SUBROUTINES - alphabetical !============================================================================= ELEMENTAL PURE SUBROUTINE aquifer_props(v_aquifer) IMPLICIT NONE TYPE(vars_aquifer), INTENT(INOUT) :: v_aquifer v_aquifer%zzero = 53.43_r_2 ! water table depth corresponding to Wa = zero v_aquifer%Sy = 0.2_r_2 ! specific yield of aquifer ! initialise water content of aquifer v_aquifer%Wa = v_aquifer%Sy*(v_aquifer%zzero-MAX(v_aquifer%zdelta,v_aquifer%zsoil)) v_aquifer%isat = 0 IF (v_aquifer%zdelta <= v_aquifer%zsoil) v_aquifer%isat = 1 v_aquifer%f = 1.25_r_2 ! multiplier in exponent of Rs (m-1) v_aquifer%Rsmax = 4.5e-7_r_2 ! maximum discharge rate from aquifer (ms-1) v_aquifer%discharge = v_aquifer%Rsmax*EXP(-v_aquifer%f*v_aquifer%zdelta) END SUBROUTINE aquifer_props !============================================================================= ELEMENTAL PURE SUBROUTINE flux(parin, v1, v2, dz, q, qya, qyb, qTa, qTb) ! VH modified 25/05/10 to include t-dep component of liquid flux (in frozen soil) IMPLICIT NONE TYPE(params), INTENT(IN) :: parin REAL(r_2), INTENT(IN) :: dz REAL(r_2), INTENT(OUT) :: q, qya, qyb, qTa, qTb TYPE(vars), INTENT(IN) :: v1, v2 ! Gets flux and partial derivs for specified flow path. ! Definitions of arguments: ! j - soil type no. ! v1 - water vars at upper end of path. ! v2 - ditto at lower end. ! dz - length of path. ! q - flux. ! qya - partial deriv of flux wrt S (if unsat) or phi (if sat) at upper end. ! qyb - ditto at lower end. REAL(r_2) :: w, rdz ! gf is gravity factor (0 to 1) assumed available in module IF (gf < zero) THEN !if ((v1%isat /= 0 .and. v2%isat /= 0) .or. v1%h-gf*(-dz) >= parin%he) then IF ((v1%isat /= 0 .AND. v2%isat /= 0) .OR. v1%h-gf*(-dz) >= v1%he) THEN w = one ELSE w = weight(parin, v1%h, v1%K*v1%macropore_factor, v1%phi, -dz) w = one-w END IF ELSE !if ((v1%isat /= 0 .and. v2%isat /= 0) .or. v2%h-gf*dz >= parin%he) then IF ((v1%isat /= 0 .AND. v2%isat /= 0) .OR. v2%h-gf*dz >= v2%he) THEN w = zero ELSE w = weight(parin, v2%h, v2%K*v2%macropore_factor, v2%phi, dz) END IF END IF rdz = one/dz q = (v1%phi-v2%phi)*rdz + gf*(w*v1%K*v1%macropore_factor+(one-w)*v2%K*v2%macropore_factor) IF (v1%isat==0) THEN qya = v1%phiS*rdz + gf*w*v1%KS*v1%macropore_factor qTa = v1%phiT*rdz + gf*w*v1%KT*v1%macropore_factor ELSE qya = rdz qTa = zero END IF IF (v2%isat==0) THEN qyb = -v2%phiS*rdz + gf*(one-w)*v2%KS*v2%macropore_factor qTb = -v2%phiT*rdz + gf*(one-w)*v2%KT*v2%macropore_factor ELSE qyb = -rdz qTb = zero END IF END SUBROUTINE flux !============================================================================= SUBROUTINE forcerestore(Tg0, Rnet0, lE0, dlEdTg, Ta, Tbar, d1, rrc, lambda, & cs, dt, iice, Tg, G, H, lE) ! method applicable to multilayer soil or soil/snow column ! derived using Eq's 3-12 in Hirota et al. JGR 2002 IMPLICIT NONE REAL(r_2), INTENT(IN) :: Tg0 ! ground surface temp of previous time-step [deg C] REAL(r_2), INTENT(IN) :: Rnet0 ! Rnet at current time step, assuming Tg of previous time-step [W m-2] REAL(r_2), INTENT(IN) :: lE0 ! latent heat flux at current time step, assuming Tg of previous time-step [W m-2] REAL(r_2), INTENT(IN) :: dlEdTg ! derivative of latent heat flux wrt Tg [W m-2 K-1] REAL(r_2), INTENT(IN) :: Ta ! air temperature [deg C] REAL(r_2), INTENT(IN) :: Tbar ! temperature at diurnal damping depth REAL(r_2), INTENT(IN) :: d1 ! diurnal damping depth (m) = sqrt(2*lambda/c/omega) (Hirota et al. eq 40) REAL(r_2), INTENT(IN) :: rrc ! resistance to sensible heat and radiation transfer at ground/air interface [m-1 s] REAL(r_2), INTENT(IN) :: lambda ! thermal conductivity of soil or snow at surface [W m-1 K-1] REAL(r_2), INTENT(IN) :: cs ! heat capacity of soil or snow at surface [J m-3 K-1] REAL(r_2), INTENT(IN) :: dt ! time step [s] INTEGER, INTENT(IN):: iice ! top layer frozen (1) or not (0) REAL(r_2), INTENT(OUT) :: Tg, G, H, lE ! local variables REAL(r_2), PARAMETER :: tau1 = 86400._r_2 ! period of diurnal forcing (seconds) REAL(r_2) :: c1, c2, omega, a, b a = Rnet0 - rhocp/rrc*(Tg0-Ta) -lE0 ! G0 at current time step, assuming surface T of previous time step omega = 2._r_2*pi/tau1 ! diurnal forcing frequency (s-1) c1 = omega * d1 / lambda c2 = omega b = -rhocp/rrc - dlEdTg !Tg = (Tg0 + dt*(c1*a+c2*Tbar)) / (one + c2 *dt) Tg = (Tg0 + dt*(c1*(a-b*Tg0)+c2*Tbar)) / (one + c2 *dt - b*c1*dt) IF (iice.EQ.1) THEN Tg = MIN(zero, Tg) ENDIF G = a + b*(Tg-Tg0) H = rhocp/rrc*(Tg-Ta) lE = lE0 + dlEdTg*(Tg-Tg0) END SUBROUTINE forcerestore !============================================================================= SUBROUTINE forcerestore_Deardorff(Tg0, Rnet0, lE0, dlEdTg, Ta, Tbar, d1, rrc, rhos, & cs, dt, iice, Tg, G, H, lE) IMPLICIT NONE REAL(r_2), INTENT(IN) :: Tg0 ! ground surface temp of previous time-step [deg C] REAL(r_2), INTENT(IN) :: Rnet0 ! Rnet at current time step, assuming Tg of previous time-step [W m-2] REAL(r_2), INTENT(IN) :: lE0 ! latent heat flux at current time step, assuming Tg of previous time-step [W m-2] REAL(r_2), INTENT(IN) :: dlEdTg ! derivative of latent heat flux wrt Tg [W m-2 K-1] REAL(r_2), INTENT(IN) :: Ta ! air temperature [deg C] REAL(r_2), INTENT(IN) :: Tbar ! temperature at diurnal damping depth REAL(r_2), INTENT(IN) :: d1 ! diurnal damping depth (m) REAL(r_2), INTENT(IN) :: rrc ! resistance to sensible heat and radiation transfer at ground/air interface [m-1 s] REAL(r_2), INTENT(IN) :: rhos ! density of soil or snow at surface [kg m-3] REAL(r_2), INTENT(IN) :: cs ! heat capacity of soil or snow at surface [J kg-1 K-1] REAL(r_2), INTENT(IN) :: dt ! time step [s] INTEGER, INTENT(IN):: iice ! top layer frozen (1) or not (0) REAL(r_2), INTENT(OUT) :: Tg, G, H, lE ! local variables REAL(r_2), PARAMETER :: c1 = 3.72_r_2 ! Deardorff JGR (1978) REAL(r_2), PARAMETER :: c2 = 7.4_r_2 ! Deardorff JGR (1978) REAL(r_2), PARAMETER :: tau1 = 86400._r_2 ! period of diurnal forcing (seconds) REAL(r_2) :: a, b a = Rnet0 - rhocp/rrc*(Tg0-Ta) -lE0 b = -rhocp/rrc - dlEdTg Tg = Tg0 + (c1*a/(rhos*cs*d1) - c2/tau1*(Tg0-Tbar))/ & (1._r_2/dt - b*c1/(rhos*cs*d1) + c2/tau1) IF (iice.EQ.1) THEN Tg = MIN(zero, Tg) ENDIF G = a + b*(Tg-Tg0) H = rhocp/rrc*(Tg-Ta) lE = lE0 + dlEdTg*(Tg-Tg0) END SUBROUTINE forcerestore_Deardorff !============================================================================= ! Surface Energy Balance SUBROUTINE SEB(n, par, vmet, vsnow, var, qprec, qprec_snow, dx, h0, Tsoil, & Tsurface, G0, lE0, Epot, qsurface, qevap, qliq, qv, & qyb, qTb, qlyb, qvyb, qlTb, qvTb, qh, qadv, qhyb, qhTb, qadvyb, qadvTb, irec) IMPLICIT NONE INTEGER(i_d), INTENT(IN) :: n TYPE(params), DIMENSION(1:n), INTENT(IN) :: par TYPE(vars_met), INTENT(IN) :: vmet TYPE(vars_snow), INTENT(IN) :: vsnow TYPE(vars), DIMENSION(1:n), INTENT(IN) :: var REAL(r_2), INTENT(IN) :: qprec REAL(r_2), INTENT(IN) :: qprec_snow INTEGER(i_d), INTENT(IN) :: irec REAL(r_2), DIMENSION(1:n), INTENT(IN) :: dx REAL(r_2), INTENT(IN) :: h0 REAL(r_2), DIMENSION(1:n), INTENT(IN) :: Tsoil REAL(r_2), INTENT(OUT) :: Tsurface, G0, lE0, Epot ! SEB REAL(r_2), INTENT(OUT) :: qsurface ! water flux into surface REAL(r_2), INTENT(OUT) :: qevap ! evaporative water flux ! liquid and vapour components of water flux from surface into soil REAL(r_2), INTENT(OUT) :: qliq, qv ! derivatives of water fluxes wrt moisture and T REAL(r_2), INTENT(OUT) :: qyb, qTb, qlyb, qvyb, qlTb, qvTb ! total and advective components of heat flux into surface REAL(r_2), INTENT(OUT) :: qh, qadv ! derivatives of heat fluxes wrt moiture and T REAL(r_2), INTENT(OUT) :: qhyb, qhTb, qadvyb, qadvTb ! local variables INTEGER(i_d) :: surface_case REAL(r_2) :: Tsurface_pot, Hpot, Gpot, dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil REAL(r_2) :: E_vap, dE_vapdT1, E_liq REAL(r_2) :: Kmin, Khmin, phimin REAL(r_2) :: Tqw, dtqwdtb, rhocp1, cs LOGICAL :: isEpot IF (vsnow%nsnow.EQ.0) surface_case = 1 IF (vsnow%nsnow>0) surface_case = 2 SELECT CASE (surface_case) CASE (1) ! no snow CALL potential_evap(vmet%Rn, vmet%rbh, vmet%rbw, vmet%Ta, vmet%rha, & Tsoil(1), var(1)%kth, half*dx(1)+h0, var(1)%lambdav, Tsurface_pot, Epot, Hpot, & Gpot, dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil) IF (var(1)%iice.EQ.1.AND.Tsurface_pot> zero) THEN Tsurface_pot = 0.0_r_2 Tsurface = 0.0_r_2 Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) vmet%cva)*rhow*var(1)%lambdav/vmet%rbw dEdTsoil = zero dGdTsoil = zero Hpot = rhocp*(Tsurface - vmet%Ta)/vmet%rbh Gpot = vmet%Rn - Hpot - Epot dEdTs = zero ENDIF IF (var(1)%isat.EQ.1) THEN ! saturated surface =>. potential evporation isEpot = .TRUE. Tsurface = Tsurface_pot lE0 = Epot G0 = Gpot E_vap = zero dE_vapdT1 = zero E_liq = lE0 ELSE ! unsaturated surface: finite vapour transfer; surface flux may be supply limited CALL hyofh(hmin, par(1)%lam, par(1)%eta, par(1)%Ke, par(1)%he, & Kmin, Khmin, phimin) ! get phi at hmin E_liq = ((var(1)%phi-phimin)/(half*dx(1))-var(1)%K)*thousand*var(1)%lambdav IF (var(1)%Dv > 1.e-12_r_2) THEN ! E = (cs-ca)/rbw = E_liq + E_vap = E_liq + (c1-cs)/Dv/dx/2 cs = ( E_liq/var(1)%lambdav + var(1)%rh*csat(Tsoil(1))*var(1)%Dv/(half*dx(1)) + & vmet%cva*thousand/vmet%rbw ) & / (one/vmet%rbw + var(1)%Dv/(half*dx(1))) E_vap = (var(1)%rh*csat(Tsoil(1)) - cs) * var(1)%Dv/(half*dx(1)) * var(1)%lambdav ! dE_vapdT1 = var(1)%Dv/(half*dx(1))*var(1)%rh*slope_csat(Tsoil(1)) & ! * (one - vmet%rbw*var(1)%Dv/(half*dx(1))/(one+vmet%rbw*var(1)%Dv/(half*dx(1)))) & ! * var(1)%lambdav dE_vapdT1 = var(1)%Dv/(half*dx(1)) * var(1)%rh*slope_csat(Tsoil(1)) * var(1)%lambdav ELSE E_vap = zero dE_vapdT1 = zero ENDIF IF (Epot <= (E_vap+E_liq)) THEN isEpot = .TRUE. lE0 = Epot cs = csat(Tsurface_pot) E_vap = (var(1)%rh*csat(Tsoil(1)) - cs) * var(1)%Dv/(half*dx(1)) * var(1)%lambdav E_liq = lE0 - E_vap dE_vapdT1 = var(1)%Dv/(half*dx(1)) * var(1)%rh*slope_csat(Tsoil(1)) * var(1)%lambdav ELSE isEpot = .FALSE. lE0 = E_vap+E_liq dEdTs = zero ENDIF ! lE0 = min(Epot, E_vap+E_liq) ! analytic approximation (See Haverd et al. 2013, Appxx) ! if (Epot .gt. (E_vap+E_liq)) dEdTs = zero Tsurface = (-half*dx(1)*lE0 + half*dx(1)*vmet%Rn + & var(1)%kth*Tsoil(1) + half*dx(1)*(one/vmet%rbh*rhocp)*vmet%Ta) & /(var(1)%kth + half*dx(1)*(one/vmet%rbh*rhocp)) G0 = var(1)%kth/(half*dx(1))*(Tsurface-Tsoil(1)) dGdTsoil = -var(1)%kth/(half*dx(1)) IF ((var(1)%iice.EQ.1) .AND. (Tsurface>zero)) THEN Tsurface = 0.0_r_2 rhocp1 = rmair*101325._r_2/rgas/(vmet%Ta+Tzero)*cpa G0 = vmet%Rn - rhocp1*(Tsurface - vmet%Ta)/vmet%rbh - lE0 dGdTsoil = 0.0_r_2 ENDIF ENDIF ! write(*,*) var(1)%phi, phimin qevap = lE0/(thousand*var(1)%lambdav) qsurface = qprec + qprec_snow - qevap ! derivatives ! q refers to moisture in numerator ! qh refers to heat in numerator ! y refers to moisture in denominator ! T refers to temperature in denominator ! a refers to the layer above ! b refers to the layer below ! initialise derivatives to zero qyb = zero qTb = zero qhyb = zero qhTb = zero ! liquid and vapour fluxes qlyb = zero qvyb = zero ! potential evap independent of S(1), dependent on T1 IF ((var(1)%isat.EQ.1 .OR. vsnow%nsnow.GT.0)) THEN qyb = zero qTb = -dEdTsoil/(thousand*var(1)%lambdav) qliq = -Epot/(thousand*var(1)%lambdav) qv = zero qlyb = zero qvyb = zero qlTb = qTb qvTb = zero ELSEIF (isEpot) THEN qTb = -dE_vapdT1/(thousand*var(1)%lambdav) qyb = zero qliq = -E_liq/(thousand*var(1)%lambdav) qv = -E_vap/(thousand*var(1)%lambdav) qlyb = zero qvyb = zero qlTb = qTb qvTb = zero ELSE ! supply limited qTb = -dE_vapdT1/(thousand*var(1)%lambdav) qyb = -(var(1)%phiS/(half*dx(1)) - var(1)%KS) !vh! include vapour component?? qliq = -E_liq/(thousand*var(1)%lambdav) qv = -E_vap/(thousand*var(1)%lambdav) qlyb = -(var(1)%phiS/(half*dx(1)) - var(1)%KS) qvyb = zero qlTb = zero qvTb = qTb ENDIF ! end of partial derivative evaluation ! advective component of heat flux qadv = rhow*cswat*qprec*(vmet%Ta) + rhow*csice*qprec_snow*(MIN(vmet%Ta,zero)) & - rhow*qprec_snow*lambdaf Tqw = MERGE(vmet%Ta, Tsoil(1), (-qevap)>zero) dTqwdTb = MERGE(zero, one, (-qevap)>zero) qadv = qadv + rhow*cswat*Tqw*(-qevap) qadvTb = dTqwdTb + rhow*cswat*Tqw*qTb qadvyb = rhow*cswat*qyb*Tqw qadvyb = 0 qadvTb = 0 ! test vh! qh = qadv + G0 qhyb = qadvyb qhTb = dGdTsoil + qadvTb CASE (2) !dedicated snow layer ! NB Only longwave component of net radiation directly affects SEB: sw component is absorbed internally ! SEB at snow/air interface IF (vsnow%hliq(1)>zero) THEN Tsurface = 0.0_r_2 !Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) ! vmet%cva)*rhow*rlambda/vmet%rbw !vh check this ! Epot = (csat(Tsurface)/thousand - vmet%cva)/vmet%rbw *rlambda*rhow ! m3 H2O (liq) m-3 (air) -> W/m2 ! write(*,*) "Epot", vmet%rha, vmet%Ta, ! esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero)*rhow*rlambda/vmet%rbw, & !vmet%cva*rhow*rlambda/vmet%rbw, Epot dEdTsoil = zero dGdTsoil = zero rhocp1 = rmair*101325._r_2/rgas/(vmet%Ta+Tzero)*cpa Hpot = rhocp1*(Tsurface - vmet%Ta)/vmet%rbh Gpot = vmet%Rn-vmet%Rnsw - Hpot - Epot dEdTs = zero qevap = Epot/(rhow*rlambda) qTb = zero ! write(*,*) "Epot2", Tsurface, vmet%Ta, Epot, Hpot, vmet%rrc, rhocp1 ELSE ! vh ! use max snow depth of 20 cm in this calculation to avoid huge resistances ! leading to large negative surface temperatures when snow-pack is thick and ! Rn is large and negative (~-100 Wm-2) CALL potential_evap(vmet%Rn-vmet%Rnsw, vmet%rbh, vmet%rbw, vmet%Ta, vmet%rha, & vsnow%tsn(1), MAX(vsnow%kth(1),0.1_r_2), half*MIN(vsnow%depth(1),0.05_r_2), & lambdas, Tsurface, Epot, Hpot, & Gpot, dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil,iice=.TRUE.) IF (Tsurface > zero) THEN ! temperature of frozen surface must be <= zero Tsurface = 0.0_r_2 Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) vmet%cva)*rhow*lambdas/vmet%rbw dEdTsoil = zero dGdTsoil = zero Hpot = rhocp*(Tsurface - vmet%Ta)/vmet%rbh Gpot = vmet%Rn-vmet%Rnsw - Hpot - Epot dEdTs= zero ! write(*,*) "Epot3", Tsurface, vmet%Ta, Epot, Hpot, vmet%rbh ENDIF !$ elseif (abs(Tsurface - vmet%Ta).gt. 20) then !$ Tsurface = min(vmet%Ta, 0.0) !$ Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) !$ vmet%cva)*rhow*lambdas/vmet%rbw !$ dEdTsoil = zero !$ dGdTsoil = zero !$ Hpot = rhocp*(Tsurface - vmet%Ta)/vmet%rbh !$ Gpot = vmet%Rn-vmet%Rnsw - Hpot - Epot !$ dEdTs= zero !$ endif qevap = Epot/(rhow*lambdas) qTb = -dEdTsoil/(thousand*lambdas) ENDIF lE0 = Epot ! moisture flux at air/snow interface qsurface = qprec_snow+qprec-qevap qyb = zero ! conductive heat flux at air/snow interface G0 = Gpot qhTb = dGdTsoil IF (vsnow%hliq(1)>zero) THEN qhTb = zero qTb = zero ENDIF ! advective heat flux at air/snow interface qadv = rhow*(qprec_snow)*(csice*(MIN(vmet%Ta,zero))-lambdaf) + & rhow*(qprec)*cswat*(MAX(vmet%Ta,zero)) Tqw = MERGE(vmet%Ta, vsnow%tsn(1), -qevap>zero) dTqwdTb = MERGE(zero,one, -qevap>zero) IF (vsnow%hliq(vsnow%nsnow)>zero) THEN qadv = qadv + rhow*(-qevap)*cswat*Tqw qadvTb = zero ELSE qadv = qadv + rhow*(-qevap)*cswat*Tqw qadvTb = rhow*cswat*(-qevap)*dTqwdTb + rhow*cswat*Tqw*qTb ENDIF ! add sw energy absorption to G0 G0 = G0 + vmet%Rnsw qh = qadv + G0 qhyb = zero qadvyb = zero qhyb = qhyb + qadvyb qhTb = qhTb + qadvTb qv = -qevap qliq = zero qvyb = qyb qvTb = qTb qlyb = zero qlTb = zero END SELECT ! surface_case ! finished all the surfaces END SUBROUTINE SEB !============================================================================= ! Surface Energy Balance SUBROUTINE SEB_FR(n, par, vmet, vsnow, var, qprec, qprec_snow, & nsteps, dx, h0, Tsoil, dt, Tsurface0, & Tsurface, G0, lE0, TsurfaceFR, G0FR, lEFR, HFR, qsurface, qevap, qliq, qv, & qyb, qTb, qlyb, qvyb, qlTb, qvTb, qh, qadv, qhyb, qhTb, qadvyb, qadvTb, irec) IMPLICIT NONE INTEGER(i_d), INTENT(IN) :: n TYPE(params), DIMENSION(1:n), INTENT(IN) :: par TYPE(vars_met), INTENT(IN) :: vmet TYPE(vars_snow), INTENT(IN) :: vsnow TYPE(vars), DIMENSION(1:n), INTENT(IN) :: var REAL(r_2), INTENT(IN) :: qprec REAL(r_2), INTENT(IN) :: qprec_snow INTEGER(i_d), INTENT(IN) :: nsteps, irec REAL(r_2), DIMENSION(1:n), INTENT(IN) :: dx REAL(r_2), INTENT(IN) :: h0 REAL(r_2), DIMENSION(1:n), INTENT(IN) :: Tsoil REAL(r_2), INTENT(IN) :: dt REAL(r_2), INTENT(IN) :: Tsurface0 REAL(r_2), INTENT(OUT) :: Tsurface, G0, lE0 ! SEB (subdiurnal, uses T in top layer) REAL(r_2), INTENT(OUT) :: TsurfaceFR, G0FR, lEFR, HFR ! SEB (Force-Restore) REAL(r_2), INTENT(OUT) :: qsurface ! water flux into surface REAL(r_2), INTENT(OUT) :: qevap ! evaporative water flux REAL(r_2), INTENT(OUT) :: qliq, qv ! liquid and vapour components of water flux from surface into soil REAL(r_2), INTENT(OUT) :: qyb, qTb, qlyb, qvyb, qlTb, qvTb ! derivatives of water fluxes wrt moisture and T REAL(r_2), INTENT(OUT) :: qh, qadv ! total and advective components of heat flux into surface REAL(r_2), INTENT(OUT) :: qhyb, qhTb, qadvyb, qadvTb ! derivatives of heat fluxes wrt moiture and T ! local variables INTEGER(i_d) :: surface_case, j REAL(r_2) :: Tsurface_pot, Epot, Hpot, Gpot, dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil REAL(r_2) :: E_vap, dE_vapdT1, E_liq REAL(r_2) :: Kmin, Khmin, phimin REAL(r_2) :: Tqw, dtqwdtb, d1, tmp1d2, Tbar, f, csnow IF (vsnow%nsnow.EQ.0) surface_case = 1 IF (vsnow%nsnow>0) surface_case = 2 SELECT CASE (surface_case) CASE (1) CALL potential_evap(vmet%Rn, vmet%rrc, vmet%rbw, vmet%Ta, vmet%rha, & Tsoil(1), var(1)%kth, half*dx(1)+h0, var(1)%lambdav, Tsurface_pot, Epot, Hpot, & Gpot, dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil) IF (var(1)%iice.EQ.1.AND.Tsurface_pot> zero) THEN Tsurface_pot = 0.0_r_2 Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) vmet%cva)*rhow*var(1)%lambdav/vmet%rbw dEdTsoil = zero dGdTsoil = zero Hpot = rhocp*(Tsurface - vmet%Ta)/vmet%rrc Gpot = vmet%Rn - Hpot - Epot dEdTs = zero ENDIF IF (var(1)%isat.EQ.1) THEN ! saturated surface =>. potential evporation Tsurface = Tsurface_pot lE0 = Epot G0 = Gpot E_vap = zero dE_vapdT1 = zero E_liq = lE0 ELSE ! unsaturated surface: finite vapour transfer; surface flux may be supply limited IF (var(1)%Dv > 1.e-12_r_2) THEN E_vap = (var(1)%rh*csat(Tsoil(1))-vmet%cva*thousand)/(vmet%rbw + half*dx(1)/var(1)%Dv)*var(1)%lambdav dE_vapdT1 = (var(1)%rh*slope_csat(Tsoil(1)))/(vmet%rbw + half*dx(1)/var(1)%Dv)*var(1)%lambdav ELSE E_vap = zero dE_vapdT1 = zero ENDIF CALL hyofh(hmin, par(1)%lam, par(1)%eta, par(1)%Ke, par(1)%he, & Kmin, Khmin, phimin) ! get phi at hmin E_liq = ((var(1)%phi-phimin)/(half*dx(1))-var(1)%K)*thousand*var(1)%lambdav lE0 = MIN(Epot,E_vap+ E_liq) ! analytic approximation (See Haverd et al. 2013, Appxx) IF (Epot.GT.(E_vap+ E_liq)) dEdTs = zero Tsurface = (-half*dx(1)*lE0 + half*dx(1)*vmet%Rn + & var(1)%kth*Tsoil(1) + half*dx(1)*(one/vmet%rrc*rhocp)*vmet%Ta) & /(var(1)%kth + half*dx(1)*(one/vmet%rrc*rhocp)) IF (var(1)%iice.EQ.1.AND.Tsurface> zero) Tsurface = 0.0_r_2 G0 = var(1)%kth/(half*dx(1))*(Tsurface-Tsoil(1)) dGdTsoil = -var(1)%kth/(half*dx(1)) ENDIF ! write(*,*) var(1)%phi, phimin qevap = lE0/(thousand*var(1)%lambdav) qsurface = qprec + qprec_snow - qevap ! derivatives ! q refers to moisture in numerator ! qh refers to heat in numerator ! y refers to moisture in denominator ! T refers to temperature in denominator ! a refers to the layer above ! b refers to the layer below ! initialise derivatives to zero qyb = zero qTb = zero qhyb = zero qhTb = zero ! liquid and vapour fluxes qlyb = zero qvyb = zero ! potential evap independent of S(1), dependent on T1 IF ((Epot<=(E_vap+ E_liq)).OR.(var(1)%isat.EQ.1.OR.vsnow%nsnow.GT.0)) THEN qyb = zero qTb = -dEdTsoil/(thousand*var(1)%lambdav) qliq = -Epot/(thousand*var(1)%lambdav) qv = zero qlyb = zero qvyb = zero qlTb = qTb qvTb = zero ELSE ! supply limited qTb = -dE_vapdT1/(thousand*var(1)%lambdav) qyb = -(var(1)%phiS/(half*dx(1)) - var(1)%KS) !vh! include vapour component?? qliq = -E_liq/(thousand*var(1)%lambdav) qv = -E_vap/(thousand*var(1)%lambdav) qlyb = -(var(1)%phiS/(half*dx(1)) - var(1)%KS) qvyb = zero qlTb = zero qvTb = qTb ENDIF ! end of partial derivative evaluation ! advective component of heat flux qadv = rhow*cswat*qprec*(vmet%Ta) + rhow*csice*qprec_snow*(MIN(vmet%Ta,zero)) & - rhow*qprec_snow*lambdaf Tqw = MERGE(vmet%Ta, Tsoil(1), (-qevap)>zero) dTqwdTb = MERGE(zero, one, (-qevap)>zero) qadv = qadv + rhow*cswat*Tqw*(-qevap) qadvTb = dTqwdTb + rhow*cswat*Tqw*qTb qadvyb = rhow*cswat*qyb*Tqw qh = qadv + G0 qhyb = qadvyb qhTb = dGdTsoil + qadvTb IF (nsteps.EQ.-1) THEN j = 1 tmp1d2 = zero DO WHILE (tmp1d2.LT.one) ! d1 = (var(j)%kth/(var(j)%csoileff)*86400./pi)**0.5 d1 = SQRT(var(j)%kth/(var(j)%csoileff)*86400./pi) tmp1d2 = tmp1d2 + dx(j)/d1 Tbar = Tsoil(j) j = j+1 ENDDO j = j-1 ! integer corresponding to lowest layer contributing to soil col above damping depth ! damping depth = :sum(dx(k)) (k=1,j-1) +f dx(j) IF (j.GT.1) THEN f = (one - (tmp1d2 -dx(j)/d1))*d1/dx(j) ! fraction of lowest layer contributing to soil col above damping depth d1 = SUM(dx(1:j-1)) + f*dx(j) ! multilayer damping depth ELSE f = d1/dx(1) d1 = f*dx(1) ENDIF WRITE(21,"(2i8,3e16.6, i8, 8e16.6)") irec, j, f, d1, Tbar, vsnow%nsnow,Tsurface0, vmet%Rn, lE0, dEdTs, & vsnow%depth(1)!, vsnow%depth(2) CALL forcerestore(Tsurface0, vmet%Rn, lE0, dEdTs, vmet%Ta, & Tbar, d1, vmet%rrc, var(1)%kth, & var(1)%csoileff, dt, var(1)%iice, TsurfaceFR, G0FR, HFR, lEFR) ! write(37,"(16e16.6)") Tsurface0, vmet%Rn, lE0, dEdTs, vmet%Ta, & ! Tbar, d1, vmet%rrc, var(1)%csoileff, dt, real(var(1)%iice), TsurfaceFR, G0FR, HFR, lEFR ENDIF CASE (2) !dedicated snow layer ! SEB at snow/air interface IF (vsnow%hliq(1)>zero) THEN Tsurface = 0.0_r_2 Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) vmet%cva)*rhow*lambdaf/vmet%rbw dEdTsoil = zero dGdTsoil = zero Hpot = rhocp*(Tsurface - vmet%Ta)/vmet%rrc Gpot = vmet%Rn - Hpot - Epot dEdTs = zero qevap = Epot/(thousand*lambdaf) qTb = zero ELSE CALL potential_evap(vmet%Rn, vmet%rrc, vmet%rbw, vmet%Ta, vmet%rha, & vsnow%tsn(1), vsnow%kth(1), half*vsnow%depth(1), & lambdas, Tsurface, Epot, Hpot, & Gpot, dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil) IF (Tsurface > zero) THEN ! temperature of frozen surface must be <= zero Tsurface = 0.0_r_2 Epot = (esat(Tsurface)*0.018_r_2/thousand/8.314_r_2/(vmet%Ta+Tzero) - & ! m3 H2O (liq) m-3 (air) vmet%cva)*rhow*lambdas/vmet%rbw dEdTsoil = zero dGdTsoil = zero Hpot = rhocp*(Tsurface - vmet%Ta)/vmet%rrc Gpot = vmet%Rn - Hpot - Epot dEdTs= zero ENDIF qevap = Epot/(thousand*lambdas) qTb = -dEdTsoil/(thousand*lambdas) ENDIF ! moisture flux at air/snow interface qsurface = qprec_snow+qprec-qevap qyb = zero ! conductive heat flux at air/snow interface G0 = Gpot qhTb = dGdTsoil IF (vsnow%hliq(1)>zero) THEN qhTb = zero qTb = zero ENDIF ! advective heat flux at air/snow interface qadv = rhow*(qprec_snow)*(csice*(MIN(vmet%Ta,zero))-lambdaf) + & rhow*(qprec)*cswat*(vmet%Ta) Tqw = MERGE(vmet%Ta, vsnow%tsn(1), -qevap>zero) dTqwdTb = MERGE(zero,one, -qevap>zero) IF (vsnow%hliq(vsnow%nsnow)>zero) THEN qadv = qadv + rhow*(-qevap)*cswat*Tqw qadvTb = zero ELSE qadv = qadv + rhow*(-qevap)*cswat*Tqw qadvTb = rhow*cswat*(-qevap)*dTqwdTb + rhow*cswat*Tqw*qTb ENDIF qh = qadv + G0 qhyb = qhyb + qadvyb qhTb = qhTb + qadvTb IF (nsteps.EQ.-1) THEN j = -vsnow%nsnow + 1 tmp1d2 = zero DO WHILE (tmp1d2.LT.one) IF (j.LT.1) THEN !write(*,*) 'chk1', -j+vsnow%nsnow, j, vsnow%nsnow csnow = (csice*(vsnow%hsnow(j+vsnow%nsnow)-vsnow%hliq(j+vsnow%nsnow))+cswat*vsnow%hliq(j+vsnow%nsnow))/ & vsnow%depth(j+vsnow%nsnow)*rhow ! d1 = (vsnow%kth(j+vsnow%nsnow)/(csnow)*86400./pi)**0.5 d1 = SQRT(vsnow%kth(j+vsnow%nsnow)/(csnow)*86400._r_2/pi) tmp1d2 = tmp1d2 + vsnow%depth(j+vsnow%nsnow)/d1 ! check snow index here! Tbar = vsnow%tsn(j+vsnow%nsnow) ELSE ! d1 = (var(j)%kth/(var(j)%csoileff)*86400./pi)**0.5 d1 = SQRT(var(j)%kth/(var(j)%csoileff)*86400._r_2/pi) tmp1d2 = tmp1d2 + dx(j)/d1 Tbar = Tsoil(j) ENDIF j = j+1 ENDDO j = j-1 ! integer corresponding to lowest layer contributing to soil col above damping depth ! damping depth = sum(dx(k)) (k=1,j-1) +f dx(j) IF (j.GT.(-vsnow%nsnow + 1).AND.j.LE.0) THEN ! damping depth within snowpack ! fraction of lowest layer contributing to soil col above damping depth f = (one - (tmp1d2 -vsnow%depth(j+vsnow%nsnow)/d1))*d1/vsnow%depth(j+vsnow%nsnow) d1 = SUM(vsnow%depth(1:j+vsnow%nsnow-1)) + f*vsnow%depth(j+vsnow%nsnow) ! multilayer damping depth ELSEIF (j.GT.0) THEN ! damping depth within soil column f = (one - (tmp1d2 -dx(j)/d1))*d1/dx(j) ! fraction of lowest layer contributing to soil col above damping depth IF (j.GT.1) THEN d1 = SUM(vsnow%depth(1:vsnow%nsnow))+SUM(dx(1:j-1)) + f*dx(j) ! multilayer damping depth ELSE d1 = SUM(vsnow%depth(1:vsnow%nsnow)) + f*dx(j) ! multilayer damping depth ENDIF ELSEIF (j.EQ.-vsnow%nsnow + 1) THEN ! damping depth within top layer of snowpack f = d1/vsnow%depth(j+vsnow%nsnow) d1 = f*vsnow%depth(j+vsnow%nsnow) ENDIF csnow = (csice*(vsnow%hsnow(1)-vsnow%hliq(1))+cswat*vsnow%hliq(1))/vsnow%depth(1)*rhow WRITE(21,"(2i8,3e16.6,i8,8e16.6)") irec,j, f, d1, Tbar, vsnow%nsnow, Tsurface0, vmet%Rn, lE0, dEdTs, & vsnow%depth(1)!, vsnow%depth(2) CALL forcerestore(Tsurface0, vmet%Rn, lE0, dEdTs, vmet%Ta, & Tbar, d1, vmet%rrc, vsnow%kth(1), & csnow, dt, 1, TsurfaceFR, G0FR, HFR, lEFR) ! write(37,"(16e16.6)") Tsurface0, vmet%Rn, lE0, dEdTs, vmet%Ta, & ! Tbar, d1, vmet%rrc, var(1)%csoileff, dt, real(var(1)%iice), TsurfaceFR, G0FR, HFR, lEFR ENDIF END SELECT ! surface_case ! finished all the surfaces END SUBROUTINE SEB_FR !============================================================================= SUBROUTINE generic_thomas_1d(n,A,B,C,r,u,err) USE sli_numbers, ONLY: one IMPLICIT NONE ! in/out INTEGER(i_d), INTENT(IN) :: n REAL(r_2), DIMENSION(1:n,1:2,1:2), INTENT(IN) :: A, B, C REAL(r_2), DIMENSION(1:n,1:2), INTENT(IN) :: r REAL(r_2), DIMENSION(1:n,1:2), INTENT(OUT) :: u INTEGER(i_d), OPTIONAL, INTENT(OUT) :: err ! 0: no error; >0: error ! local REAL(r_2), DIMENSION(1:n,1:2,1:2) :: G REAL(r_2), DIMENSION(1:2,1:2) :: bet REAL(r_2), DIMENSION(1:2) :: d REAL(r_2) :: detbet, detbet1 INTEGER(i_d) :: j IF (PRESENT(err)) err = 0 ! j=1 bet(1:2,1:2) = B(1,1:2,1:2) detbet = bet(1,1)*bet(2,2) - bet(1,2)*bet(2,1) IF (ABS(detbet) < EPSILON(detbet)) THEN WRITE(*,*) 'generic_thomas_1d error1: det = 0' IF (PRESENT(err)) THEN err = 1 RETURN ELSE STOP 1 ENDIF ENDIF detbet1 = one/detbet d(1:2) = r(1,1:2) u(1,1) = (bet(2,2)*d(1) - bet(1,2)*d(2))*detbet1 u(1,2) = (bet(1,1)*d(2) - bet(2,1)*d(1))*detbet1 ! j=2, n DO j=2, n d(1:2) = C(j-1,1:2,1) G(j,1,1) = (bet(2,2)*d(1) - bet(1,2)*d(2))*detbet1 G(j,2,1) = (bet(1,1)*d(2) - bet(2,1)*d(1))*detbet1 d(1:2) = C(j-1,1:2,2) G(j,1,2) = (bet(2,2)*d(1) - bet(1,2)*d(2))*detbet1 G(j,2,2) = (bet(1,1)*d(2) - bet(2,1)*d(1))*detbet1 bet(1:2,1:2) = B(j,1:2,1:2) - MATMUL(A(j,1:2,1:2),G(j,1:2,1:2)) detbet = bet(1,1)*bet(2,2) - bet(1,2)*bet(2,1) IF (ABS(detbet) < EPSILON(detbet)) THEN WRITE(*,*) 'generic_thomas_1d error2: det = 0 at j=', j IF (PRESENT(err)) THEN err = 1 RETURN ELSE STOP 1 ENDIF ENDIF detbet1 = one/detbet d(1:2) = r(j,1:2) - MATMUL(A(j,1:2,1:2),u(j-1,1:2)) u(j,1) = (bet(2,2)*d(1) - bet(1,2)*d(2))*detbet1 u(j,2) = (bet(1,1)*d(2) - bet(2,1)*d(1))*detbet1 END DO ! back substitution DO j=n-1, 1, -1 u(j,1:2) = u(j,1:2) - MATMUL(G(j+1,1:2,1:2),u(j+1,1:2)) END DO ! END SUBROUTINE generic_thomas_1d SUBROUTINE generic_thomas_2d(mp, n,A,B,C,r,u,err) USE sli_numbers, ONLY: one IMPLICIT NONE ! in/out INTEGER(i_d), INTENT(IN) :: mp INTEGER(i_d), INTENT(IN) :: n REAL(r_2), DIMENSION(1:mp,1:n,1:2,1:2), INTENT(IN) :: A, B, C REAL(r_2), DIMENSION(1:mp,1:n,1:2), INTENT(IN) :: r REAL(r_2), DIMENSION(1:mp,1:n,1:2), INTENT(OUT) :: u INTEGER(i_d), OPTIONAL, INTENT(OUT) :: err ! 0: no error; >0: error ! local REAL(r_2), DIMENSION(1:mp,1:n,1:2,1:2) :: G REAL(r_2), DIMENSION(1:mp,1:2,1:2) :: bet REAL(r_2), DIMENSION(1:mp,1:2) :: d REAL(r_2), DIMENSION(1:mp) :: detbet, detbet1 INTEGER(i_d) :: j REAL(r_2), DIMENSION(1:mp,1:2) :: tmp1d REAL(r_2), DIMENSION(1:mp,1:2,1:2) :: tmp2d IF (PRESENT(err)) err = 0 ! j=1 bet(1:mp,1:2,1:2) = B(1:mp,1,1:2,1:2) detbet(1:mp) = bet(1:mp,1,1)*bet(1:mp,2,2) - bet(1:mp,1,2)*bet(1:mp,2,1) IF (ANY(ABS(detbet(1:mp)) < EPSILON(detbet))) THEN WRITE(*,*) 'generic_thomas_2d error1: det = 0' IF (PRESENT(err)) THEN err = 1 RETURN ELSE STOP 1 ENDIF ENDIF detbet1(1:mp) = one/detbet(1:mp) d(1:mp,1:2) = r(1:mp,1,1:2) u(1:mp,1,1) = (bet(1:mp,2,2)*d(1:mp,1) - bet(1:mp,1,2)*d(1:mp,2))*detbet1(1:mp) u(1:mp,1,2) = (bet(1:mp,1,1)*d(1:mp,2) - bet(1:mp,2,1)*d(1:mp,1))*detbet1(1:mp) ! j=2, n DO j=2, n d(1:mp,1:2) = C(1:mp,j-1,1:2,1) G(1:mp,j,1,1) = (bet(1:mp,2,2)*d(1:mp,1) - bet(1:mp,1,2)*d(1:mp,2))*detbet1(1:mp) G(1:mp,j,2,1) = (bet(1:mp,1,1)*d(1:mp,2) - bet(1:mp,2,1)*d(1:mp,1))*detbet1(1:mp) d(1:mp,1:2) = C(1:mp,j-1,1:2,2) G(1:mp,j,1,2) = (bet(1:mp,2,2)*d(1:mp,1) - bet(1:mp,1,2)*d(1:mp,2))*detbet1(1:mp) G(1:mp,j,2,2) = (bet(1:mp,1,1)*d(1:mp,2) - bet(1:mp,2,1)*d(1:mp,1))*detbet1(1:mp) tmp2d(1:mp,1,1) = A(1:mp,j,1,1)*G(1:mp,j,1,1) + A(1:mp,j,1,2)*G(1:mp,j,2,1) tmp2d(1:mp,1,2) = A(1:mp,j,1,1)*G(1:mp,j,1,2) + A(1:mp,j,1,2)*G(1:mp,j,2,2) tmp2d(1:mp,2,1) = A(1:mp,j,2,1)*G(1:mp,j,1,1) + A(1:mp,j,2,2)*G(1:mp,j,2,1) tmp2d(1:mp,2,2) = A(1:mp,j,2,1)*G(1:mp,j,1,2) + A(1:mp,j,2,2)*G(1:mp,j,2,2) bet(1:mp,1:2,1:2) = B(1:mp,j,1:2,1:2) - tmp2d(1:mp,1:2,1:2) detbet(1:mp) = bet(1:mp,1,1)*bet(1:mp,2,2) - bet(1:mp,1,2)*bet(1:mp,2,1) IF (ANY(ABS(detbet(1:mp)) < EPSILON(detbet))) THEN WRITE(*,*) 'generic_thomas_2d error2: det = 0 at j=', j IF (PRESENT(err)) THEN err = 1 RETURN ELSE STOP 1 ENDIF ENDIF detbet1(1:mp) = one/detbet(1:mp) tmp1d(1:mp,1) = A(1:mp,j,1,1)*u(1:mp,j-1,1) + A(1:mp,j,1,2)*u(1:mp,j-1,2) tmp1d(1:mp,2) = A(1:mp,j,2,1)*u(1:mp,j-1,1) + A(1:mp,j,2,2)*u(1:mp,j-1,2) d(1:mp,1:2) = r(1:mp,j,1:2) - tmp1d(1:mp,1:2) u(1:mp,j,1) = (bet(1:mp,2,2)*d(1:mp,1) - bet(1:mp,1,2)*d(1:mp,2))*detbet1(1:mp) u(1:mp,j,2) = (bet(1:mp,1,1)*d(1:mp,2) - bet(1:mp,2,1)*d(1:mp,1))*detbet1(1:mp) END DO ! back substitution DO j=n-1, 1, -1 tmp1d(1:mp,1) = G(1:mp,j+1,1,1)*u(1:mp,j+1,1) + G(1:mp,j+1,1,2)*u(1:mp,j+1,2) tmp1d(1:mp,2) = G(1:mp,j+1,2,1)*u(1:mp,j+1,1) + G(1:mp,j+1,2,2)*u(1:mp,j+1,2) u(1:mp,j,1:2) = u(1:mp,j,1:2) - tmp1d(1:mp,1:2) END DO ! END SUBROUTINE generic_thomas_2d !============================================================================= SUBROUTINE getfluxes_vp_1d(n, dx, vtop, vbot, parin, var, hint, phimin, q, qya, qyb, qTa, qTb, & ql, qlya, qlyb, qv, qvT, qvh, qvya, qvyb, & iflux, init, getq0, getqn, Tsoil, T0, nsat, nsatlast) IMPLICIT NONE INTEGER(i_d), INTENT(IN) :: n REAL(r_2), DIMENSION(1:n), INTENT(IN) :: dx TYPE(vars), INTENT(IN) :: vtop TYPE(vars), INTENT(IN) :: vbot TYPE(params), DIMENSION(1:n), INTENT(IN) :: parin TYPE(vars), DIMENSION(1:n), INTENT(INOUT) :: var REAL(r_2), DIMENSION(1:n), INTENT(INOUT) :: hint REAL(r_2), DIMENSION(1:n), INTENT(INOUT) :: phimin REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: q REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qya REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qyb REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qTa REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qTb REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: ql REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qlya REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qlyb REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qv REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qvT REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qvh REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qvya REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qvyb INTEGER(i_d), INTENT(IN) :: iflux LOGICAL, INTENT(IN) :: init LOGICAL, INTENT(IN) :: getq0 LOGICAL, INTENT(IN) :: getqn REAL(r_2), DIMENSION(1:n), INTENT(IN) :: Tsoil REAL(r_2), INTENT(IN) :: T0 INTEGER(i_d), INTENT(IN) :: nsat INTEGER(i_d), INTENT(IN) :: nsatlast ! Gets fluxes q and partial derivs qya, qyb wrt S (if unsat) or phi (if sat). ! Fluxes at top and bottom of profile, and fluxes due to plant extraction of ! water are included. ! Definitions of arguments: ! k - land point ! n - no. of soil layers. ! jt(1:n) - layer soil type nos. ! dx(1:n) - layer thicknesses. ! dz(1:n-1) - distances between layer centres. ! vtop - water vars at soil surface. ! vbot - water vars at bottom of profile. ! var(1:n) - water vars at layer centres. ! hint(1:n) - values of h at interfaces are stored sequentially in hint. ! phimin(1:n) - similarly for phi at hmin in layers above interfaces. ! q(0:n) - fluxes; q(i), i=1,...,n-1 is flux from layer i to layer i+1. ! q(0) is surface flux and q(n) is flux at bottom of profile. ! qya(0:n) - partial deriv of q(i), i=0,...,n, wrt the variable to be solved ! for (S, phi or h) at upper end of flow path. ! qyb(0:n) - ditto for var at lower end. ! iflux - if iflux/=1, get only fluxes involving sat layers. ! init - true if hint and phimin to be initialised. ! getq0 - true if q(0) required. ! getqn - true if q(n) required. LOGICAL :: flag, limit INTEGER(i_d) :: i, itmp, l REAL(r_2) :: dphii1, dhi, h1, h2, hi, Khi1, Khi2, phii1, q2, qya2, qyb2, y, y1, y2 REAL(r_2) :: qTa2, qTb2 TYPE(vars) :: vi1, vi2 REAL(r_2), DIMENSION(1:n-1) :: dz dz(:) = half*(dx(1:n-1)+dx(2:n)) vi1 = zerovars() vi2 = zerovars() IF ((iflux==1) .OR. (var(1)%isat /= 0)) THEN ! get top flux if required IF (getq0) THEN CALL flux(parin(1), vtop, var(1), half*dx(1), q(0), qya(0), qyb(0), qTa(0), qTb(0)) q(0) = q(0)+(T0-Tsoil(1))*(var(1)%kE)/thousand/var(1)%lambdav/dx(1)*two qTa(0) = zero qTb(0) = -(var(1)%kE)/thousand/var(1)%lambdav/dx(1)*two qv(0) = (vtop%phiv-var(1)%phiv)/dx(1)*two +(T0-Tsoil(1))*(var(1)%kE)/thousand/var(1)%lambdav/dx(1)*two ql(0) = q(0) - qv(0) qvya(0) = zero IF (vtop%isat==0) THEN qvyb(0) = vtop%phivS/dx(1)*two ELSE qvyb(0) = zero END IF qlya(0) = qya(0) - qvya(0) qlyb(0) = qyb(0) - qvyb(0) END IF END IF ! otherwise undefined qvh(0) = zero qvT(0) = zero ! get other fluxes l = 0 DO i=1, n-1 IF (iflux==1 .OR. var(i)%isat/=0 .OR. var(i+1)%isat/=0 .OR. nsat/=nsatlast) THEN ! get flux IF (parin(i)%ishorizon == parin(i+1)%ishorizon) THEN ! same soil type, no interface CALL flux(parin(i), var(i), var(i+1), dz(i), q(i), qya(i), qyb(i), qTa(i), qTb(i)) ELSE ! interface l = l+1 IF (init) THEN ! initialise CALL hyofh(hmin, parin(i)%lam, parin(i)%eta, parin(i)%Ke, parin(i)%he, vi1%K, Khi1, phimin(l)) ! get phi at hmin h1 = var(i)%h h2 = var(i+1)%h y1 = var(i)%K*dx(i+1) y2 = var(i+1)%K*dx(i) ! equate fluxes (K constant) to get initial estimate of h at interface hint(l) = (y1*h1+y2*h2+half*gf*(var(i)%K-var(i+1)%K)*dx(i)*dx(i+1))/(y1+y2) END IF hi = hint(l) flag = .TRUE. itmp = 0 ! iterate to get hi at interface for equal fluxes using Newton's method ! get dphii1 at interface in upper layer, because of better linearity, ! then convert to dhi DO WHILE (flag) itmp = itmp+1 IF (itmp>100) THEN !write(*,*) "getfluxes: too many iterations finding interface h" !stop CALL flux(parin(i), var(i), var(i+1), dz(i), q(i), qya(i), qyb(i), qTa(i), qTb(i)) GOTO 111 END IF IF (hi<parin(i)%he) THEN vi1%isat = 0 CALL hyofh(hi, parin(i)%lam, parin(i)%eta, parin(i)%Ke, parin(i)%he, vi1%K, Khi1, phii1) vi1%KS = Khi1/vi1%K ! use dK/dphi, not dK/dS ELSE vi1%isat = 1 vi1%K = var(i)%Ksat phii1 = var(i)%phie+(hi-parin(i)%he)*var(i)%Ksat vi1%KS = zero END IF vi1%h = hi vi1%phi = phii1 vi1%phiS = one ! use dphi/dphi not dphi/dS ! define phiT=0, KT=0 to be consistent with undefined version vi1%phiT = zero vi1%KT = zero ! macropore_factor was not defined but is used in flux(), set to factor of upper layer vi1%macropore_factor = var(i)%macropore_factor CALL flux(parin(i), var(i), vi1, half*dx(i), q(i), qya(i), qyb(i), qTa(i), qTb(i)) IF (hi<parin(i+1)%he) THEN vi2%isat = 0 CALL hyofh(hi, parin(i+1)%lam, parin(i+1)%eta, parin(i+1)%Ke, parin(i+1)%he, vi2%K, Khi2, vi2%phi) vi2%KS = Khi2/vi2%K ! dK/dphi ELSE vi2%isat = 1 vi2%K = var(i+1)%Ksat vi2%phi = var(i+1)%phie+(hi-parin(i+1)%he)*var(i+1)%Ksat END IF vi2%h = hi vi2%phiS = one ! dphi/dphi ! define phiT=0, KT=0 to be consitent with undefined version vi2%phiT = zero vi2%KT = zero ! macropore_factor was not defined but is used in flux(), set to factor of lower layer vi2%macropore_factor = var(i+1)%macropore_factor CALL flux(parin(i+1), vi2, var(i+1), half*dx(i+1), q2, qya2, qyb2, qTa2, qTb2) qya2 = qya2*vi2%K/vi1%K ! partial deriv wrt phii1 ! adjust for equal fluxes dphii1 = -(q(i)-q2)/(qyb(i)-qya2) limit = .FALSE. IF (phii1+dphii1<=phimin(l)) THEN ! out of range limit = .TRUE. dphii1 = -half*(phii1-phimin(l)) END IF phii1 = phii1+dphii1 dhi = dphii1/(vi1%K+half*vi1%KS*dphii1) ! 2nd order Pade approx IF (-vi1%KS*dphii1 > 1.5_r_2*vi1%K) THEN ! use 1st order approx for dhi dhi = dphii1/vi1%K END IF hi = hi+dhi ! check for convergence - dphi/(mean phi)<=dpmaxr IF (.NOT.(limit .OR. ABS(dphii1/(phii1-half*dphii1))>dpmaxr)) THEN flag = .FALSE. END IF END DO ! while flag q(i) = q(i) + qyb(i)*dphii1 hint(l) = hi ! adjust derivs y = one/(qya2-qyb(i)) qya(i) = qya(i)*qya2*y qyb(i) = -qyb2*qyb(i)*y END IF END IF 111 ql(i) = q(i) qTa(i) = qTa(i)+(var(i)%kE+var(i+1)%kE)/thousand/var(1)%lambdav/two/dz(i) qTb(i) = qTb(i)-(var(i)%kE+var(i+1)%kE)/thousand/var(1)%lambdav/two/dz(i) qvT(i) = (Tsoil(i)-Tsoil(i+1))*(var(i)%kE+var(i+1)%kE)/thousand/var(1)%lambdav/two/dz(i) !MC The full description (next two lines) gave problems before -> third line ! Try again original !VH Do both formulations come to the same? I found that I could not reproduce ! Barnes and Allison semi-analytic solution with original !MC This should be re-checked ! qvh(i) = ((((Tsoil(i)+Tzero)/Tzero)**1.88+((Tsoil(i+1)+Tzero)/Tzero)**1.88)/two) & ! * ((var(i)%cvsat+var(i+1)%cvsat)/two)*(var(i)%phiv-var(i+1)%phiv)/dz(i) qvh(i) = (half*(EXP(1.88_r_2*LOG((Tsoil(i)+Tzero)/Tzero))+EXP(1.88_r_2*LOG((Tsoil(i+1)+Tzero)/Tzero)))) & * (half*(var(i)%cvsat+var(i+1)%cvsat)) * (var(i)%phiv-var(i+1)%phiv)/dz(i) ! qvh(i) = ((var(i)%Dv+var(i+1)%Dv)/two)* ((var(i)%cvsat+var(i+1)%cvsat)/two)*(var(i)%rh-var(i+1)%rh)/dz(i) qv(i) = qvh(i) + qvT(i) ! whole vapour flux has one part from humidity (qvh) and one part from temp diff (qvT) q(i) = qv(i) + ql(i) IF (var(i)%isat==0) THEN ! qvya(i) = var(i)%phivS/dz(i) *((((Tsoil(i)+Tzero)/Tzero)**1.88_r_2+ & ! ((Tsoil(i+1)+Tzero)/Tzero)**1.88_r_2)/two) & ! * ((var(i)%cvsat+var(i+1)%cvsat)/two) qvya(i) = var(i)%phivS/dz(i) *(half*(EXP(1.88_r_2*LOG((Tsoil(i)+Tzero)/Tzero)) + & EXP(1.88_r_2*LOG((Tsoil(i+1)+Tzero)/Tzero)))) & * (half*(var(i)%cvsat+var(i+1)%cvsat)) ELSE qvya(i) = zero END IF IF (var(i)%isat==0) THEN ! qvyb(i) = -var(i+1)%phivS/dz(i) *((((Tsoil(i)+Tzero)/Tzero)**1.88_r_2+ & ! ((Tsoil(i+1)+Tzero)/Tzero)**1.88_r_2)/two) & ! * ((var(i)%cvsat+var(i+1)%cvsat)/two) qvyb(i) = -var(i+1)%phivS/dz(i) *(half*(EXP(1.88_r_2*LOG((Tsoil(i)+Tzero)/Tzero)) + & EXP(1.88_r_2*LOG((Tsoil(i+1)+Tzero)/Tzero)))) & * (half*(var(i)%cvsat+var(i+1)%cvsat)) ELSE qvyb(i) = zero END IF qlya(i) = qya(i) qlyb(i) = qyb(i) qya(i) = qya(i) + qvya(i) qyb(i) = qyb(i) + qvyb(i) END DO IF (iflux==1 .OR. var(n)%isat/=0) THEN ! get bottom flux if required IF (getqn) THEN CALL flux(parin(n), var(n), vbot, half*dx(n), q(n), qya(n), qyb(n), qTa(n), qTb(n)) qvya(n) = zero qvyb(n) = zero qlya(n) = qya(n) qlyb(n) = zero ELSE qvya(n) = zero qvyb(n) = zero qlya(n) = qya(n) qlyb(n) = zero END IF ELSE qvya(n) = zero qvyb(n) = zero qlya(n) = qya(n) qlyb(n) = zero END IF ! otherwise undefined ql(n) = q(n) qv(n) = zero qvh(n) = zero qvT(n) = zero DO i=1, n-1 IF (var(i)%Dv == zero .OR. var(i+1)%Dv == zero) THEN q(i) = q(i) - qv(i) qya(i) = qya(i) - qvya(i) qyb(i) = qyb(i) - qvyb(i) qv(i) = zero !qTa(i) = zero !qTb(i) = zero qvya(i) = zero qvyb(i) = zero ENDIF ENDDO END SUBROUTINE getfluxes_vp_1d SUBROUTINE getfluxes_vp_2d(dx, vtop, vbot, parin, var, hint, phimin, i_q, i_qya, i_qyb, i_qTa, i_qTb, & i_ql, i_qlya, i_qlyb, i_qv, i_qvT, i_qvh, i_qvya, i_qvyb, iflux, init, getq0, getqn, Tsoil, T0, nsat, nsatlast) IMPLICIT NONE REAL(r_2), DIMENSION(:,:), INTENT(IN) :: dx ! 1:n TYPE(vars), DIMENSION(:), INTENT(IN) :: vtop TYPE(vars), DIMENSION(:), INTENT(IN) :: vbot TYPE(params), DIMENSION(:,:), INTENT(IN) :: parin ! 1:n TYPE(vars), DIMENSION(:,:), INTENT(IN) :: var ! 1:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: hint ! 1:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: phimin ! 1:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_q ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qya ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qyb ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qTa ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qTb ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_ql ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qlya ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qlyb ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qv ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qvT ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qvh ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qvya ! 0:n REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: i_qvyb ! 0:n INTEGER(i_d), DIMENSION(:), INTENT(IN) :: iflux LOGICAL, DIMENSION(:), INTENT(IN) :: init LOGICAL, DIMENSION(:), INTENT(IN) :: getq0 LOGICAL, DIMENSION(:), INTENT(IN) :: getqn REAL(r_2), DIMENSION(:,:), INTENT(IN) :: Tsoil ! 1:n REAL(r_2), DIMENSION(:), INTENT(IN) :: T0 INTEGER(i_d), DIMENSION(:), INTENT(IN) :: nsat INTEGER(i_d), DIMENSION(:), INTENT(IN) :: nsatlast ! Gets fluxes q and partial derivs qya, qyb wrt S (if unsat) or phi (if sat). ! Fluxes at top and bottom of profile, and fluxes due to plant extraction of ! water are included. ! Definitions of arguments: ! k - land point ! n - no. of soil layers. ! jt(1:n) - layer soil type nos. ! dx(1:n) - layer thicknesses. ! dz(1:n-1) - distances between layer centres. ! vtop - water vars at soil surface. ! vbot - water vars at bottom of profile. ! var(1:n) - water vars at layer centres. ! hint(1:n) - values of h at interfaces are stored sequentially in hint. ! phimin(1:n) - similarly for phi at hmin in layers above interfaces. ! q(0:n) - fluxes; q(i), i=1,...,n-1 is flux from layer i to layer i+1. ! q(0) is surface flux and q(n) is flux at bottom of profile. ! qya(0:n) - partial deriv of q(i), i=0,...,n, wrt the variable to be solved ! for (S, phi or h) at upper end of flow path. ! qyb(0:n) - ditto for var at lower end. ! iflux - if iflux/=1, get only fluxes involving sat layers. ! init - true if hint and phimin to be initialised. ! getq0 - true if q(0) required. ! getqn - true if q(n) required. LOGICAL, DIMENSION(1:SIZE(dx,1)) :: limit, l1, l2, l3 REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: dphii1, dhi, h1, h2, hi, Khi1, Khi2, phii1 REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: q2, qya2, qyb2, y, y1, y2 REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: htmp REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: ztmp1, ztmp2, ztmp3, ztmp4, ztmp5 TYPE(vars), DIMENSION(1:SIZE(dx,1)) :: vi1, vi2 REAL(r_2), DIMENSION(1:SIZE(dx,1),1:SIZE(dx,2)-1) :: dz REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: q REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qya REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qyb REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qTa REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qTb REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: ql REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qlya REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qlyb REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qv REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qvT REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qvh REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qvya REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qvyb TYPE(vars) :: vtmp INTEGER(i_d) :: i, n, mp, itmp mp = SIZE(dx,1) n = SIZE(dx,2) dz(:,1:n-1) = half*(dx(:,1:n-1)+dx(:,2:n)) q(:,0:n) = i_q(:,1:n+1) qya(:,0:n) = i_qya(:,1:n+1) qyb(:,0:n) = i_qyb(:,1:n+1) qTa(:,0:n) = i_qTa(:,1:n+1) qTb(:,0:n) = i_qTb(:,1:n+1) ql(:,0:n) = zero qlya(:,0:n) = zero qlyb(:,0:n) = zero qv(:,0:n) = zero qvT(:,0:n) = zero qvh(:,0:n) = zero qvya(:,0:n) = zero qvyb(:,0:n) = zero vtmp = zerovars() vi1 = SPREAD(vtmp,1,mp) vi2 = SPREAD(vtmp,1,mp) ztmp1(:) = zero ztmp2(:) = zero ztmp3(:) = zero ztmp4(:) = zero ztmp5(:) = zero l1(:) = ((iflux(:)==1) .OR. (var(:,1)%isat /= 0)) .AND. getq0(:) IF (ANY(l1(:))) & CALL flux(parin(:,1), vtop(:), var(:,1), half*dx(:,1), ztmp1(:), ztmp2(:), ztmp3(:), ztmp4(:), ztmp5(:)) WHERE (l1(:)) ! get top flux if required q(:,0) = ztmp1(:) qya(:,0) = ztmp2(:) qyb(:,0) = ztmp3(:) qTa(:,0) = ztmp4(:) qTb(:,0) = ztmp5(:) q(:,0) = q(:,0)+(T0(:)-Tsoil(:,1))*(var(:,1)%kE)/thousand/var(:,1)%lambdav/dx(:,1)*two qTa(:,0) = zero qTb(:,0) = -(var(:,1)%kE)/thousand/var(:,1)%lambdav/dx(:,1)*two qv(:,0) = (vtop(:)%phiv-var(:,1)%phiv)/dx(:,1)*two +(T0(:)-Tsoil(:,1))* & (var(:,1)%kE)/thousand/var(:,1)%lambdav/dx(:,1)*two ql(:,0) = q(:,0) - qv(:,0) qvya(:,0) = zero qvyb(:,0) = zero ! isat==0 below qlya(:,0) = qya(:,0) - qvya(:,0) qlyb(:,0) = qyb(:,0) - qvyb(:,0) endwhere WHERE (l1(:) .AND. (vtop(:)%isat==0)) qvyb(:,0) = vtop(:)%phivS/dx(:,1)*two qlyb(:,0) = qyb(:,0) - qvyb(:,0) endwhere ! otherwise undefined qvh(:,0) = zero qvT(:,0) = zero ! get other fluxes DO i=1, n-1 l1(:) = (iflux(:)==1 .OR. var(:,i)%isat/=0 .OR. var(:,i+1)%isat/=0 .OR. nsat(:)/=nsatlast(:)) IF (ANY(l1(:) .AND. (parin(:,i)%ishorizon==parin(:,i+1)%ishorizon))) & CALL flux(parin(:,i), var(:,i), var(:,i+1), dz(:,i), ztmp1(:), ztmp2(:), ztmp3(:), ztmp4(:), ztmp5(:)) WHERE (l1(:) .AND. (parin(:,i)%ishorizon==parin(:,i+1)%ishorizon)) ! same soil type, no interface q(:,i) = ztmp1(:) qya(:,i) = ztmp2(:) qyb(:,i) = ztmp3(:) qTa(:,i) = ztmp4(:) qTb(:,i) = ztmp5(:) endwhere l2(:) = l1(:) .AND. (parin(:,i)%ishorizon /= parin(:,i+1)%ishorizon) ! interface htmp(:) = hmin IF (ANY(l2(:) .AND. init(:))) & CALL hyofh(htmp(:), parin(:,i)%lam, parin(:,i)%eta, parin(:,i)%Ke, parin(:,i)%he, ztmp1(:), ztmp2(:), ztmp3(:)) WHERE (l2(:) .AND. init(:)) ! interface & get phi at hmin vi1(:)%K = ztmp1(:) Khi1(:) = ztmp2(:) phimin(:,i) = ztmp3(:) h1(:) = var(:,i)%h h2(:) = var(:,i+1)%h y1(:) = var(:,i)%K*dx(:,i+1) y2(:) = var(:,i+1)%K*dx(:,i) ! equate fluxes (K constant) to get initial estimate of h at interface hint(:,i) = (y1(:)*h1(:) + y2(:)*h2(:) + half*gf*(var(:,i)%K-var(:,i+1)%K)*dx(:,i)*dx(:,i+1)) / (y1(:)+y2(:)) endwhere !where ((.not. l2(:)) .and. init(:)) hint(:,i) = zero hi(:) = hint(:,i) ! iterate to get hi at interface for equal fluxes using Newton's method ! get dphii1 at interface in upper layer, because of better linearity, ! then convert to dhi IF (ANY(l2(:))) THEN l3(:) = l2(:) limit(:) = .FALSE. DO itmp=1, 100 IF (ANY(l3(:) .AND. (hi(:)<parin(:,i)%he))) & CALL hyofh(hi(:), parin(:,i)%lam, parin(:,i)%eta, parin(:,i)%Ke, parin(:,i)%he, ztmp1(:), ztmp2(:), ztmp3(:)) WHERE (l3(:) .AND. (hi(:)<parin(:,i)%he)) vi1(:)%isat = 0 vi1(:)%K = ztmp1(:) Khi1(:) = ztmp2(:) phii1(:) = ztmp3(:) vi1(:)%KS = Khi1(:)/vi1(:)%K ! use dK/dphi, not dK/dS endwhere WHERE (l3(:) .AND. (hi(:)>=parin(:,i)%he)) vi1(:)%isat = 1 vi1(:)%K = var(:,i)%Ksat phii1(:) = var(:,i)%phie+(hi(:)-parin(:,i)%he)*var(:,i)%Ksat vi1(:)%KS = zero endwhere WHERE (l3(:)) vi1(:)%h = hi(:) vi1(:)%phi = phii1(:) vi1(:)%phiS = one ! use dphi/dphi not dphi/dS !define phiT=0, KT=0 to be consitent with undefined version vi1(:)%phiT = zero vi1(:)%KT = zero ! macropore_factor was not defined but is used in flux(), set to factor of upper layer vi1(:)%macropore_factor = var(:,i)%macropore_factor endwhere IF (ANY(l3(:))) & CALL flux(parin(:,i), var(:,i), vi1(:), half*dx(:,i), ztmp1(:), ztmp2(:), ztmp3(:), ztmp4(:), ztmp5(:)) WHERE (l3(:)) q(:,i) = ztmp1(:) qya(:,i) = ztmp2(:) qyb(:,i) = ztmp3(:) qTa(:,i) = ztmp4(:) qTb(:,i) = ztmp5(:) endwhere IF (ANY(l3(:) .AND. (hi(:)<parin(:,i+1)%he))) & CALL hyofh(hi(:), parin(:,i+1)%lam, parin(:,i+1)%eta, parin(:,i+1)%Ke, parin(:,i+1)%he, & ztmp1(:), ztmp2(:), ztmp3(:)) WHERE (l3(:) .AND. (hi(:)<parin(:,i+1)%he)) vi2(:)%K = ztmp1(:) Khi2(:) = ztmp2(:) vi2(:)%phi = ztmp3(:) vi2(:)%isat = 0 vi2(:)%KS = Khi2(:)/vi2(:)%K ! dK/dphi endwhere WHERE (l3(:) .AND. (hi(:)>=parin(:,i+1)%he)) vi2(:)%isat = 1 vi2(:)%K = var(:,i+1)%Ksat vi2(:)%phi = var(:,i+1)%phie+(hi(:)-parin(:,i+1)%he)*var(:,i+1)%Ksat endwhere WHERE (l3(:)) vi2(:)%h = hi(:) vi2(:)%phiS = one ! dphi/dphi ! define phiT=0, KT=0 to be consitent with undefined version vi2(:)%phiT = zero vi2(:)%KT = zero ! macropore_factor was not defined but is used in flux(), set to factor of lower layer vi2(:)%macropore_factor = var(:,i+1)%macropore_factor endwhere IF (ANY(l3(:))) & CALL flux(parin(:,i+1), vi2(:), var(:,i+1), half*dx(:,i+1), ztmp1(:), ztmp2(:), ztmp3(:), ztmp4(:), ztmp5(:)) WHERE (l3(:)) q2(:) = ztmp1(:) qya2(:) = ztmp2(:) qyb2(:) = ztmp3(:) qya2(:) = qya2(:)*vi2(:)%K/vi1(:)%K ! partial deriv wrt phii1 ! adjust for equal fluxes dphii1(:) = -(q(:,i)-q2(:))/(qyb(:,i)-qya2(:)) limit(:) = .FALSE. endwhere WHERE (l3(:) .AND. (phii1(:)+dphii1(:)<=phimin(:,i))) ! out of range limit(:) = .TRUE. dphii1(:) = -half*(phii1(:)-phimin(:,i)) endwhere WHERE (l3(:)) phii1(:) = phii1(:)+dphii1(:) dhi(:) = dphii1(:)/(vi1(:)%K+half*vi1(:)%KS*dphii1(:)) ! 2nd order Pade approx endwhere WHERE (l3(:) .AND. (-vi1%KS*dphii1 > 1.5_r_2*vi1%K)) ! use 1st order approx for dhi dhi(:) = dphii1(:)/vi1(:)%K endwhere WHERE (l3(:)) hi(:) = hi(:)+dhi(:) endwhere ! check for convergence - dphi/(mean phi)<=dpmaxr WHERE (l3(:) .AND. & .NOT. (limit(:) .OR. (ABS(dphii1(:)/(phii1(:)-half*dphii1(:)))>dpmaxr))) l3(:) = .FALSE. IF (.NOT. ANY(l3(:))) EXIT END DO ! do itmp=1, 100 IF (itmp>=100) THEN !write(*,*) "getfluxes: too many iterations finding interface h" !stop IF (ANY(l2(:) .AND. l3(:))) & CALL flux(parin(:,i), var(:,i), var(:,i+1), dz(:,i), ztmp1(:), ztmp2(:), ztmp3(:), ztmp4(:), ztmp5(:)) WHERE (l2(:) .AND. l3(:)) q(:,i) = ztmp1(:) qya(:,i) = ztmp2(:) qyb(:,i) = ztmp3(:) qTa(:,i) = ztmp4(:) qTb(:,i) = ztmp5(:) endwhere ELSE WHERE (l2(:) .AND. (.NOT. l3(:))) q(:,i) = q(:,i) + qyb(:,i)*dphii1(:) hint(:,i) = hi(:) ! adjust derivs y(:) = one/(qya2(:)-qyb(:,i)) qya(:,i) = qya(:,i)*qya2(:)*y(:) qyb(:,i) = -qyb2(:)*qyb(:,i)*y(:) endwhere END IF END IF ql(:,i) = q(:,i) qTa(:,i) = qTa(:,i)+(var(:,i)%kE+var(:,i+1)%kE)/thousand/var(:,1)%lambdav/two/dz(:,i) qTb(:,i) = qTb(:,i)-(var(:,i)%kE+var(:,i+1)%kE)/thousand/var(:,1)%lambdav/two/dz(:,i) qvT(:,i) = (Tsoil(:,i)-Tsoil(:,i+1))*(var(:,i)%kE+var(:,i+1)%kE)/thousand/var(:,1)%lambdav/two/dz(:,i) !MC The full description (next two lines) gave problems before -> third line ! Try again original !VH Do both formulations come to the same? I found that I could not reproduce ! Barnes and Allison semi-analytic solution with original !MC This should be re-checked ! qvh(:,i) = ((((Tsoil(:,i)+Tzero)/Tzero)**1.88+((Tsoil(:,i+1)+Tzero)/Tzero)**1.88)/two) & ! * ((var(:,i)%cvsat+var(:,i+1)%cvsat)/two)*(var(:,i)%phiv-var(:,i+1)%phiv)/dz(:,i) qvh(:,i) = (half*(EXP(1.88_r_2*LOG((Tsoil(:,i)+Tzero)/Tzero))+EXP(1.88_r_2*LOG((Tsoil(:,i+1)+Tzero)/Tzero)))) & * (half*(var(:,i)%cvsat+var(:,i+1)%cvsat)) * (var(:,i)%phiv-var(:,i+1)%phiv)/dz(:,i) ! qvh(:,i) = ((var(:,i)%Dv+var(:,i+1)%Dv)/two) * & ! ((var(:,i)%cvsat+var(:,i+1)%cvsat)/two)*(var(:,i)%rh-var(:,i+1)%rh)/dz(:,i) qv(:,i) = qvh(:,i) + qvT(:,i) ! whole vapour flux has one part from humidity (qvh) and one part from temp diff (qvT) q(:,i) = qv(:,i) + ql(:,i) WHERE (var(:,i)%isat==0) ! qvya(:,i) = var(:,i)%phivS/dz(:,i) *((((Tsoil(:,i)+Tzero)/Tzero)**1.88_r_2+ & ! ((Tsoil(:,i+1)+Tzero)/Tzero)**1.88_r_2)/two) & ! * ((var(:,i)%cvsat+var(:,i+1)%cvsat)/two) qvya(:,i) = var(:,i)%phivS/dz(:,i) *(half*(EXP(1.88_r_2*LOG((Tsoil(:,i)+Tzero)/Tzero)) + & EXP(1.88_r_2*LOG((Tsoil(:,i+1)+Tzero)/Tzero)))) & * (half*(var(:,i)%cvsat+var(:,i+1)%cvsat)) ELSEWHERE qvya(:,i) = zero endwhere WHERE (var(:,i)%isat==0) ! qvyb(:,i) = -var(:,i+1)%phivS/dz(:,i) *((((Tsoil(:,i)+Tzero)/Tzero)**1.88_r_2+ & ! ((Tsoil(:,i+1)+Tzero)/Tzero)**1.88_r_2)/two) & ! * ((var(:,i)%cvsat+var(:,i+1)%cvsat)/two) qvyb(:,i) = -var(:,i+1)%phivS/dz(:,i) *(half*(EXP(1.88_r_2*LOG((Tsoil(:,i)+Tzero)/Tzero)) + & EXP(1.88_r_2*LOG((Tsoil(:,i+1)+Tzero)/Tzero)))) & * (half*(var(:,i)%cvsat+var(:,i+1)%cvsat)) ELSEWHERE qvyb(:,i) = zero endwhere qlya(:,i) = qya(:,i) qlyb(:,i) = qyb(:,i) qya(:,i) = qya(:,i) + qvya(:,i) qyb(:,i) = qyb(:,i) + qvyb(:,i) END DO l1(:) = (iflux(:)==1) .OR. (var(:,n)%isat/=0) IF (ANY(l1(:) .AND. getqn(:))) & ! get bottom flux if required CALL flux(parin(:,n), var(:,n), vbot(:), half*dx(:,n), ztmp1(:), ztmp2(:), ztmp3(:), ztmp4(:), ztmp5(:)) WHERE (l1(:) .AND. getqn(:)) q(:,n) = ztmp1(:) qya(:,n) = ztmp2(:) qyb(:,n) = ztmp3(:) qTa(:,n) = ztmp4(:) qTb(:,n) = ztmp5(:) qvya(:,n) = zero qvyb(:,n) = zero qlya(:,n) = qya(:,n) qlyb(:,n) = zero endwhere WHERE (l1(:) .AND. (.NOT. getqn(:))) qvya(:,n) = zero qvyb(:,n) = zero qlya(:,n) = qya(:,n) qlyb(:,n) = zero endwhere WHERE (.NOT. l1(:)) qvya(:,n) = zero qvyb(:,n) = zero qlya(:,n) = qya(:,n) qlyb(:,n) = zero endwhere ! otherwise undefined ql(:,n) = q(:,n) qv(:,n) = zero qvh(:,n) = zero qvT(:,n) = zero DO i=1, n-1 WHERE (var(:,i)%Dv == zero .OR. var(:,i+1)%Dv == zero) q(:,i) = q(:,i) - qv(:,i) qya(:,i) = qya(:,i) - qvya(:,i) qyb(:,i) = qyb(:,i) - qvyb(:,i) qv(:,i) = zero !qTa(:,i) = zero !qTb(:,i) = zero qvya(:,i) = zero qvyb(:,i) = zero endwhere ENDDO i_q(:,1:n+1) = q(:,0:n) i_qya(:,1:n+1) = qya(:,0:n) i_qyb(:,1:n+1) = qyb(:,0:n) i_qTa(:,1:n+1) = qTa(:,0:n) i_qTb(:,1:n+1) = qTb(:,0:n) i_ql(:,1:n+1) = ql(:,0:n) i_qlya(:,1:n+1) = qlya(:,0:n) i_qlyb(:,1:n+1) = qlyb(:,0:n) i_qv(:,1:n+1) = qv(:,0:n) i_qvT(:,1:n+1) = qvT(:,0:n) i_qvh(:,1:n+1) = qvh(:,0:n) i_qvya(:,1:n+1) = qvya(:,0:n) i_qvyb(:,1:n+1) = qvyb(:,0:n) END SUBROUTINE getfluxes_vp_2d !============================================================================= SUBROUTINE getheatfluxes_1d(n, dx, dxL, qh, qhya, qhyb, qhTa, qhTb, var, vlit, T, TL, litter, & q, qya, qyb, qTa, qTb,qadv, qadvya, qadvyb, qadvTa, qadvTb,advection) ! modified 25/05/10 to include contribution to heat flux from liquid water flux in the presence of ice IMPLICIT NONE INTEGER(i_d), INTENT(IN) :: n REAL(r_2), DIMENSION(1:n), INTENT(IN) :: dx REAL(r_2), INTENT(IN) :: dxL REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qh, q, qadv REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qhya, qya, qadvya REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qhyb, qyb, qadvyb REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qhTa, qTa, qadvTa REAL(r_2), DIMENSION(0:n), INTENT(INOUT) :: qhTb, qTb, qadvTb TYPE(vars), DIMENSION(1:n), INTENT(IN) :: var TYPE(vars), INTENT(IN) :: vlit REAL(r_2), DIMENSION(1:n), INTENT(IN) :: T REAL(r_2), INTENT(IN) :: TL LOGICAL, INTENT(IN) :: litter INTEGER(i_d), INTENT(IN) :: advection ! Gets heat fluxes qh and partial derivs qhya, qhyb wrt T and S (if unsat) or phi (if sat). INTEGER(i_d) :: i REAL(r_2) :: rdz, keff , w REAL(r_2), DIMENSION(1:n-1) :: dz REAL(r_2) :: dTqwdTa, dTqwdTb, Tqw dz(:) = half*(dx(1:n-1)+dx(2:n)) DO i=1, n-1 rdz = one/dz(i) keff = 2_r_2*(var(i)%kth*var(i+1)%kth)/(var(i)%kth*dx(i)+var(i+1)%kth*dx(i+1)) ! qh(i) = keff*(T(i)-T(i+1)) +(var(i)%phiv-var(i+1)%phiv)*var(i)%lambdav*thousand*rdz & ! * ((((T(i)+Tzero) /Tzero)**1.88_r_2+((T(i+1)+Tzero) /Tzero)**1.88_r_2)/two) & ! *((var(i)%cvsat+var(i+1)%cvsat)/two) qh(i) = keff*(T(i)-T(i+1)) +(var(i)%phiv-var(i+1)%phiv)*var(i)%lambdav*thousand*rdz & * (half*(EXP(1.88_r_2*LOG((T(i)+Tzero) /Tzero))+EXP(1.88_r_2*LOG((T(i+1)+Tzero) /Tzero)))) & *(half*(var(i)%cvsat+var(i+1)%cvsat)) IF (var(i)%isat==0) THEN ! qhya(i) = rdz*var(i)%lambdav*thousand*var(i)%phivS*((((T(i)+Tzero) /Tzero)**1.88_r_2+ & ! ((T(i+1)+Tzero) /Tzero)**1.88_r_2)/two) & ! * ((var(i)%cvsat+var(i+1)%cvsat)/two) qhya(i) = rdz*var(i)%lambdav*thousand*var(i)%phivS*(half*(EXP(1.88_r_2*LOG((T(i)+Tzero) /Tzero)) + & EXP(1.88_r_2*LOG((T(i+1)+Tzero) /Tzero)))) & * (half*(var(i)%cvsat+var(i+1)%cvsat)) ELSE qhya(i) = zero END IF IF (var(i+1)%isat==0) THEN ! qhyb(i) = -rdz*var(i)%lambdav*thousand*var(i+1)%phivS*((((T(i)+Tzero) /Tzero)**1.88_r_2+ & ! ((T(i+1)+Tzero) /Tzero)**1.88_r_2)/two) & ! * ((var(i)%cvsat+var(i+1)%cvsat)/two) qhyb(i) = -rdz*var(i)%lambdav*thousand*var(i+1)%phivS*(half*(EXP(1.88_r_2*LOG((T(i)+Tzero) /Tzero)) + & EXP(1.88_r_2*LOG((T(i+1)+Tzero) /Tzero)))) & * (half*(var(i)%cvsat+var(i+1)%cvsat)) ELSE qhyb(i) = zero END IF qhTa(i) = keff qhTb(i) = -keff ! add advective terms IF (advection==1) THEN !$ if (q(i) > zero) then !$ w = (var(i)%kth/dx(i))/(var(i)%kth/dx(i)+var(i+1)%kth/dx(i+1)) !$ else !$ w = (var(i)%kth/dx(i))/(var(i)%kth/dx(i)+var(i+1)%kth/dx(i+1)) !$ endif !$ qadv(i) = rhow*cswat*q(i)*(w*(T(i)+zero)+(one-w)*(T(i+1)+zero)) !$ qadvya(i) = rhow*cswat*qya(i)*(w*(T(i)+zero)+(one-w)*(T(i+1)+zero)) !$ qadvyb(i) = rhow*cswat*qyb(i)*(w*(T(i)+zero)+(one-w)*(T(i+1)+zero)) !$ qadvTa(i) = rhow*cswat*q(i)*w !$ qadvTb(i) = rhow*cswat*q(i)*(one-w) Tqw = MERGE(T(i), T(i+1), q(i)>zero) +zero dTqwdTa = MERGE(one, zero, (q(i)>zero)) dTqwdTb = MERGE(zero,one, q(i)>zero) qadv(i) = rhow*cswat*q(i)*Tqw qadvya(i) = rhow*cswat*qya(i)*Tqw qadvyb(i) = rhow*cswat*qyb(i)*Tqw qadvTa(i) = rhow*cswat*q(i)*dTqwdTa + rhow*cswat*Tqw*qTa(i) qadvTb(i) = rhow*cswat*q(i)*dTqwdTb + rhow*cswat*Tqw*qTb(i) qh(i) = qh(i) + qadv(i) qhya(i) = qhya(i) + qadvya(i) qhyb(i) = qhyb(i) + qadvyb(i) qhTa(i) = qhTa(i) + qadvTa(i) qhTb(i) = qhTb(i) + qadvTb(i) ENDIF ENDDO qh(n) = zero qhya(n) = zero qhyb(n) = zero qhTa(n) = zero qhTb(n) = zero END SUBROUTINE getheatfluxes_1d SUBROUTINE getheatfluxes_2d(dx, dxL, i_qh, i_qhya, i_qhyb, i_qhTa, i_qhTb, var, vlit, T, TL, & litter, i_q,i_qya,i_qyb,i_qTa,i_qTb, & i_qadv,i_qadvya, i_qadvyb, i_qadvTa, i_qadvTb, advection) ! modified 25/05/10 to include contribution to heat flux from liquid water flux in the presence of ice IMPLICIT NONE REAL(r_2), DIMENSION(:,:), INTENT(IN) :: dx ! :,1:n REAL(r_2), DIMENSION(:), INTENT(IN) :: dxL REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qh ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qhya ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qhyb ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qhTa ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qhTb ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_q ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qya ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qyb ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qTa ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qTb ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qadv ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qadvya ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qadvyb ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qadvTa ! :,0:n => :,1:n+1 REAL(r_2), DIMENSION(:,:), INTENT(INOUT) :: i_qadvTb ! :,0:n => :,1:n+1 TYPE(vars), DIMENSION(:,:), INTENT(IN) :: var TYPE(vars), DIMENSION(:), INTENT(IN) :: vlit REAL(r_2), DIMENSION(:,:), INTENT(IN) :: T ! :,1:n REAL(r_2), DIMENSION(:), INTENT(IN) :: TL LOGICAL, INTENT(IN) :: litter INTEGER(i_d), INTENT(IN) :: advection ! Gets heat fluxes qh and partial derivs qhya, qhyb wrt T and S (if unsat) or phi (if sat). INTEGER(i_d) :: i, n REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: rdz REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: keff REAL(r_2), DIMENSION(1:SIZE(dx,1),1:SIZE(dx,2)-1) :: dz REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qh ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qhya ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qhyb ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qhTa ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qhTb ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: q ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qya ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qyb ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qTa ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qTb ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qadv ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qadvya ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qadvyb ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qadvTa ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1),0:SIZE(dx,2)) :: qadvTb ! :,0:n REAL(r_2), DIMENSION(1:SIZE(dx,1)) :: dTqwdTa, dTqwdTb, Tqw n = SIZE(dx,2) dz(:,1:n-1) = half*(dx(:,1:n-1)+dx(:,2:n)) qh(:,0:n) = i_qh(:,1:n+1) qhya(:,0:n) = i_qhya(:,1:n+1) qhyb(:,0:n) = i_qhyb(:,1:n+1) qhTa(:,0:n) = i_qhTa(:,1:n+1) qhTb(:,0:n) = i_qhTb(:,1:n+1) q(:,0:n) = i_q(:,1:n+1) qya(:,0:n) = i_qya(:,1:n+1) qyb(:,0:n) = i_qyb(:,1:n+1) qTa(:,0:n) = i_qTa(:,1:n+1) qTb(:,0:n) = i_qTb(:,1:n+1) qadv(:,0:n) = i_qadv(:,1:n+1) qadvya(:,0:n) = i_qadvya(:,1:n+1) qadvyb(:,0:n) = i_qadvyb(:,1:n+1) qadvTa(:,0:n) = i_qadvTa(:,1:n+1) qadvTb(:,0:n) = i_qadvTb(:,1:n+1) DO i=1, n-1 rdz = one/dz(:,i) keff = 2_r_2*(var(:,i)%kth*var(:,i+1)%kth)/(var(:,i)%kth*dx(:,i)+var(:,i+1)%kth*dx(:,i+1)) ! qh(:,i) = keff*(T(:,i)-T(:,i+1)) & ! + (var(:,i)%phiv-var(:,i+1)%phiv)*var(:,i)%lambdav*thousand*rdz(:) & ! * ((((T(:,i)+Tzero) /Tzero)**1.88_r_2+((T(:,i+1)+Tzero) /Tzero)**1.88_r_2)/two) & ! *((var(:,i)%cvsat+var(:,i+1)%cvsat)/two) qh(:,i) = keff*(T(:,i)-T(:,i+1)) & + (var(:,i)%phiv-var(:,i+1)%phiv)*var(:,i)%lambdav*thousand*rdz(:) & * (half*EXP(1.88_r_2*LOG(((T(:,i)+Tzero) /Tzero)) + EXP(1.88_r_2*LOG((T(:,i+1)+Tzero) /Tzero)))) & *(half*(var(:,i)%cvsat+var(:,i+1)%cvsat)) WHERE (var(:,i)%isat == 0) ! qhya(:,i) = rdz(:)*var(:,i)%lambdav*thousand*var(:,i)%phivS*((((T(:,i)+Tzero) /Tzero)**1.88_r_2+ & ! ((T(:,i+1)+Tzero) /Tzero)**1.88_r_2)/two) & ! * ((var(:,i)%cvsat+var(:,i+1)%cvsat)/two) qhya(:,i) = rdz(:)*var(:,i)%lambdav*thousand*var(:,i)%phivS*(half*(EXP(1.88_r_2*LOG((T(:,i)+Tzero) /Tzero)) + & EXP(1.88_r_2*LOG((T(:,i+1)+Tzero) /Tzero)))) & * (half*(var(:,i)%cvsat+var(:,i+1)%cvsat)) ELSEWHERE qhya(:,i) = zero endwhere WHERE (var(:,i+1)%isat == 0) ! qhyb(:,i) = -rdz(:)*var(:,i)%lambdav*thousand*var(:,i+1)%phivS*((((T(:,i)+Tzero) /Tzero)**1.88_r_2+ & ! ((T(:,i+1)+Tzero) /Tzero)**1.88_r_2)/two) & ! * ((var(:,i)%cvsat+var(:,i+1)%cvsat)/two) qhyb(:,i) = -rdz(:)*var(:,i)%lambdav*thousand*var(:,i+1)%phivS*(half*EXP(1.88_r_2*LOG(((T(:,i)+Tzero) /Tzero)) + & EXP(1.88_r_2*LOG((T(:,i+1)+Tzero) /Tzero)))) & * (half*(var(:,i)%cvsat+var(:,i+1)%cvsat)) ELSEWHERE qhyb(:,i) = zero endwhere qhTa(:,i) = keff qhTb(:,i) = -keff ! add advective terms IF (advection==1) THEN Tqw = MERGE(T(:,i), T(:,i+1), q(:,i)>zero) +zero dTqwdTa = MERGE(one, zero, q(:,i)>zero) dTqwdTb = MERGE(zero,one, q(:,i)>zero) qadv(:,i) = rhow*cswat*q(:,i)*Tqw qadvya(:,i) = rhow*cswat*qya(:,i)*Tqw qadvyb(:,i) = rhow*cswat*qyb(:,i)*Tqw qadvTa(:,i) = rhow*cswat*q(:,i)*dTqwdTa + rhow*cswat*Tqw*qTa(:,i) qadvTb(:,i) = rhow*cswat*q(:,i)*dTqwdTb + rhow*cswat*Tqw*qTb(:,i) qh(:,i) = qh(:,i) + qadv(:,i) qhya(:,i) = qhya(:,i) + qadvya(:,i) qhyb(:,i) = qhyb(:,i) + qadvyb(:,i) qhTa(:,i) = qhTa(:,i) + qadvTa(:,i) qhTb(:,i) = qhTb(:,i) + qadvTb(:,i) ENDIF ENDDO qh(:,n) = zero qhya(:,n) = zero qhyb(:,n) = zero qhTa(:,n) = zero qhTb(:,n) = zero i_qh(:,1:n+1) = qh(:,0:n) i_qhya(:,1:n+1) = qhya(:,0:n) i_qhyb(:,1:n+1) = qhyb(:,0:n) i_qhTa(:,1:n+1) = qhTa(:,0:n) i_qhTb(:,1:n+1) = qhTb(:,0:n) IF (advection==1) THEN i_qadv(:,1:n+1) = qadv(:,0:n) i_qadvya(:,1:n+1) = qadvya(:,0:n) i_qadvyb(:,1:n+1) = qadvyb(:,0:n) i_qadvTa(:,1:n+1) = qadvTa(:,0:n) i_qadvTb(:,1:n+1) = qadvTb(:,0:n) ENDIF END SUBROUTINE getheatfluxes_2d !============================================================================= ELEMENTAL PURE SUBROUTINE hyofh(h, lam, eta, Ke, he, K, Kh, phi) IMPLICIT NONE REAL(r_2), INTENT(IN) :: h REAL(r_2), INTENT(IN) :: lam REAL(r_2), INTENT(IN) :: eta REAL(r_2), INTENT(IN) :: Ke REAL(r_2), INTENT(IN) :: he REAL(r_2), INTENT(OUT) :: K REAL(r_2), INTENT(OUT) :: Kh REAL(r_2), INTENT(OUT) :: phi ! Get soil water variables from h. ! Definitions of arguments: ! h - matric head. ! K - hydraulic conductivity. ! Kh - derivative dK/dh. ! phi - matric flux potential (MFP). REAL(r_2) :: a a = -lam * eta K = Ke * EXP(a*LOG(h/he)) Kh = a * K / h phi = K * h / (one+a) END SUBROUTINE hyofh !============================================================================= ! For debug: remove elemental pure ELEMENTAL PURE SUBROUTINE hyofS(S, Tsoil, parin, var) ! SUBROUTINE hyofS(S, Tsoil, parin, var) IMPLICIT NONE REAL(r_2), INTENT(IN) :: S REAL(r_2), INTENT(IN) :: Tsoil TYPE(params), INTENT(INOUT) :: parin TYPE(vars), INTENT(INOUT) :: var ! Get soil water variables from S. ! Definitions of arguments: ! S(1:ms) - degree of saturation ("effective satn") of layers. ! ms - no. of soil layers. ! jt(1:ms) - layer soil type nos. ! var(1:ms) - other water vars of layers. REAL(r_2) :: lnS, theta, c, v3, v4 REAL(r_2) :: dhdS, lambda, crh, int REAL(r_2) :: thetal_max REAL(r_2) :: Sliq REAL(r_2) :: A, B, D, C1 INTEGER(i_d) :: E REAL(r_2) :: F1, F2, F ! REAL(r_2) :: macropore_modifier REAL(r_2) :: cdry, tmp_thetai REAL(r_2), PARAMETER :: tol = 1.e-6_r_2 theta = S*(parin%thre) + (parin%the - parin%thre) var%lambdav = rlambda ! latent heat of vaporisation var%lambdav = 1.91846e6_r_2*((Tsoil+Tzero)/((Tsoil+Tzero)-33.91_r_2))**2 ! Henderson-Sellers, QJRMS, 1984 var%lambdaf = lambdaf ! latent heat of fusion !var%Tfrz = Tfrz(S,parin%he,one/parin%lam) var%Tfrz = Tfrz(S, parin%he, one/(parin%lambc*freezefac)) ! freezefac for test of steep freezing curve IF ((Tsoil < var%Tfrz-tol) .AND. (experiment/=184)) THEN ! ice parin%lam = parin%lambc * freezefac ! freezefac>1 -> steeper freezing curve parin%eta = two/parin%lam + two + one ! freezefac>1 -> steeper freezing curve thetal_max = thetalmax(Tsoil,S,parin%he,one/parin%lam,parin%thre,parin%the) var%dthetaldT = dthetalmaxdT(Tsoil,S,parin%he,one/parin%lam,parin%thre,parin%the) var%iice = 1 var%thetai = MAX((theta - thetal_max),zero) ! volumetric ice content (m3(liq H2O)/m3 soil) tmp_thetai = MAX(MIN(theta, parin%thre) - thetal_max,zero) var%thetal = thetal_max ! liquid water content, relative to saturation ! Sliq = (var%thetal - (parin%the-parin%thre))/parin%thre IF ((parin%thre-tmp_thetai) .LE. MAX(parin%thr,1.e-5_r_2)) THEN Sliq = MAX(parin%thr,1.e-5_r_2) ELSE ! Sliq = min((var%thetal-(parin%the-parin%thre))/(parin%thre-var%thetai), one) Sliq = MIN((var%thetal-(parin%the-parin%thre))/(parin%thre-tmp_thetai), one) ENDIF var%Sliq = Sliq ! saturated liquid water content (< 1 for frozen soil) ! air entry potential, flux matric potential and hydraulic conductivity for saturated frozen soil var%he = parin%he lnS = LOG(Sliq) v3 = EXP(-lnS/parin%lam) v4 = EXP(parin%eta*lnS) var%phie = parin%phie*v3*v4 var%Ksat = parin%Ke*v4 var%h = parin%he*v3 ! matric potential dhdS = zero var%K = var%Ksat var%KS = zero ! var%KT = var%dthetaldT * parin%Ke * parin%eta * exp(lnS*(parin%eta-one))/parin%thre var%KT = var%dthetaldT * parin%Ke * parin%eta * EXP(lnS*(parin%eta-one))/MAX(parin%thre-var%thetai,MAX(parin%thr,1e-5_r_2)) IF (S < one) var%phi = var%phie ! var%phiT = parin%phie * exp(lnS*(parin%eta-one/parin%lam-one)) * var%dthetaldT * & ! (parin%eta-one/parin%lam)/(parin%thre) ! var%phiT = parin%phie * exp(lnS*(parin%eta-one/parin%lam-one)) * var%dthetaldT * & ! (parin%eta-one/parin%lam)/(parin%thre-var%thetai) ! var%phiT = -(((-1 +parin%eta*parin%lam)*parin%phie*(theta - parin%thre)* & ! (var%thetal/(-theta + var%thetal + parin%thre))**(-1 + parin%eta - 1./parin%lam))/ & ! (parin%lam*(-theta + var%thetal + parin%thre)**2))* var%dthetaldT IF (var%isat==0) THEN var%phiT = -( ( (-one +parin%eta*parin%lam) * parin%phie * (MIN(theta,parin%thre) - parin%thre) & * EXP( (-one + parin%eta - one/parin%lam) & * LOG(var%thetal/(-MIN(theta,parin%thre) + var%thetal + parin%thre)) ) ) & / (parin%lam*(-MIN(theta,parin%thre) + var%thetal + parin%thre)**2) ) * var%dthetaldT ! if (S < one) then ! var%phiS = ((parin%eta - 1./parin%lam)*parin%phie* & ! (Sliq)**(1 + parin%eta - 1/parin%lam))/var%thetal*parin%thre ! var%phiS =((-1. + parin%eta*parin%lam)*parin%phie* & ! (var%thetal/(-theta + var%thetal + parin%thre))**(1. + parin%eta - 1./parin%lam))/ & ! (parin%lam*var%thetal)*parin%thre var%phiS =((-one + parin%eta*parin%lam)*parin%phie* & EXP((one + parin%eta - one/parin%lam)*LOG(var%thetal/(-MIN(theta,parin%thre) + var%thetal + parin%thre)))) / & (parin%lam*var%thetal)*parin%thre ELSE var%phiS = zero var%phiT = zero ENDIF !if ( ((parin%thre-var%thetai)<=max(parin%thr,1e-5_r_2)) .or. (Sliq==one) ) then ! var%phiS = zero ! var%phiT = zero !endif var%rh = MAX(EXP(Mw*gravity*var%h/Rgas/(Tsoil+Tzero)),rhmin) ELSE ! no ice parin%lam = parin%lambc ! freezefac>1 -> steeper freezing curve parin%eta = two/parin%lam + two + one ! freezefac>1 -> steeper freezing curve var%he = parin%he var%phie = parin%phie var%Ksat = parin%Ke var%iice = 0 var%thetai = zero var%thetal = S*(parin%thre) + (parin%the - parin%thre) var%dthetaldT = zero dhdS = zero Sliq = S ! liquid water content, relative to saturation lnS = LOG(Sliq) v3 = EXP(-lnS/parin%lam) v4 = EXP(parin%eta*lnS) ENDIF IF (var%thetal < 1.e-12_r_2) THEN ! completely frozen var%dthetaldT = zero var%lambdav = lambdas ! latent heat of sublimation var%KT = zero ENDIF IF (((Tsoil >= var%Tfrz-tol) .OR. (experiment==184)) .AND. (S < one)) THEN ! unsaturated var%h = parin%he*v3 ! matric potential ! dhdS = -parin%he/parin%lam*S**(-one/parin%lam-one) dhdS = -parin%he/parin%lam*EXP((-one/parin%lam-one)*LOG(S)) var%K = parin%Ke*v4 var%phi = parin%phie*v3*v4 var%KS = parin%eta*var%K/S var%KT = zero var%phiS = (parin%eta-one/parin%lam)*var%phi/S var%phiT = zero var%rh = MAX(EXP(Mw*gravity*var%h/Rgas/(Tsoil+Tzero)),rhmin) ENDIF IF (((Tsoil >= var%Tfrz-tol) .OR. (experiment==184)) .AND. (S >= one)) THEN ! saturated var%h = parin%he dhdS = zero var%K = parin%Ke var%phi = parin%phie var%KS = parin%eta*parin%Ke var%phiS = (parin%eta-one/parin%lam)*parin%phie var%phiT = zero var%rh = one var%KT = zero ENDIF ! variables required for vapour phase transfer theta = MIN(S,one)*(parin%thre) + (parin%the - parin%thre) ! if (z.lt.3.and.theta>parin%thfc) then ! macropore_modifier = exp(-parin%zeta*(z-3)) ! ! if (theta > parin%thfc) then ! ! macropore_modifier = 1000._r_2*exp(-var%h/(-4._r_2)) ! var%macropore_factor = macropore_modifier ! else ! var%macropore_factor = one ! endif var%macropore_factor = one var%sl = slope_esat(Tsoil) * Mw/thousand/Rgas/(Tsoil+Tzero) ! m3 m-3 K-1 c = Mw*gravity/Rgas/(Tsoil+Tzero) lambda = parin%lam IF (S < one) THEN var%rh = MAX(EXP(Mw*gravity*var%h/Rgas/(Tsoil+Tzero)),rhmin) ELSE var%rh = one ENDIF crh = c*var%rh var%hS = dhdS var%rhS = crh*dhdS var%cvsat = esat(Tsoil)*Mw/thousand/Rgas/(Tsoil+Tzero) ! m3 m-3 IF (S < one) THEN var%cv = var%rh*var%cvsat var%cvS = var%rhS *var%cvsat var%cvsatT = slope_esat(Tsoil)*Mw/thousand/Rgas/(Tsoil+Tzero) ! m3 m-3 K-1 ELSE var%cv = zero var%cvS = zero var%cvsatT = zero ENDIF ! Penman (1940): tortuosity*theta ! var%Dv = Dva*parin%tortuosity*(parin%the-theta) * ((Tsoil+Tzero)/Tzero)**1.88_r_2 ! m2 s-1 var%Dv = Dva*parin%tortuosity*(parin%the-theta) * EXP(1.88_r_2*LOG((Tsoil+Tzero)/Tzero)) ! m2 s-1 ! ! Moldrup et al. (2003), Table 2 repacked soil (from Moldrup et al. (2000) ! ! (thetas-theta)**1.5 * (thetas-theta)/thetas ! var%Dv = Dva* (parin%the-theta)**1.5_r_2 * (parin%the-theta)/parin%the * ((Tsoil+Tzero)/Tzero)**1.88_r_2 ! m2 s-1 ! int = (-c*parin%he)**lambda * igamma(one-lambda,-c*var%h) int = EXP(lambda*LOG(-c*parin%he)) * igamma(one-lambda,-c*var%h) var%phiv = Dva*parin%tortuosity * (parin%thre*EXP(c*var%h) -parin%thre*int) ! var%phivS = dhdS * Dva*parin%tortuosity * & ! ( (parin%thre)*c*exp(c*var%h) - parin%thre*c*(var%h/parin%he)**(-parin%lam)*exp(c*var%h) ) var%phivS = dhdS * Dva*parin%tortuosity * & ( (parin%thre)*c*EXP(c*var%h) - parin%thre*c*EXP(-parin%lam*LOG(var%h/parin%he))*EXP(c*var%h) ) var%kv = var%Dv * var%cvsat *c * var%rh SELECT CASE (ithermalcond) CASE (0) SELECT CASE (experiment) CASE (11,17) ! Hansson et al. (2004) - special soil properties ! Hansson et al. (2004) - Eq. 13b ! A=C1, B=C2, C1=C3, D=C4, E=C5 A = 0.55_r_2 B = 0.8_r_2 C1 = 3.07_r_2 D = 0.13_r_2 E = 4 CASE default ! calculate v%kH as in Campbell (1985) p.32 eq. 4.20 A = 0.65_r_2 - 0.78_r_2*parin%rho/thousand + 0.60_r_2*(parin%rho/thousand)**2 ! (4.27) B = 2.8_r_2 * (one-parin%thre) ! *theta ! (4.24) !MC if (parin%clay > 0.05) then IF (parin%clay > 0.001_r_2) THEN ! clay=0.001 -> C1=9.2 C1 = one + 2.6_r_2/SQRT(parin%clay*100._r_2) ! (4.28) ELSE C1 = one ENDIF D = 0.03_r_2 + 0.7_r_2*(one-parin%thre)**2 ! (4.22) E = 4 END SELECT SELECT CASE (experiment) CASE (15,17,182:184) ! ice = liquid thermal conductivity ! calculate v%kH as in Campbell (1985) p.32 eq. 4.20 var%kH = A + B*theta-(A-D)*EXP(-(C1*theta)**E) ! (4.20) CASE default ! different thermal conductivity in ice and water ! Hansson et al. (2004) - Eq. 15 F1 = 13.05_r_2 F2 = 1.06_r_2 IF (Tsoil < var%Tfrz-tol .AND. var%thetai.GT.zero ) THEN ! ice ! F = one + F1*var%thetai**F2 F = one + F1*EXP(F2*LOG(var%thetai)) IF ((C1*(theta+F*var%thetai))**E > 100._r_2) THEN var%kH = A + B*(theta+F*var%thetai) ELSE var%kH = A + B*(theta+F*var%thetai)-(A-D)*EXP(-(C1*(theta+F*var%thetai))**E) ENDIF ELSE var%kH = A + B*(theta)-(A-D)*EXP(-(C1*(theta))**E) ENDIF END SELECT CASE(1) ! van de Griend & O'Neill (1986) cdry = 2.0e6_r_2 * (one-parin%thre) ! Sispat Manual (2000), Eq. 2.21 var%kH = one/(cdry + 4.18e6_r_2*theta) * ( (parin%LambdaS + 2300._r_2*theta - 1890._r_2)/0.654_r_2 )**2 END SELECT var%eta_th = one var%kE = var%Dv*var%rh*var%sl*thousand*var%lambdav*var%eta_th var%kth = var%kE + var%kH ! thermal conductivity of soil (includes contribution from vapour phase) var%csoil = parin%css*parin%rho + rhow*cswat*var%thetal + MERGE(0._r_2, rhow*csice*var%thetai, experiment==183) IF (Tsoil < var%Tfrz) THEN ! increase effective heat capacity due to presence of ice var%csoileff = var%csoil + MERGE(0._r_2, rhow*lambdaf*var%dthetaldT, experiment==183) ELSE var%csoileff = var%csoil ENDIF END SUBROUTINE hyofS !============================================================================= SUBROUTINE isosub(iso, c, p, f, fd) IMPLICIT NONE CHARACTER(LEN=2), INTENT(IN) :: iso REAL(r_2), INTENT(IN) :: c REAL(r_2), DIMENSION(:), INTENT(INOUT) :: p REAL(r_2), INTENT(OUT) :: f, fd ! Subroutine to get adsorbed solute (units/g soil) from concn in soil water ! according to chosen isotherm code ("Fr" for Freundlich, "La" for Langmuir ! and "Ll" for Langmuir-linear). ! Definitions of arguments: ! iso - 2 character code. ! c - concn in soil water. ! p(:) - isotherm parameters. ! f - adsorbed mass/g soil. ! fc - deriv of f wrt c (slope of isotherm curve). REAL(r_2) :: x SELECT CASE (iso) CASE ("Fr") IF (p(3)==zero) THEN ! linearise near zero p(3) = (0.01_r_2*dsmmax/p(1))**(one/p(2)) ! concn at 0.01*dsmmax p(4) = p(1)*p(3)**(p(2)-one) ! slope END IF IF (c < p(3)) THEN fd = p(4) f = fd*c ELSE x = p(1)*EXP((p(2)-one)*LOG(c)) f = x*c fd = p(2)*x END IF CASE ("La") x = one/(one+p(2)*c) f = p(1)*c*x fd = p(1)*(x-p(2)*c*x**2) CASE ("Ll") x = one/(one+p(2)*c) f = p(1)*c*x+p(3)*c fd = p(1)*(x-p(2)*c*x**2)+p(3) CASE default WRITE(*,*) "isosub: illegal isotherm type" STOP 2 END SELECT END SUBROUTINE isosub !============================================================================= SUBROUTINE massman_sparse_1d(aa, aah, bb, bbh, cc, cch, dd, ddh, ee, eeh, ff, ffh, gg, ggh, dy, dT, condition, err) USE cable_def_types_mod, ONLY: r_2, i_d USE sli_numbers, ONLY: zero, one IMPLICIT NONE ! in/out REAL(r_2), DIMENSION(:), INTENT(IN) :: aa, aah, bb, bbh, ee, eeh, ff, ffh REAL(r_2), DIMENSION(:), INTENT(IN) :: cc, cch, dd, ddh, gg, ggh REAL(r_2), DIMENSION(:), INTENT(OUT) :: dy, dT INTEGER(i_d), OPTIONAL, INTENT(IN) :: condition INTEGER(i_d), OPTIONAL, INTENT(OUT) :: err ! 0: no error; >0: error ! local INTEGER(i_d) :: n, n2 REAL(r_2), DIMENSION(SIZE(cc),2,2) :: A, B, C REAL(r_2), DIMENSION(SIZE(cc),2) :: d, x ! for conditioning REAL(r_2), DIMENSION(2*SIZE(cc)) :: lST, cST REAL(r_2), DIMENSION(SIZE(cc)) :: lS, lT, cS, cT REAL(r_2), DIMENSION(2*SIZE(cc)*2*SIZE(cc)) :: allvec REAL(r_2), DIMENSION(2*SIZE(cc),2*SIZE(cc)) :: allmat REAL(r_2) :: eps INTEGER(i_d) :: docond ! 0: no conditioning, 1: columns, 2: lines, 3: both INTEGER(i_d) :: ierr ! error code for generic_thomas ! CHARACTER(LEN=20) :: form1 ! integer :: i, nn ! ! check input sizes IF (.NOT. ALL((/SIZE(aa)+1,SIZE(bb)+1,SIZE(dd),SIZE(ee)+1,SIZE(ff)+1,SIZE(gg)/) == SIZE(cc))) THEN WRITE(*,*) 'massman_sparse_1d error1: unequal humidity coeffs.' STOP 2 END IF IF (.NOT. ALL((/SIZE(aah)+1,SIZE(bbh)+1,SIZE(ddh),SIZE(eeh)+1,SIZE(ffh)+1,SIZE(ggh)/) == SIZE(cch))) THEN WRITE(*,*) 'massman_sparse_1d error2: unequal temperature coeffs.' STOP 2 END IF IF (SIZE(cc) /= SIZE(cch)) THEN WRITE(*,*) 'massman_sparse_1d error3: unequal temperature and humidity coeffs.' STOP 2 END IF n = SIZE(cc) IF (PRESENT(condition)) THEN docond = condition ELSE docond = 0 ENDIF IF (PRESENT(err)) err = 0 ierr = 0 ! ! Overall matrix IF (docond >= 1 .AND. docond <= 3) THEN eps = EPSILON(one) n2 = 2*n*2*n allvec(:) = zero allvec(1:n2:4*n+2) = cc(1:n) allvec(2:n2:4*n+2) = dd(1:n) allvec(3:n2-4*n:4*n+2) = ee(1:n-1) allvec(4:n2-4*n:4*n+2) = ff(1:n-1) allvec(2*n+1:n2:4*n+2) = cch(1:n) allvec(2*n+2:n2:4*n+2) = ddh(1:n) allvec(2*n+3:n2-4*n:4*n+2) = eeh(1:n-1) allvec(2*n+4:n2-4*n:4*n+2) = ffh(1:n-1) allvec(4*n+1:n2:4*n+2) = aa(1:n-1) allvec(4*n+2:n2:4*n+2) = bb(1:n-1) allvec(6*n+1:n2:4*n+2) = aah(1:n-1) allvec(6*n+1:n2:4*n+2) = bbh(1:n-1) allmat = RESHAPE(allvec,shape=(/2*n,2*n/),order=(/2,1/)) ENDIF ! Get conditioning numbers SELECT CASE (docond) CASE (1) cST = MAXVAL(ABS(allmat),1) WHERE (cST < eps) cST = one cST = one / cST cS(1:n) = cST(1:2*n-1:2) cT(1:n) = cST(2:2*n:2) lS(1:n) = one lT(1:n) = one CASE (2) lST = MAXVAL(ABS(allmat),2) WHERE (lST < eps) lST = one lST = one / lST lS(1:n) = lST(1:2*n-1:2) lT(1:n) = lST(2:2*n:2) cS(1:n) = one cT(1:n) = one CASE (3) lST = MAXVAL(ABS(allmat),2) WHERE (lST < eps) lST = one lST = one / lST lS(1:n) = lST(1:2*n-1:2) lT(1:n) = lST(2:2*n:2) allmat = SPREAD(lST,2,2*n)*allmat cST = MAXVAL(ABS(allmat),1) WHERE (cST < eps) cST = one cST = one / cST cS(1:n) = cST(1:2*n-1:2) cT(1:n) = cST(2:2*n:2) CASE default cS(1:n) = one cT(1:n) = one lS(1:n) = one lT(1:n) = one END SELECT ! ! fill matrices of generic thomas algorithm A(1,1:2,1:2) = zero A(2:n,1,1) = aa(1:n-1) * lS(2:n) * cS(1:n-1) A(2:n,1,2) = bb(1:n-1) * lS(2:n) * cT(1:n-1) A(2:n,2,1) = aah(1:n-1) * lT(2:n) * cS(1:n-1) A(2:n,2,2) = bbh(1:n-1) * lT(2:n) * cT(1:n-1) B(1:n,1,1) = cc(1:n) * lS(1:n) * cS(1:n) B(1:n,1,2) = dd(1:n) * lS(1:n) * cT(1:n) B(1:n,2,1) = cch(1:n) * lT(1:n) * cS(1:n) B(1:n,2,2) = ddh(1:n) * lT(1:n) * cT(1:n) C(1:n-1,1,1) = ee(1:n-1) * lS(1:n-1) * cS(2:n) C(1:n-1,1,2) = ff(1:n-1) * lS(1:n-1) * cT(2:n) C(1:n-1,2,1) = eeh(1:n-1) * lT(1:n-1) * cS(2:n) C(1:n-1,2,2) = ffh(1:n-1) * lT(1:n-1) * cT(2:n) C(n,1:2,1:2) = zero d(1:n,1) = gg(1:n) * lS(1:n) d(1:n,2) = ggh(1:n) * lT(1:n) ! ! Call Generic Thomas algorithm CALL generic_thomas(n,A,B,C,d,x,ierr) IF (ierr /= 0) THEN IF (PRESENT(err)) THEN err = 1 RETURN ELSE STOP 1 ENDIF ENDIF dy(1:n) = x(1:n,1) * cS(1:n) dT(1:n) = x(1:n,2) * cT(1:n) ! END SUBROUTINE massman_sparse_1d SUBROUTINE massman_sparse_2d(aa, aah, bb, bbh, cc, cch, dd, ddh, ee, eeh, ff, ffh, gg, ggh, dy, dT, condition, err) USE cable_def_types_mod, ONLY: r_2, i_d USE sli_numbers, ONLY: zero, one IMPLICIT NONE ! in/out REAL(r_2), DIMENSION(:,:), INTENT(IN) :: aa, aah, bb, bbh, ee, eeh, ff, ffh REAL(r_2), DIMENSION(:,:), INTENT(IN) :: cc, cch, dd, ddh, gg, ggh REAL(r_2), DIMENSION(:,:), INTENT(OUT) :: dy, dT INTEGER(i_d), OPTIONAL, INTENT(IN) :: condition INTEGER(i_d), OPTIONAL, INTENT(OUT) :: err ! 0: no error; >0: error ! local INTEGER(i_d) :: n, mp INTEGER(i_d) :: n2 REAL(r_2), DIMENSION(1:SIZE(cc,1),1:SIZE(cc,2),2,2) :: A, B, C REAL(r_2), DIMENSION(1:SIZE(cc,1),1:SIZE(cc,2),2) :: d, x ! for conditioning REAL(r_2), DIMENSION(1:SIZE(cc,1),2*SIZE(cc,2)) :: lST, cST REAL(r_2), DIMENSION(1:SIZE(cc,1),SIZE(cc,2)) :: lS, lT, cS, cT REAL(r_2), DIMENSION(1:SIZE(cc,1),2*SIZE(cc,2)*2*SIZE(cc,2)) :: allvec REAL(r_2), DIMENSION(1:SIZE(cc,1),2*SIZE(cc,2),2*SIZE(cc,2)) :: allmat REAL(r_2) :: eps INTEGER(i_d) :: docond ! 0: no conditioning, 1: columns, 2: lines, 3: both INTEGER(i_d) :: ierr ! error code for generic_thomas ! CHARACTER(LEN=20) :: form1 ! integer :: i, k, nn ! ! check input sizes IF (.NOT. ALL((/SIZE(aa,1)+1,SIZE(bb,1)+1,SIZE(dd,1),SIZE(ee,1)+1,SIZE(ff,1)+1,SIZE(gg,1)/) == SIZE(cc,1))) THEN WRITE(*,*) 'massman_sparse_2d error1: unequal humidity coeffs (1st dim).' STOP 2 END IF IF (.NOT. ALL((/SIZE(aah,1)+1,SIZE(bbh,1)+1,SIZE(ddh,1),SIZE(eeh,1)+1,SIZE(ffh,1)+1,SIZE(ggh,1)/) == SIZE(cch,1))) THEN WRITE(*,*) 'massman_sparse_2d error2: unequal temperature coeffs (1st dim).' STOP 2 END IF IF (SIZE(cc,1) /= SIZE(cch,1)) THEN WRITE(*,*) 'massman_sparse_2d error3: unequal temperature and humidity coeffs (1st dim).' STOP 2 END IF IF (.NOT. ALL((/SIZE(aa,2)+1,SIZE(bb,2)+1,SIZE(dd,2),SIZE(ee,2)+1,SIZE(ff,2)+1,SIZE(gg,2)/) == SIZE(cc,2))) THEN WRITE(*,*) 'massman_sparse_2d error4: unequal humidity coeffs (2nd dim).' STOP 2 END IF IF (.NOT. ALL((/SIZE(aah,2)+1,SIZE(bbh,2)+1,SIZE(ddh,2),SIZE(eeh,2)+1,SIZE(ffh,2)+1,SIZE(ggh,2)/) == SIZE(cch,2))) THEN WRITE(*,*) 'massman_sparse_2d error5: unequal temperature coeffs (2nd dim).' STOP 2 END IF IF (SIZE(cc,2) /= SIZE(cch,2)) THEN WRITE(*,*) 'massman_sparse_2d error6: unequal temperature and humidity coeffs (2nd dim).' STOP 2 END IF mp = SIZE(cc,1) n = SIZE(cc,2) IF (PRESENT(condition)) THEN docond = condition ELSE docond = 0 ENDIF IF (PRESENT(err)) err = 0 ierr = 0 ! ! Overall matrix IF (docond >= 1 .AND. docond <= 3) THEN eps = EPSILON(one) n2 = 2*n*2*n allvec(1:mp,:) = zero allvec(1:mp,1:n2:4*n+2) = cc(1:mp,1:n) allvec(1:mp,2:n2:4*n+2) = dd(1:mp,1:n) allvec(1:mp,3:n2-4*n:4*n+2) = ee(1:mp,1:n-1) allvec(1:mp,4:n2-4*n:4*n+2) = ff(1:mp,1:n-1) allvec(1:mp,2*n+1:n2:4*n+2) = cch(1:mp,1:n) allvec(1:mp,2*n+2:n2:4*n+2) = ddh(1:mp,1:n) allvec(1:mp,2*n+3:n2-4*n:4*n+2) = eeh(1:mp,1:n-1) allvec(1:mp,2*n+4:n2-4*n:4*n+2) = ffh(1:mp,1:n-1) allvec(1:mp,4*n+1:n2:4*n+2) = aa(1:mp,1:n-1) allvec(1:mp,4*n+2:n2:4*n+2) = bb(1:mp,1:n-1) allvec(1:mp,6*n+1:n2:4*n+2) = aah(1:mp,1:n-1) allvec(1:mp,6*n+1:n2:4*n+2) = bbh(1:mp,1:n-1) allmat = RESHAPE(allvec,shape=(/mp,2*n,2*n/),order=(/1,3,2/)) ENDIF ! Get conditioning numbers SELECT CASE (docond) CASE (1) cST = MAXVAL(ABS(allmat),2) WHERE (cST < eps) cST = one cST = one / cST cS(1:mp,1:n) = cST(1:mp,1:2*n-1:2) cT(1:mp,1:n) = cST(1:mp,2:2*n:2) lS(1:mp,1:n) = one lT(1:mp,1:n) = one CASE (2) lST = MAXVAL(ABS(allmat),3) WHERE (lST < eps) lST = one lST = one / lST lS(1:mp,1:n) = lST(1:mp,1:2*n-1:2) lT(1:mp,1:n) = lST(1:mp,2:2*n:2) cS(1:mp,1:n) = one cT(1:mp,1:n) = one CASE (3) lST = MAXVAL(ABS(allmat),3) WHERE (lST < eps) lST = one lST = one / lST lS(1:mp,1:n) = lST(1:mp,1:2*n-1:2) lT(1:mp,1:n) = lST(1:mp,2:2*n:2) allmat = SPREAD(lST,2,2*n)*allmat cST = MAXVAL(ABS(allmat),2) WHERE (cST < eps) cST = one cST = one / cST cS(1:mp,1:n) = cST(1:mp,1:2*n-1:2) cT(1:mp,1:n) = cST(1:mp,2:2*n:2) CASE default cS(1:mp,1:n) = one cT(1:mp,1:n) = one lS(1:mp,1:n) = one lT(1:mp,1:n) = one END SELECT ! ! fill matrices of generic thomas algorithm A(1:mp,1,1:2,1:2) = zero A(1:mp,2:n,1,1) = aa(1:mp,1:n-1) * lS(1:mp,2:n) * cS(1:mp,1:n-1) A(1:mp,2:n,1,2) = bb(1:mp,1:n-1) * lS(1:mp,2:n) * cT(1:mp,1:n-1) A(1:mp,2:n,2,1) = aah(1:mp,1:n-1) * lT(1:mp,2:n) * cS(1:mp,1:n-1) A(1:mp,2:n,2,2) = bbh(1:mp,1:n-1) * lT(1:mp,2:n) * cT(1:mp,1:n-1) B(1:mp,1:n,1,1) = cc(1:mp,1:n) * lS(1:mp,1:n) * cS(1:mp,1:n) B(1:mp,1:n,1,2) = dd(1:mp,1:n) * lS(1:mp,1:n) * cT(1:mp,1:n) B(1:mp,1:n,2,1) = cch(1:mp,1:n) * lT(1:mp,1:n) * cS(1:mp,1:n) B(1:mp,1:n,2,2) = ddh(1:mp,1:n) * lT(1:mp,1:n) * cT(1:mp,1:n) C(1:mp,1:n-1,1,1) = ee(1:mp,1:n-1) * lS(1:mp,1:n-1) * cS(1:mp,2:n) C(1:mp,1:n-1,1,2) = ff(1:mp,1:n-1) * lS(1:mp,1:n-1) * cT(1:mp,2:n) C(1:mp,1:n-1,2,1) = eeh(1:mp,1:n-1) * lT(1:mp,1:n-1) * cS(1:mp,2:n) C(1:mp,1:n-1,2,2) = ffh(1:mp,1:n-1) * lT(1:mp,1:n-1) * cT(1:mp,2:n) C(1:mp,n,1:2,1:2) = zero d(1:mp,1:n,1) = gg(1:mp,1:n) * lS(1:mp,1:n) d(1:mp,1:n,2) = ggh(1:mp,1:n) * lT(1:mp,1:n) ! ! Call Generic Thomas algorithm CALL generic_thomas(mp,n,A,B,C,d,x,ierr) IF (ierr /= 0) THEN IF (PRESENT(err)) THEN err = 1 RETURN ELSE STOP 1 ENDIF ENDIF dy(1:mp,1:n) = x(1:mp,1:n,1) * cS(1:mp,1:n) dT(1:mp,1:n) = x(1:mp,1:n,2) * cT(1:mp,1:n) ! END SUBROUTINE massman_sparse_2d !============================================================================= ELEMENTAL PURE SUBROUTINE litter_props(S, Tsoil, vlit, plit, h0) IMPLICIT NONE REAL(r_2), INTENT(IN) :: S REAL(r_2), INTENT(IN) :: Tsoil TYPE(vars), INTENT(INOUT) :: vlit TYPE(params), INTENT(IN) :: plit REAL(r_2), INTENT(IN) :: h0 ! Get soil water variables from S. ! Definitions of arguments: ! S(1:ms) - degree of saturation ("effective satn") of layers. ! ms - no. of soil layers. ! jt(1:ms) - layer soil type nos. ! var(1:ms) - other water vars of layers. REAL(r_2), PARAMETER :: u = 1.0 ! wind speed at litter surface REAL(r_2) :: theta, c REAL(r_2) :: dhdS, sl REAL(r_2) :: rhoL, copo, c_w ! params for thermal vapour transfer enhancement factor (Campbell 1985) REAL(r_2) :: chi, DT0 sl = slope_esat(Tsoil)*Mw/thousand/Rgas/(Tsoil+Tzero) ! m3 m-3 K-1 vlit%sl = sl theta = S*(plit%thre) + (plit%the - plit%thre) c = Mw*gravity/Rgas/(Tsoil+Tzero) rhoL = plit%rho ! dhdS = -plit%he/plit%lam*S**(-one/plit%lam-one)*(thousand/rhoL*plit%the)**(-one/plit%lam) dhdS = -plit%he/plit%lam*EXP((-one/plit%lam-one)*LOG(S))*EXP(-one/plit%lam*LOG(thousand/rhoL*plit%the)) ! vlit%hS = dhdS ! Mathews (2006), A process-based model of offine fuel moisture, ! International Journal of Wildland Fire 15,155-168 chi = 2.08_r_2+u*2.38_r_2 ! (Eq. 45, Tab. 1) DT0 = Dva*EXP(u*2.6_r_2) ! (Eq. 46, Tab. 1) vlit%Dv = DT0*EXP(-half*chi) ! heat and vapour diffusivity half-way through layer (Eq. 11) !write(*,*) 'litter_props', vlit%Dv IF (S < one) THEN ! vlit%h = plit%he*(thousand/rhoL*S*plit%the)**(-one/plit%lam) vlit%h = plit%he*EXP(-one/plit%lam*LOG(thousand/rhoL*S*plit%the)) vlit%K = zero ELSE vlit%K = vlit%Ke vlit%h = plit%he ENDIF vlit%rh = MAX(EXP(Mw*gravity*vlit%h/Rgas/(Tsoil+Tzero)),rhmin) vlit%cvsat = esat(Tsoil)*Mw/thousand/Rgas/(Tsoil+Tzero) ! m3 m-3 vlit%cv = vlit%cvsat*vlit%rh vlit%cvS = vlit%rhS *vlit%cvsat vlit%KS = zero vlit%phiv = vlit%Dv*vlit%cvsat*vlit%rh vlit%phivS = dhdS*vlit%phiv*c vlit%phi = vlit%phiv vlit%phiS = vlit%phivS vlit%rhS = dhdS*c*vlit%rh ! sensible heat conductivity of wet soil W m-1 K-1 (Matthews 2006) vlit%kH = 0.2_r_2 + 0.14_r_2*theta*thousand/rhoL copo = 1932._r_2 c_w = 4.185e6_r_2 vlit%csoil = copo*rhoL + c_w*theta + c_w*h0 !volumetric heat capacity of wet soil vlit%eta_th = one !enhancement factor for transport of water vapour due to a temperature gradient vlit%kE = vlit%Dv*vlit%rh*sl*thousand*rlambda vlit%kth = vlit%kE + vlit%kH ! thermal conductivity of soil END SUBROUTINE litter_props !============================================================================= ! ELEMENTAL PURE SUBROUTINE potential_evap(Rn, rbh, rbw, Ta, rha, Tsoil, k, dz,lambdav, & Ts, E, H, G, & dEdrha, dEdTs, dEdTsoil, dGdTa, dGdTsoil,iice) ! Pennman-Monteith equation, with additional account for heat flux into the surface IMPLICIT NONE REAL(r_2), INTENT(IN) :: Rn, rbh, rbw, Ta, rha, Tsoil, k, dz,