! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/soil_hyd.F SUBROUTINE SOIL_HYD (ETRAN_L, DT, DRAINAGE, W_FLUX, SUBSOIL_L) IMPLICIT NONE include "size.h" include "mtlm.h" ! Description: ! Increments the layer soil moisture contents and calculates ! calculates gravitational runoff. Calls the following: ! ! HYD_CON - to calculate the hydraulic conductivity ! (Cox, 6/95) ! ! DARCY - to calculate the Darcian fluxes between soil layers ! (Cox, 6/95) ! ! ! Documentation : UM Documentation Paper 25 ! !~Input / output variables~ ! DT ! Timestep (s) ! ETRAN_L ! Extraction of water from each soil layer (kg/m2/s) ! DRAINAGE ! Drainage from the base of the soil profile (kg/m2/s). ! W_FLUX ! The fluxes of water between layers (kg/m2/s) REAL ETRAN_L(POINTS, NGND), DT REAL DRAINAGE(POINTS) REAL W_FLUX(POINTS, 0:NGND) !~Optional ponding parameters~ ! Ponding parameters ! POND_MAX ! Maximum amount of ponded water (kg/m2) ! POND_FRAC_MAX ! Maximum ponded fraction of the gridbox REAL POND_MAX, POND_FRAC_MAX PARAMETER(POND_MAX = 0.0, POND_FRAC_MAX = 0.5) !~Other variables~ ! I,J,N ! Loop counters. ! DM ! The transfer of soil moisture within a layer(kg/m2/timestep). ! DZ_STRUCT ! Effective thickness of soil layers when using structural excess ice (m) ! MMAX ! The maximum moisture content of each layer (kg/m2) ! MSAT ! The saturation moisture content of each layer (kg/m2). ! MU ! Unfrozen soil moisture contents of each layer (kg/m2) ! SUK ! Fractional saturation of lowest layer INTEGER I, L, N REAL DM(POINTS), MMAX(POINTS, NGND) REAL MU(POINTS, NGND), SUK(POINTS), MSAT(POINTS, NGND) REAL DZ_STRUCT(POINTS,NGND), SUBSOIL_L(POINTS, NGND) REAL MF(POINTS, NGND) !----------------------------------------------------------------------- ! Calculate the unfrozen soil moisture contents and the saturation ! total soil moisture for each layer. !----------------------------------------------------------------------- DZ_STRUCT(:,:) = 0. DO I = 1, LAND_PTS L = LAND_INDEX(I) DO N = 1, NSOIL(L) MSAT(L,N)=RHO_W*DZ_GND(N)*V_SAT(L,N) MU(L,N)=SU(L,N)*MSAT(L,N) MF(L,N) = MAX(M(L,N) - MU(L,N), 0.) ENDDO ENDDO MMAX(:,:) = MSAT(:,:) !----------------------------------------------------------------------- ! Allow for some ponding of water by permitting excess moisture to ! remain in the top layer. !----------------------------------------------------------------------- DO I=1,LAND_PTS L=LAND_INDEX(I) MMAX(L,1)=MSAT(L,1)+POND_MAX ENDDO !---------------------------------------------------------------------------- ! Calculate the soil (unfrozen) moisture fluxes in absence of any restriction ! ie. don't worry about layers having too much or too little water. !---------------------------------------------------------------------------- ! Bottom Layer, bottom flux DO I=1,LAND_PTS L=LAND_INDEX(I) !---------------------------------------------------------------------------- ! Free drainage from the bottom of the soil column so W = K !---------------------------------------------------------------------------- SUK(L)=SU(L,NSOIL(L)) CALL HYD_CON(B(L,NSOIL(L)),KS(L,NSOIL(L)),SUK(L) &, W_FLUX(L,NSOIL(L))) !---------------------------------------------------------------------------- ! Calculate the Darcian fluxes between layers !---------------------------------------------------------------------------- IF (NSOIL(L).gt. 1) then DO N=NSOIL(L),2,-1 CALL DARCY(B(L,N-1),KS(L,N-1),PHIS(L,N-1) &, SU(L,N-1), DZ_GND(N-1) &, B(L,N), KS(L,N), PHIS(L,N) &, SU(L,N), DZ_GND(N) &, W_FLUX(L,N-1)) ENDDO ENDIF !----------------------------------------------------------------------- ! Now restrict the fluxes. Make sure that: ! 1) The flux is not too large so as to saturate whatever layer it goes into ! 2) Enough moisture is present in the source layer to satisfy the flux !----------------------------------------------------------------------- ! Bottom layer, bottom flux DM(L)=(W_FLUX(L,NSOIL(L))+ETRAN_L(L,NSOIL(L)))*DT IF (DM(L).GT.(0.25*MU(L,NSOIL(L)))) THEN DM(L)=0.25*MU(L,NSOIL(L)) W_FLUX(L,NSOIL(L))=DM(L)/DT-ETRAN_L(L,NSOIL(L)) ENDIF M(L,NSOIL(L))= M(L,NSOIL(L))-DT*(W_FLUX(L,NSOIL(L)) & + ETRAN_L(L,NSOIL(L))) MU(L,NSOIL(L)) = MU(L,NSOIL(L))-DT*(W_FLUX(L,NSOIL(L)) & + ETRAN_L(L,NSOIL(L))) ! Next deal with the middle layers IF (NSOIL(L) .gt. 1) then DO N = NSOIL(L)-1, 1, -1 IF (W_FLUX(L,N) .gt. 0.) then ! Flux is directed downwards. Check layer N+1 to ensure no supersaturation IF (W_FLUX(L,N) .gt. (MMAX(L,N+1) - M(L,N+1))/DT) then W_FLUX(L,N) = (MMAX(L,N+1) - M(L,N+1))/DT ENDIF ! Check layer N to ensure adequate supply of liquid moisture. DM(L)=(W_FLUX(L,N)+ETRAN_L(L,N))*DT IF (DM(L).GT.0.25*MU(L,N)) THEN DM(L)=0.25*MU(L,N) W_FLUX(L,N)=DM(L)/DT-ETRAN_L(L,N) ENDIF ! Update layer below with the moisture flux now. M(L,N+1) = M(L,N+1) + DT*W_FLUX(L,N) MU(L,N+1) = MU(L,N+1) + DT*W_FLUX(L,N) ! Remove mositure from layer N M(L,N) = M(L,N) - DT*(W_FLUX(L,N) + ETRAN_L(L,N)) MU(L,N) = MU(L,N) - DT*(W_FLUX(L,N) + ETRAN_L(L,N)) ELSE ! Flux is directed upwards. Check layer N to ensure no supersaturation IF ((-W_FLUX(L,N)) .gt. (MMAX(L,N) - M(L,N))/DT) then W_FLUX(L,N) = -(MMAX(L,N) - M(L,N))/DT ENDIF ! Check layer N+1 to ensure adequate supply of liquid moisture. DM(L)=(-W_FLUX(L,N))*DT IF (DM(L).GT. 0.25*MU(L,N+1)) THEN DM(L)=0.25*MU(L,N+1) W_FLUX(L,N)=-DM(L)/DT ENDIF ! Remove moisture from layer N+1 (W_FLUX < 0). M(L,N+1) = M(L,N+1) + DT*W_FLUX(L,N) MU(L,N+1) = MU(L,N+1) + DT*W_FLUX(L,N) ! Add moisture to layer N M(L,N) = M(L,N) - DT*(W_FLUX(L,N)+ETRAN_L(L,N)) MU(L,N) = MU(L,N) - DT*(W_FLUX(L,N) + ETRAN_L(L,N)) ENDIF ENDDO ENDIF ! Lastly, deal with the water flux from the surface (layer 0) IF (W_FLUX(L,0) .gt. (MMAX(L,1) - M(L,1))/DT) then SURF_ROFF(L) = SURF_ROFF(L) & + W_FLUX(L,0) - (MMAX(L,1)-M(L,1))/DT W_FLUX(L,0) = (MMAX(L,1) - M(L,1))/DT ENDIF M(L,1) = M(L,1) + DT*(W_FLUX(L,0)) MU(L,1) = MU(L,1) + DT*(W_FLUX(L,0)) SUBSOIL_L(L,:) = 0. SUBSOIL_L(L,1) = SUBSOIL(L) if (SUBSOIL(L) .gt. 0.) then if (SUBSOIL(L) .lt. 0.25*MF(L,1)/DT) then MF(L,1) = MF(L,1) - SUBSOIL(L)*DT M(L,1) = M(L,1) - SUBSOIL(L)*DT else do N = 1, NSOIL(L) if (SUBSOIL_L(L,N) .gt. 0.25*MF(L,N)/DT) then SUBSOIL_L(L,N+1) = SUBSOIL_L(L,N+1) + SUBSOIL_L(L,N) & - 0.25*MF(L,N)/DT SUBSOIL_L(L,N) = 0.25*MF(L,N)/DT MF(L,N) = MF(L,N) - SUBSOIL_L(L,N)*DT M(L,N) = M(L,N) - SUBSOIL_L(L,N)*DT endif enddo endif endif !---------------------------------------------------------------------- ! Update frozen and unfrozen moisture contents & drainage !---------------------------------------------------------------------- DO N = 1, NSOIL(L) SF(L,N) = MAX(MF(L,N)/MSAT(L,N), 0.) SU(L,N) = MU(L,N) / MSAT(L,N) ENDDO ENDDO DRAINAGE(:) = 0.0 DO I=1,LAND_PTS L=LAND_INDEX(I) DRAINAGE(L)=W_FLUX(L,NSOIL(L)) ENDDO RETURN END !----------------------------------------------------------------------- SUBROUTINE DARCY (B1, KS1, PHIS1, SU1, DZ1 &, B2, KS2, PHIS2, SU2, DZ2 &, W_FLUX) IMPLICIT NONE ! ! Description: ! Calculates the Darcian fluxes between adjacent soil layers. ! (Cox, 6/95) ! ! Documentation : UM Documentation Paper 25 ! Arguments with intent(IN) : ! B1 / B2 ! Clapp-Hornberger exponent. ! DZ1 / DZ2 ! Layer thickness (m). ! DZSTRUCT1 / DZSTRUCT2 ! Thickness of structural excess ice (m) ! DZEFF1 / DZEFF2 ! Effective thickness of soil layers (m) ! K ! Actual hydraulic conductivity ! KS ! Saturated hydraulic conductivity (kg/m2/s) ! PHIS1 / PHIS2 ! Saturated soil water pressure (m). ! PSI1 / PSI2 ! Actural soil water pressure (suction) (m) ! SU1 / SU2 ! Unfrozen soil moisture content as a fraction of saturation. ! WFLUX ! The flux of water between layers (kg/m2/s). REAL B1, DZ1, DZ2, KS1, PHIS1, SU1, SU2 REAL B2, KS2, PHIS2, KSK, BK REAL DZSTRUCT1, DZSTRUCT2 REAL PSI1, PSI2, K1, K2, K REAL DZEFF1, DZEFF2, W_FLUX DZEFF1 = DZ1 DZEFF2 = DZ2 !----------------------------------------------------------------------- ! Calculate the soil water suction of the layers !----------------------------------------------------------------------- IF (SU1.LE.0.01) THEN ! Prevent blow up for dry soil. PSI1=PHIS1/(0.01**B1) ELSEIF (SU1.GT.0.01.AND.SU1.LE.1.0) THEN PSI1=PHIS1/(SU1**B1) ELSE PSI1=PHIS1 ENDIF IF (SU2.LE.0.01) THEN ! Prevent blow up for dry soil. PSI2=PHIS2/(0.01**B2) ELSEIF (SU2.GT.0.01.AND.SU2.LE.1.0) THEN PSI2=PHIS2/(SU2**B2) ELSE PSI2=PHIS2 ENDIF CALL HYD_CON(B1,KS1,SU1 &, K1) CALL HYD_CON(B2,KS2,SU2 &, K2) !----------------------------------------------------------------------- ! Determine hydraulic conductivity and soil moisture fluxes !----------------------------------------------------------------------- K = (DZEFF2*K1+DZEFF1*K2)/(DZEFF2+DZEFF1) W_FLUX=K*(2.0*(PSI2-PSI1)/(DZEFF2+DZEFF1)+1) RETURN END !----------------------------------------------------------------------- SUBROUTINE HYD_CON (B,KS,SU &, K) IMPLICIT NONE ! ! Description: ! Calculates the hydraulic conductivity (Cox, 6/95) ! ! ! Documentation : UM Documentation Paper 25 ! !~Variables~ ! B ! Exponent in conductivity and soil water suction fits ! DZ ! Soil layer thickness (m) ! DZ_STRUCT ! Thickness of structural excess ice (m) ! K ! The hydraulic conductivity (kg/m2/s). ! KICE ! Hydraulic conductivity of ice (kg/m2) ! KS ! The saturated hydraulic conductivity (kg/m2) ! SU ! Fractional saturation of unfrozen moisture REAL B, KS, SU, KICE, K PARAMETER (KICE = 0.) IF (SU.GE.0.0.AND.SU.LT.1.0) THEN K=KS*SU**(2*B+3) ELSEIF (SU.LT.0.0) THEN K=0.0 ELSE K=KS ENDIF ! If structural excess ice, determine hydraulic conductivity as the volume weighted ! average of the conductivity of soil and of ice RETURN END