From fba0e4a30e02a374679e39b0478a35e012dccf55 Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh Date: Sun, 7 Apr 2024 21:14:31 +1000 Subject: [PATCH 1/8] Replace cable_gw_hydro.F90 with its new version --- src/offline/cable_define_types.F90 | 18 +- src/offline/cable_parameters.F90 | 574 +++++- src/science/gw_hydro/cable_gw_hydro.F90 | 2333 ++++++++++++++--------- src/science/soilsnow/cbl_thermal.F90 | 12 +- src/util/cable_common.F90 | 67 +- 5 files changed, 2006 insertions(+), 998 deletions(-) diff --git a/src/offline/cable_define_types.F90 b/src/offline/cable_define_types.F90 index 46f8977d91..9d01874758 100644 --- a/src/offline/cable_define_types.F90 +++ b/src/offline/cable_define_types.F90 @@ -162,6 +162,8 @@ MODULE cable_def_types_mod rhosoil_vec,& !soil density [kg/m3] ssat_vec, & !volumetric water content at saturation [mm3/mm3] watr, & !residual water content of the soil [mm3/mm3] + smpc_vec, & ! Hutson Cass SWC potential cutoff ! 2 lines inserted by rk4417 - phase2 + wbc_vec, & ! Hutson Cass SWC volumetric water cutoff sfc_vec, & !field capcacity (hk = 1 mm/day) swilt_vec ! wilting point (hk = 0.02 mm/day) @@ -187,6 +189,8 @@ MODULE cable_def_types_mod GWssat_vec, & !saturated water content of the aquifer [mm3/mm3] GWwatr, & !residual water content of the aquifer [mm3/mm3] GWz, & !node depth of the aquifer [m] + smpc_GW, & ! Hutson Cass SWC potential cutoff ! 2 lines inserted by rk4417 - phase2 + wbc_GW, & ! Hutson Cass SWC volumetric water cutoff GWdz, & !thickness of the aquifer [m] GWrhosoil_vec !density of the aquifer substrate [kg/m3] @@ -240,7 +244,7 @@ MODULE cable_def_types_mod owetfac, & ! surface wetness fact. at previous time step t_snwlr, & ! top snow layer depth in 3 layer snowpack tggav, & ! mean soil temperature in K - otgg, & ! soil temperature in K +! otgg, & ! soil temperature in K ! moved below by rk4417 - phase2 otss, & ! surface temperature (weighted soil, snow) otss_0, & ! surface temperature (weighted soil, snow) tprecip, & @@ -267,6 +271,7 @@ MODULE cable_def_types_mod sdepth, & ! snow depth smass, & ! snow mass ssdn, & ! snow densities + otgg, & ! soil temperature in K ! moved here from above by rk4417 - phase2 tgg, & ! soil temperature in K tggsn, & ! snow temperature in K dtmlt, & ! water flux to the soil @@ -880,6 +885,10 @@ SUBROUTINE alloc_soil_parameter_type(var, mp) ALLOCATE( var%ssat_vec(mp,ms) ) ALLOCATE( var%watr(mp,ms) ) var%watr(:,:) = 0.05 + ALLOCATE( var%wbc_GW(mp) ) ! block inserted by rk4417 - phase2 + ALLOCATE( var%smpc_GW(mp) ) + ALLOCATE( var%wbc_vec(mp,ms) ) + ALLOCATE( var%smpc_vec(mp,ms) ) ! end - rk4417 - phase ALLOCATE( var%sfc_vec(mp,ms) ) ALLOCATE( var%swilt_vec(mp,ms) ) ALLOCATE( var%sand_vec(mp,ms) ) @@ -977,7 +986,8 @@ SUBROUTINE alloc_soil_snow_type(var, mp) ALLOCATE( var%t_snwlr(mp) ) ALLOCATE( var%wbfice(mp,ms) ) ALLOCATE( var%tggav(mp) ) - ALLOCATE( var%otgg(mp) ) +! ALLOCATE( var%otgg(mp) ) ! replaced by below - rk4417 - phase2 + ALLOCATE( var%otgg(mp,ms) ) ALLOCATE( var%otss(mp) ) ALLOCATE( var%otss_0(mp) ) ALLOCATE( var%tprecip(mp) ) @@ -1520,6 +1530,10 @@ SUBROUTINE dealloc_soil_parameter_type(var) DEALLOCATE( var% cnsd_vec ) DEALLOCATE( var%hyds_vec ) DEALLOCATE( var%sucs_vec ) + DEALLOCATE( var%wbc_GW ) ! block inserted by rk4417 - phase2 + DEALLOCATE( var%smpc_GW ) + DEALLOCATE( var%wbc_vec ) + DEALLOCATE( var%smpc_vec ) ! end - rk4417 - phase DEALLOCATE( var%bch_vec ) DEALLOCATE( var%ssat_vec ) DEALLOCATE( var%watr ) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 534f73d43e..e69ead88c9 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -1757,6 +1757,20 @@ SUBROUTINE write_cnp_params(veg, casaflux, casamet) END SUBROUTINE write_cnp_params !============================================================================ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) + + ! ______ MMY: derived_parameters is changed by MMY @ 13 Oct 2022 ______ + ! MMY required input sand_vec, silt_vec, clay_vec ,rhosoil_vec and org_vec + ! MMY units: soil%hyds_vec : mm/s + ! soil%sand_vec, soil%clay_vec, soil%silt_vec : fraction + ! soil%sucs_vec : mm, positive + ! soil%bch_vec : - + ! soil%watr, soil%swilt_vec, soil%sfc_vec, soil%ssat_vec: m3/m3 + ! soil%org_vec : m3/m3 + ! soil%rhosoil_vec : kg/m3 + ! soil%css: "J/kg/K" + ! soil%cnsd: "W/m/K" + ! soil%wbc_vec, soil%smpc_vec + ! Gives values to parameters that are derived from other parameters. TYPE (soil_snow_type), INTENT(INOUT) :: ssnow TYPE (veg_parameter_type), INTENT(IN) :: veg @@ -1765,122 +1779,466 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) TYPE (balances_type), INTENT(INOUT) :: bal TYPE (roughness_type), INTENT(INOUT) :: rough - INTEGER :: j,i,klev ! do loop counter +! INTEGER :: j,i,klev ! do loop counter + INTEGER :: j,i,klev,k ! do loop counter ! replaces line above - rk4417 - phase2 REAL(r_2) :: temp(mp) REAL :: tmp2(mp) REAL(r_2), DIMENSION(mp,ms) :: perc_frac - REAL(r_2), DIMENSION(17) :: psi_o,psi_c - REAL(r_2), DIMENSION(mp,ms) :: psi_tmp - REAL(r_2), DIMENSION(ms) :: soil_depth + REAL(r_2), DIMENSION(17) :: psi_o,psi_c ! psi_o is not needed - rk4417 - phase2 - soil_depth(1) = REAL(soil%zse(1),r_2) - DO klev=2,ms - soil_depth(klev) = soil_depth(klev-1) + REAL(soil%zse(klev),r_2) - END DO +! REAL(r_2), DIMENSION(mp,ms) :: psi_tmp +! REAL(r_2), DIMENSION(ms) :: soil_depth + +! 2 lines above replaced by below - rk4417 - phase2 + + REAL(r_2), DIMENSION(mp,ms) :: psi_tmp, sst_tmp ! MMY add sst_tmp + REAL(r_2), DIMENSION(mp,ms) :: soil_depth ! MMY,rhosoil_temp + + REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: ssat_bounded,rho_soil_bulk ! added by rk4417 - phase2 + +! line below inserted by rk4417 - phase2 + where(veg%iveg .eq. 17) soil%isoilm = 9 ! MMY@13April it says where iveg = ice, isoilm = permanent ice + +! soil_depth(1) = REAL(soil%zse(1),r_2) +! DO klev=2,ms +! soil_depth(klev) = soil_depth(klev-1) + REAL(soil%zse(klev),r_2) +! END DO + +! block above replaced by block below - rk4417 - phase2 - psi_o(1:3) = -66000._r_2 - psi_o(4) = -35000._r_2 - psi_o(5) = -83000._r_2 - psi_o(6:17) = -74000._r_2 - psi_c(1:3) = -2550000._r_2 + soil_depth(:,1) = soil%zse_vec(:,1) + do klev=2,ms + soil_depth(:,klev) = soil_depth(:,klev-1) + soil%zse_vec(:,klev) + end do + +! psi_o(1:3) = -66000._r_2 ! not needed - rk4417 - phase2 +! psi_o(4) = -35000._r_2 +! psi_o(5) = -83000._r_2 +! psi_o(6:17) = -74000._r_2 + + psi_c(1:3) = -2550000._r_2 ! used to be in the old cable_common.F90 - rk4417 - phase2 psi_c(4) = -2240000._r_2 psi_c(5) = -4280000._r_2 psi_c(6:17) = -2750000._r_2 + + ! _______ MMY@19May2023 move back since moving below will cause gw_off issue ________ ! Construct derived parameters and zero initialisations, ! regardless of where parameters and other initialisations ! have loaded from: soil%zshh(1) = 0.5 * soil%zse(1) ! distance between consecutive layer - ! midpoints: + ! midpoints: soil%zshh(ms + 1) = 0.5 * soil%zse(ms) soil%zshh(2:ms) = 0.5 * (soil%zse(1:ms-1) + soil%zse(2:ms)) !MD aquifer node depth soil%GWz = 0.5*soil%GWdz + SUM(soil%zse) !node is halfway through aquifer depth +! MMY note: GWz doesn't function, kept for later + +! IF (cable_user%GW_MODEL) THEN +! +! DO klev=1,ms +! soil%hyds_vec(:,klev) = 0.0070556*10.0**(-0.884 + 0.0153*soil%Sand_Vec(:,klev)*100.0)* & +! EXP(-gw_params%hkrz*(MAX(0.,soil_depth(klev)-gw_params%zdepth))) +! soil%sucs_vec(:,klev) = 10.0 * 10.0**(1.88 -0.0131*soil%Sand_Vec(:,klev)*100.0) +! soil%bch_vec(:,klev) = 2.91 + 0.159*soil%Clay_Vec(:,klev)*100.0 +! soil%ssat_vec(:,klev) = 0.489 - 0.00126*soil%Sand_Vec(:,klev)*100.0 +! soil%watr(:,klev) = 0.02 + 0.00018*soil%Clay_Vec(:,klev)*100.0 +! ENDDO +! !aquifer share non-organic with last layer if not found in param file +! IF (found_explicit_gw_parameters .EQV. .FALSE.) THEN +! soil%GWhyds_vec(:) = soil%hyds_vec(:,ms) +! soil%GWsucs_vec(:) = soil%sucs_vec(:,ms) +! soil%GWbch_vec(:) = soil%bch_vec(:,ms) +! soil%GWssat_vec(:) = soil%ssat_vec(:,ms) +! soil%GWwatr(:) = soil%watr(:,ms) +! ENDIF +! !include organin impact. fraction of grid cell where percolation through +! !organic macropores dominates +! soil%Org_Vec = MAX(0._r_2,soil%Org_Vec) +! soil%Org_Vec = MIN(1._r_2,soil%Org_Vec) +! DO klev=1,3 !0-23.3 cm, data really is to 30cm +! soil%hyds_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%hyds_vec(:,klev) + & +! soil%Org_Vec(:,klev)*gw_params%org%hyds_vec_organic +! soil%sucs_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%sucs_vec(:,klev) + & +! soil%Org_Vec(:,klev)*gw_params%org%sucs_vec_organic +! soil%bch_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%bch_vec(:,klev) +& +! soil%Org_Vec(:,klev)*gw_params%org%clappb_organic +! soil%ssat_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%ssat_vec(:,klev) + & +! soil%Org_Vec(:,klev)*gw_params%org%ssat_vec_organic +! soil%watr(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%watr(:,klev) + & +! soil%Org_Vec(:,klev)*gw_params%org%watr_organic +! END DO +! +! !vegetation dependent field capacity (point plants get stressed) and +! !wilting point +! DO i=1,mp +! psi_tmp(i,:) = -psi_c(veg%iveg(i)) +! END DO +! soil%sfc_vec = (soil%ssat_vec-soil%watr) * (ABS(psi_tmp)/(ABS(soil%sucs_vec)))**(-1.0/soil%bch_vec)+& +! soil%watr +! DO i=1,mp +! psi_tmp(i,:) = -psi_c(veg%iveg(i)) +! END DO +! soil%swilt_vec = (soil%ssat_vec-soil%watr) * (ABS(psi_tmp)/(ABS(soil%sucs_vec)))**(-1.0/soil%bch_vec)+& +! soil%watr +! +! !set the non-vectored values to srf value +! soil%sfc(:) = REAL(soil%sfc_vec(:,1)) +! soil%swilt(:) = REAL(soil%swilt_vec(:,1)) +! +! !convert the units back to what default uses and GW only uses the +! !vectored versions +! soil%hyds = REAL(soil%hyds_vec(:,1))/1000.0 +! soil%sucs = REAL(soil%sucs_vec(:,1))/1000.0 +! soil%ssat = REAL(soil%ssat_vec(:,1)) +! soil%bch = REAL(soil%bch_vec(:,1)) +! +! DO i=1,mp +! soil%slope(i) = MIN(0.9,MAX(1e-9,soil%slope(i))) +! soil%slope_std(i) = MIN(0.9,MAX(1e-9,soil%slope_std(i))) +! END DO +! +! IF ((gw_params%MaxSatFraction .LT. -9999.9) .AND. (mp .EQ. 1)) soil%slope(:) = 0.01 +! +! ELSE +! +! soil%sfc_vec = REAL(SPREAD(soil%sfc(:),2,ms),r_2) +! soil%swilt_vec = REAL(SPREAD(soil%swilt(:),2,ms),r_2) +! !These are not used when gw_model == false +! soil%watr = 0._r_2 +! soil%GWwatr = 0._r_2 +! +! END IF +! +! +! IF ( .NOT. soilparmnew) THEN ! Q,Zhang @ 12/20/2010 +! soil%cnsd = soil%sand * 0.3 + soil%clay * 0.25 & +! + soil%silt * 0.265 ! set dry soil thermal conductivity +! ! [W/m/K] +! END IF +! +! soil%hsbh = soil%hyds*ABS(soil%sucs) * soil%bch ! difsat*etasat +! soil%ibp2 = NINT(soil%bch) + 2 +! ! Ticket #66 +! WHERE( soil%ssat > 0.) & ! Avoid divide by +! soil%pwb_min = (soil%swilt/soil%ssat)**soil%ibp2 +! soil%i2bp3 = 2 * NINT(soil%bch) + 3 +! rough%hruff = MAX(0.01, veg%hc - 1.2 * ssnow%snowd/MAX(ssnow%ssdnn, 100.)) +! rough%hruff_grmx = rough%hruff +! ! owetfac introduced by EAK apr2009 +! ssnow%owetfac = MAX(0.0, MIN(1.0, & +! (REAL(ssnow%wb(:, 1)) - soil%swilt) / & +! (soil%sfc - soil%swilt))) +! temp(:) = 0.0 +! tmp2(:) = 0.0 +! WHERE ( ssnow%wbice(:, 1) > 0. ) ! Prevents divide by zero at glaciated +! ! points where wb and wbice=0. +! temp(:) = ssnow%wbice(:, 1) / ssnow%wb(:, 1) +! tmp2(:) = REAL(temp(:)) +! ssnow%owetfac = ssnow%owetfac * (1.0 - tmp2(:)) ** 2 +! +! END WHERE +! ssnow%pudsto = 0.0 +! ssnow%pudsmx = 0.0 +! +! ! Initialise sum flux variables: +! sum_flux%sumpn = 0.0 +! sum_flux%sumrp = 0.0 +! sum_flux%sumrpw = 0.0 +! sum_flux%sumrpr = 0.0 +! sum_flux%sumrs = 0.0 +! sum_flux%sumrd = 0.0 +! sum_flux%dsumpn = 0.0 +! sum_flux%dsumrp = 0.0 +! sum_flux%dsumrd = 0.0 +! ! Initialise conservation variables: +! bal%precip_tot = 0.0 +! bal%rnoff_tot = 0.0 +! bal%evap_tot = 0.0 +! bal%wbal_tot = 0.0 +! bal%ebal_tot = 0.0 +! bal%ebal_tot_cncheck = 0.0 +! bal%drybal = 0.0 +! bal%wetbal = 0.0 +! bal%wbtot0 = 0.0 +! bal%RadbalSum = 0.0 +! DO j=1, ms +! bal%wbtot0 = bal%wbtot0 + REAL(ssnow%wb(:, j)) * soil%zse(j) & +! * 1000.0 +! END DO +! bal%osnowd0 = ssnow%osnowd +! +! ! vh_js ! comment out hide% condition +! ! IF (hide%Ticket49Bug6) THEN +! +! IF(cable_user%SOIL_STRUC=='sli') THEN +! ! Only 1 horizon by default ! +! soil%nhorizons = 1 +! soil%ishorizon = 1 +! END IF +! ! END IF + +! above commented-out portion replaced by remainder of subroutine below - rk4417 - phase2 + + IF (cable_user%GW_MODEL) then + soil%qhz_max(:) = real(gw_params%MaxHorzDrainRate,r_2) !enable distributed values + soil%hkrz(:) = real(gw_params%hkrz,r_2) + soil%zdepth(:) = real(gw_params%zdepth,r_2) - IF (cable_user%GW_MODEL) THEN + !placeholder + soil%srf_frac_ma(:) = 0._r_2 + soil%edepth_ma(:) = 0._r_2 + + soil%qhz_efold(:) = real(gw_params%EfoldHorzDrainScale,r_2)*soil%drain_dens(:)& + + real(gw_params%EfoldHorzDrainRate,r_2) - DO klev=1,ms - soil%hyds_vec(:,klev) = 0.0070556*10.0**(-0.884 + 0.0153*soil%Sand_Vec(:,klev)*100.0)* & - EXP(-gw_params%hkrz*(MAX(0.,soil_depth(klev)-gw_params%zdepth))) - soil%sucs_vec(:,klev) = 10.0 * 10.0**(1.88 -0.0131*soil%Sand_Vec(:,klev)*100.0) - soil%bch_vec(:,klev) = 2.91 + 0.159*soil%Clay_Vec(:,klev)*100.0 - soil%ssat_vec(:,klev) = 0.489 - 0.00126*soil%Sand_Vec(:,klev)*100.0 - soil%watr(:,klev) = 0.02 + 0.00018*soil%Clay_Vec(:,klev)*100.0 - ENDDO - !aquifer share non-organic with last layer if not found in param file - IF (found_explicit_gw_parameters .EQV. .FALSE.) THEN - soil%GWhyds_vec(:) = soil%hyds_vec(:,ms) - soil%GWsucs_vec(:) = soil%sucs_vec(:,ms) - soil%GWbch_vec(:) = soil%bch_vec(:,ms) - soil%GWssat_vec(:) = soil%ssat_vec(:,ms) - soil%GWwatr(:) = soil%watr(:,ms) - ENDIF !include organin impact. fraction of grid cell where percolation through !organic macropores dominates - soil%Org_Vec = MAX(0._r_2,soil%Org_Vec) - soil%Org_Vec = MIN(1._r_2,soil%Org_Vec) - DO klev=1,3 !0-23.3 cm, data really is to 30cm - soil%hyds_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%hyds_vec(:,klev) + & - soil%Org_Vec(:,klev)*gw_params%org%hyds_vec_organic - soil%sucs_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%sucs_vec(:,klev) + & - soil%Org_Vec(:,klev)*gw_params%org%sucs_vec_organic - soil%bch_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%bch_vec(:,klev) +& - soil%Org_Vec(:,klev)*gw_params%org%clappb_organic - soil%ssat_vec(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%ssat_vec(:,klev) + & - soil%Org_Vec(:,klev)*gw_params%org%ssat_vec_organic - soil%watr(:,klev) = (1.-soil%Org_Vec(:,klev))*soil%watr(:,klev) + & - soil%Org_Vec(:,klev)*gw_params%org%watr_organic - END DO + soil%org_vec = max(0._r_2,soil%org_vec) + soil%org_vec = min(1._r_2,soil%org_vec) - !vegetation dependent field capacity (point plants get stressed) and - !wilting point - DO i=1,mp - psi_tmp(i,:) = -psi_c(veg%iveg(i)) - END DO - soil%sfc_vec = (soil%ssat_vec-soil%watr) * (ABS(psi_tmp)/(ABS(soil%sucs_vec)))**(-1.0/soil%bch_vec)+& - soil%watr - DO i=1,mp - psi_tmp(i,:) = -psi_c(veg%iveg(i)) + DO klev=1,ms + do i=1,mp + if (abs(soil%sand_vec(i,klev) + soil%clay_vec(i,klev) +& + soil%silt_vec(i,klev)-1.0) .gt. 0.1) then + soil%sand_vec(i,klev) = 0.4 + soil%clay_vec(i,klev) = 0.2 + soil%silt_vec(i,klev) = 0.4 + endif + END DO END DO - soil%swilt_vec = (soil%ssat_vec-soil%watr) * (ABS(psi_tmp)/(ABS(soil%sucs_vec)))**(-1.0/soil%bch_vec)+& - soil%watr - !set the non-vectored values to srf value - soil%sfc(:) = REAL(soil%sfc_vec(:,1)) - soil%swilt(:) = REAL(soil%swilt_vec(:,1)) + DO klev=1,ms;do i=1,mp + ! MMY use single var (uni) or two var (multi) regression from Cosby 1984 + ! MMY or use Hutson-Cass SWC + if (gw_params%cosby_univariate) then + soil%hyds_vec(i,klev) = 0.0070556*10.0**(-0.884 + 1.53*soil%sand_vec(i,klev))* & + exp(-soil%hkrz(i)*(MAX(0.,soil_depth(i,klev)-soil%zdepth(i)))) ! MMY add MAX(0.,...) + ! MMY add MAX(0.,...) for stopping hyds increase in shallow soil + ! MMY NEED to figure out the source of this equation + soil%sucs_vec(i,klev) = 10.0 * 10.0**(1.88 -1.31*soil%sand_vec(i,klev)) + soil%bch_vec(i,klev) = 2.91 + 15.9*soil%clay_vec(i,klev) + soil%ssat_vec(i,klev) = min(0.489,max(0.1, 0.489 - 0.126*soil%sand_vec(i,klev) ) ) + !forgot source but not from cosby and not for BC characteristic function + soil%watr(i,klev) = 0.02 + 0.018*soil%clay_vec(i,klev) !forgot + soil%wbc_vec(i,klev) = 0.0 + soil%smpc_vec(i,klev) = 0.0 + elseif (gw_params%cosby_multivariate) then + soil%hyds_vec(i,klev) = 0.00706*(10.0**(-0.60 + 1.26*soil%sand_vec(i,klev) + & + -0.64*soil%clay_vec(i,klev) ) )*& + exp(-soil%hkrz(i)*(MAX(0.,soil_depth(i,klev)-soil%zdepth(i)))) ! MMY add MAX(0.,...) + soil%sucs_vec(i,klev) = 10.0 * 10.0**(1.54 - 0.95*soil%sand_vec(i,klev) + & + 0.63*soil%silt_vec(i,klev) ) + soil%bch_vec(i,klev) = 3.1 + 15.4*soil%clay_vec(i,klev) - & + 0.3*soil%sand_vec(i,klev) + soil%ssat_vec(i,klev) = 0.505 - 0.142*soil%sand_vec(i,klev) - & + 0.037*soil%clay_vec(i,klev) + !forgot source but not from cosby and not for BC characteristic function + soil%watr(i,klev) = 0.02 + 0.018*soil%clay_vec(i,klev) + soil%wbc_vec(i,klev) = 0.0 + soil%smpc_vec(i,klev) = 0.0 + + elseif (gw_params%HC_SWC) THEN + ! Hutson-Cass SWC : seperate dry/wet avoid discont in derv at smp=sucs + ! pedotransfer from T. Mayra, N.J. Jarvis, Geoderma 91, 1999, + ! https://doi.org/10.1016/S0016-7061(98)00129-3 + soil%sucs_vec(i,klev) = 10.0 * 10.0** ( -4.98403 +& + 5.0922*soil%sand_vec(i,klev) + & + 15.751*soil%silt_vec(i,klev) + & + 0.124090*soil%rhosoil_vec(i,klev) - & + 16.4000*soil%org_vec(i,klev) - & + 21.76*(soil%silt_vec(i,klev)**2.0) + & + 14.382*(soil%silt_vec(i,klev)**3.0) + & + 8.0407*(soil%clay_vec(i,klev)**2.0) + & + 44.06*(soil%org_vec(i,klev)**2.0) ) + + soil%bch_vec(i,klev) = 10.0**(1.0 / (-0.84669 - & + 0.4680*soil%sand_vec(i,klev) & + +0.9246*soil%silt_vec(i,klev) & + -0.4543*soil%rhosoil_vec(i,klev) & + -0.04979*soil%org_vec(i,klev) & + +3.2947*(soil%sand_vec(i,klev)**2.0) & + -1.689*(soil%sand_vec(i,klev)**3.0) & + +11.2*(soil%org_vec(i,klev)**3.0) )) + + soil%ssat_vec(i,klev) = 0.234597 + & + 0.466142*soil%sand_vec(i,klev) + & + 0.88163*soil%silt_vec(i,klev) + & + 0.643386*soil%clay_vec(i,klev) - & + 0.3028160*soil%rhosoil_vec(i,klev) + & + 0.179762*(soil%sand_vec(i,klev)**2.0) - & + 0.03134631*(soil%silt_vec(i,klev)**2.0) + soil%hyds_vec(i,klev) = 0.00706*(10.0**(-0.60 & + +1.26*soil%sand_vec(i,klev) & + -0.64*soil%clay_vec(i,klev)))* & + exp(-soil%hkrz(i)* & + (soil_depth(i,klev)-soil%zdepth(i))) + soil%watr(i,klev) = 0.0 + if (klev .eq. 1) soil%GWwatr(i) = 0.0 + else + print *, "MMY Please choose the pedotransfer function from cosby_univariate, cosby_multivariate and HC_SWC" + end if + end do; end do + + ! MMY take soil organic material impact into consideration + if (.not.gw_params%HC_SWC) then + DO klev=1,ms !0-23.3 cm, data really is to 30cm + do i=1,mp + soil%hyds_vec(i,klev) = (1.-soil%org_vec(i,klev))*soil%hyds_vec(i,klev) + & + soil%org_vec(i,klev)*gw_params%org%hyds_organic* & + exp(-soil%hkrz(i)*(soil_depth(i,klev)-soil%zdepth(i))) + soil%sucs_vec(i,klev) = (1.-soil%org_vec(i,klev))*soil%sucs_vec(i,klev) + & + soil%org_vec(i,klev)*gw_params%org%sucs_organic + soil%bch_vec(i,klev) = (1.-soil%org_vec(i,klev))*soil%bch_vec(i,klev) +& + soil%org_vec(i,klev)*gw_params%org%bch_organic + soil%ssat_vec(i,klev) = (1.-soil%org_vec(i,klev))*soil%ssat_vec(i,klev) + & + soil%org_vec(i,klev)*gw_params%org%ssat_organic + soil%watr(i,klev) = (1.-soil%org_vec(i,klev))*soil%watr(i,klev) + & + soil%org_vec(i,klev)*gw_params%org%watr_organic + END DO + END DO + end if + + ! MMY sfc_vec (field capacity ): "The volumetric water content at field capacity is derived by assuming a + ! hydraulic conductivity of 0.1 mm day-1 (0.1 mm day-1 / 86400s = 1.157407e-06 mm/s = gw_params%sfc_vec_hk) + ! and inverting the hydraulic conductivity function" from CLM 4.5 tech note Page 85 + ! MMY swilt_vec (wilting point): vegetation dependent wilting point (point plants get stressed), inverting the water retention curve + + do klev=1,ms + do i=1,mp + if (soil%isoilm(i) .ne. 9 .and. veg%iveg(i) .le. 16) then + + psi_tmp(i,klev) = abs(psi_c(veg%iveg(i))) + + soil%swilt_vec(i,klev) = (soil%ssat_vec(i,klev)-soil%watr(i,klev)) * & + (psi_tmp(i,klev)/soil%sucs_vec(i,klev))& + **(-1.0/soil%bch_vec(i,klev))+soil%watr(i,klev) + soil%sfc_vec(i,klev) = (soil%ssat_vec(i,klev)-soil%watr(i,klev)) * & + (gw_params%sfc_vec_hk/soil%hyds_vec(i,klev)) & + **(1.0/(2.0*soil%bch_vec(i,klev)+3.0)) + & + soil%watr(i,klev) + soil%swilt_vec(i,klev) = min(0.95*soil%sfc_vec(i,klev),soil%swilt_vec(i,klev)) + + else + + soil%swilt_vec(i,klev) = soil%swilt(i) + soil%sfc_vec(i,klev) = soil%sfc(i) + + end if + end do + end do + ! MMY GW does not use non-vectored values so they are simply set as the 3rd soil layer's value rather + ! than derived here. Choosing 3rd layer since the top layer is thin and lack of representation. + soil%bch = real(soil%bch_vec(:,3)) + soil%sand = real(soil%sand_vec(:,3)) + soil%clay = real(soil%clay_vec(:,3)) + soil%silt = real(soil%silt_vec(:,3)) + soil%sfc = real(soil%sfc_vec(:,3)) + soil%swilt = real(soil%swilt_vec(:,3)) + soil%ssat = real(soil%ssat_vec(:,3)) + soil%rhosoil = real(soil%rhosoil_vec(:,3)) + !convert the units back to what default uses and GW only uses the !vectored versions - soil%hyds = REAL(soil%hyds_vec(:,1))/1000.0 - soil%sucs = REAL(soil%sucs_vec(:,1))/1000.0 - soil%ssat = REAL(soil%ssat_vec(:,1)) - soil%bch = REAL(soil%bch_vec(:,1)) - - DO i=1,mp - soil%slope(i) = MIN(0.9,MAX(1e-9,soil%slope(i))) - soil%slope_std(i) = MIN(0.9,MAX(1e-9,soil%slope_std(i))) - END DO + soil%hyds = real(soil%hyds_vec(:,3))/1000.0 ! MMY mm/s -> m/s + soil%sucs = -1.* ABS(real(soil%sucs_vec(:,3))/1000.0) ! MMY mm -> m + + ! MMY assume aquifer share soil params with last layer + soil%GWhyds_vec(:) = soil%hyds_vec(:,ms) + soil%GWsucs_vec(:) = soil%sucs_vec(:,ms) + soil%GWbch_vec(:) = soil%bch_vec(:,ms) + soil%GWssat_vec(:) = soil%ssat_vec(:,ms) + soil%GWwatr(:) = soil%watr(:,ms) + + ! MMY adjust slope and slope_std + do i=1,mp + soil%slope(i) = min(0.9,max(1e-5,soil%slope(i))) + soil%slope_std(i) = min(0.9,max(1e-5,soil%slope_std(i))) + end do - IF ((gw_params%MaxSatFraction .LT. -9999.9) .AND. (mp .EQ. 1)) soil%slope(:) = 0.01 + if ((gw_params%MaxSatFraction .lt. -9999.9) .and. (mp .eq. 1)) soil%slope(:) = 0.01 - ELSE + ELSE !not gw model - soil%sfc_vec = REAL(SPREAD(soil%sfc(:),2,ms),r_2) - soil%swilt_vec = REAL(SPREAD(soil%swilt(:),2,ms),r_2) - !These are not used when gw_model == false - soil%watr = 0._r_2 - soil%GWwatr = 0._r_2 + !These are not used when gw_model == false + soil%watr = 0._r_2 + soil%GWwatr = 0._r_2 - END IF + END IF + IF (cable_user%soil_thermal_fix) then + IF ((.not.(gw_params%cosby_univariate.or.gw_params%cosby_multivariate)) .or. .not.cable_user%gw_model) THEN + WRITE(logn,*) 'OVER WRITING CSS_VEC, RHOSOIL_VEC, CNSD_VEC read from gw_type input' + WRITE(logn,*) 'Forcing values consistant with soil_thermal_fix = true ' + WRITE(*,*) 'OVER WRITING CSS_VEC, RHOSOIL_VEC, CNSD_VEC read from gw_type input' + WRITE(*,*) 'Forcing values consistant with soil_thermal_fix = true ' + END IF - IF ( .NOT. soilparmnew) THEN ! Q,Zhang @ 12/20/2010 - soil%cnsd = soil%sand * 0.3 + soil%clay * 0.25 & - + soil%silt * 0.265 ! set dry soil thermal conductivity - ! [W/m/K] - END IF + if (allocated(ssat_bounded)) deallocate(ssat_bounded) + if (allocated(rho_soil_bulk)) deallocate(rho_soil_bulk) + + allocate(ssat_bounded(size(soil%ssat_vec,dim=1), size(soil%ssat_vec,dim=2) ) ) + ssat_bounded(:,:) = min( 0.8, max(0.1, soil%ssat_vec(:,:) ) ) + + allocate(rho_soil_bulk(size(soil%rhosoil_vec,dim=1), size(soil%rhosoil_vec,dim=2) ) ) + rho_soil_bulk(:,:) = min(2500.0, max(500.0 , (2700.0*(1.0 - ssat_bounded(:,:)) ) ) ) + + do klev=1,ms + do i=1,mp + soil%cnsd_vec(i,klev) = ( (0.135*(1.0-ssat_bounded(i,klev))) +& + (64.7/soil%rhosoil_vec(i,klev)) ) / & + (1.0 - 0.947*(1.0-ssat_bounded(i,klev))) + ! ??????????? MMY ??????????? + soil%rhosoil_vec(i,klev) = soil%rhosoil_vec(i,klev)/(1.0-soil%ssat_vec(i,klev)) + ! ??????????????????????????? + !took avg of results from A New Perspective on Soil Thermal Properties Ochsner, Horton,Tucheng + !Soil Sci Soc America 2001 + !to find what silt (1.0-sand-clay) is !simply regress to his means !in J/kg/K + soil%css_vec(i,klev) = max(910.6479*soil%silt_vec(i,klev) +& + 916.4438 * soil%clay_vec(i,klev) +& + 740.7491*soil%sand_vec(i,klev), 800.0) + ! ??????????????? MMY which eq are right ??????????????? + ! sst_tmp(i,klev) = soil%ssat_vec(i,klev) + ! sst_tmp(i,klev) = 1.0- MAX(0.15, MIN(0.85, sst_tmp(i,klev))) + ! soil%css_vec(i,klev) = (1.0-soil%org_vec(i,klev)) * & + ! (850*(1.0 - soil%sand_vec(i,klev) - soil%clay_vec(i,klev)) & + ! +865.0*soil%clay_vec(i,klev) + 750.0*soil%sand_vec(i,klev))& + ! + soil%org_vec(i,klev)*950.0 + ! soil%cnsd_vec(i,klev) = (1.0-soil%org_vec(i,klev)) * & + ! (0.135*sst_tmp(i,klev) + 0.0239/sst_tmp(i,klev) )/ & + ! (1.0 - 0.947*sst_tmp(i,klev)) & + ! + soil%org_vec(i,klev)*0.05 + ! ???????????????????????????????????????????????????????? + end do + end do + + if (allocated(ssat_bounded)) deallocate(ssat_bounded) + if (allocated(rho_soil_bulk)) deallocate(rho_soil_bulk) + + IF (cable_user%gw_model) then !organic correction? + do klev=1,ms + do i=1,mp + soil%css_vec(i,klev) = (1.-soil%org_vec(i,klev))*soil%css_vec(i,klev) + & + real(soil%org_vec(i,klev)*gw_params%org%css_organic) + soil%cnsd_vec(i,klev) = (1.-soil%org_vec(i,klev))*soil%cnsd_vec(i,klev) + & + soil%org_vec(i,klev)*gw_params%org%cnsd_organic + end do + end do + END IF + ! MMY to unify with hydraulic parameters - using 3rd layer's value + soil%rhosoil = real(soil%rhosoil_vec(:,3)) + soil%css = real(soil%css_vec(:,3)) + soil%cnsd = real(soil%cnsd_vec(:,3)) + ! _____ MMY @Oct2022 comment out as soilparmnew is default now _____ + ! ELSEIF ( .NOT. soilparmnew) THEN ! Q,Zhang @ 12/20/2010 + ! soil%cnsd = soil%sand * 0.3 + soil%clay * 0.25 + soil%silt * 0.265 + ! !set dry soil thermal conductivity [W/m/K] + ! soil%cnsd_vec = spread(soil%cnsd,2,ms) + ! ________________________________________________________________ + END IF soil%hsbh = soil%hyds*ABS(soil%sucs) * soil%bch ! difsat*etasat soil%ibp2 = NINT(soil%bch) + 2 @@ -1891,9 +2249,7 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) rough%hruff = MAX(0.01, veg%hc - 1.2 * ssnow%snowd/MAX(ssnow%ssdnn, 100.)) rough%hruff_grmx = rough%hruff ! owetfac introduced by EAK apr2009 - ssnow%owetfac = MAX(0.0, MIN(1.0, & - (REAL(ssnow%wb(:, 1)) - soil%swilt) / & - (soil%sfc - soil%swilt))) + ssnow%owetfac = MAX(0.0, MIN(1.0, (REAL(ssnow%wb(:, 1)) - soil%swilt) / (soil%sfc - soil%swilt))) temp(:) = 0.0 tmp2(:) = 0.0 WHERE ( ssnow%wbice(:, 1) > 0. ) ! Prevents divide by zero at glaciated @@ -1933,7 +2289,7 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) END DO bal%osnowd0 = ssnow%osnowd - ! vh_js ! comment out hide% condition + !! vh_js !! comment out hide% condition ! IF (hide%Ticket49Bug6) THEN IF(cable_user%SOIL_STRUC=='sli') THEN @@ -1943,6 +2299,54 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) END IF ! END IF + where(soil%ssat_vec .gt. 0.0) + ssnow%wblf = max(0.01_r_2,ssnow%wbliq/soil%ssat_vec) + ssnow%wbfice = ssnow%wbice / soil%ssat_vec + elsewhere + ssnow%wblf =0.01 + ssnow%wbfice =0.99 + endwhere + + ! ____________________ MMY comment out as we don't use hys ___________________________ + do k=1,ms + do i=1,mp + if (ssnow%wb_hys(i,k) .lt. 0._r_2) then + ssnow%wb_hys(i,k) = ssnow%wb(i,k) + end if + ssnow%wb_hys(i,k) = max(soil%watr(i,k) ,min(soil%ssat_vec(i,k), ssnow%wb_hys(i,k))) + + if (ssnow%smp_hys(i,k) .lt. -1.0e+30_r_2) then !set to missing, calc + ssnow%smp_hys(i,k) = -soil%sucs_vec(i,k)* & + ( (ssnow%wb_hys(i,k)-ssnow%watr_hys(i,k))/& + (ssnow%ssat_hys(i,k)-ssnow%watr_hys(i,k)) )**& + (-1._r_2/soil%bch_vec(i,k) ) + end if + ssnow%smp_hys(i,k) = max(-1.0e10,min(-soil%sucs_vec(i,k),ssnow%smp_hys(i,k) )) + end do + end do + + ! MMY note that after commenting out, wb_hys and smp_hys are not defined any more + + ! if (cable_user%gw_model .and. gw_params%bc_hysteresis) then + ! do klev=1,ms + ! do i=1,mp + ! if (soil%isoilm(i) .ne. 9 .and. veg%iveg(i) .le. 16) then + + ! psi_tmp(i,klev) = abs(psi_c(veg%iveg(i))) + + ! soil%swilt_vec(i,klev) = (ssnow%ssat_hys(i,klev)-ssnow%watr_hys(i,klev)) * & + ! (psi_tmp(i,klev)/soil%sucs_vec(i,klev))& + ! **(-1.0/soil%bch_vec(i,klev))+& + ! ssnow%watr_hys(i,klev) + ! soil%sfc_vec(i,klev) = (gw_params%sfc_vec_hk/soil%hyds_vec(i,klev))& + ! **(1.0/(2.0*soil%bch_vec(i,klev)+3.0)) *& + ! (ssnow%ssat_hys(i,klev)-ssnow%watr_hys(i,klev)) + ssnow%watr_hys(i,klev) + ! end if + ! end do + ! end do + + ! end if + END SUBROUTINE derived_parameters !============================================================================ SUBROUTINE check_parameter_values(soil, veg, ssnow) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index e7cf339188..ead9055936 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -1,6 +1,6 @@ !============================================================================== ! This source code is part of the -! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model. +! Austrian Community Atmosphere Biosphere Land Exchange (CABLE) model. ! This work is licensed under the CABLE Academic User Licence Agreement ! (the "Licence"). ! You may not use this file except in compliance with the Licence. @@ -41,27 +41,37 @@ MODULE cable_gw_hydro_module USE cable_common_module, ONLY : gw_params,cable_user,& cable_runtime,& - max_glacier_snowd - -!distribute these per sbr -USE cable_phys_constants_mod, ONLY : CTFRZ => TFRZ -USE cable_phys_constants_mod, ONLY : CHL => HL -USE cable_phys_constants_mod, ONLY : CHLF => HLF -USE cable_phys_constants_mod, ONLY : CHLS => HLS -USE cable_phys_constants_mod, ONLY : Cdensity_liq => density_liq -USE cable_phys_constants_mod, ONLY : Cdensity_ice => density_ice -USE cable_phys_constants_mod, ONLY : Ccgsnow => cgsnow -USE cable_phys_constants_mod, ONLY : Ccswat => cswat -USE cable_phys_constants_mod, ONLY : Ccs_rho_wat => cs_rho_wat -USE cable_phys_constants_mod, ONLY : Ccs_rho_ice => cs_rho_ice -USE cable_math_constants_mod, ONLY : CPI => PI - + max_glacier_snowd,& + psi_c ! psi_c added by rk4417 - phase2 + +! USE cbl_soil_snow_subrs_module, ONLY : trimb, snow_processes_soil_thermal ! replaced below by rk4417 - phase2 + +! line below commented out by rk4417 - phase2 +! USE cable_data_module, only: C=>PHYS ! all constants used in this module belong to PHYS + + +!distribute these per sbr ! added by rk4417 - phase2 + USE cable_phys_constants_mod, ONLY : CTFRZ => TFRZ + USE cable_phys_constants_mod, ONLY : CHL => HL + USE cable_phys_constants_mod, ONLY : CHLF => HLF + USE cable_phys_constants_mod, ONLY : CHLS => HLS ! not needed - rk4417 - phase2 + USE cable_phys_constants_mod, ONLY : Cdensity_liq => density_liq + USE cable_phys_constants_mod, ONLY : Cdensity_ice => density_ice + USE cable_phys_constants_mod, ONLY : Ccgsnow => cgsnow + USE cable_phys_constants_mod, ONLY : Ccswat => cswat + USE cable_phys_constants_mod, ONLY : Ccs_rho_wat => cs_rho_wat ! not needed - rk4417 - phase2 + USE cable_phys_constants_mod, ONLY : Ccs_rho_ice => cs_rho_ice ! not needed - rk4417 - phase2 + USE cable_phys_constants_mod, ONLY : Ccsice => csice ! inserted by rk4417 - phase2 + USE cable_math_constants_mod, ONLY : CPI => PI + IMPLICIT NONE PRIVATE !mrd561 GW params - REAL(r_2), SAVE :: smp_cor = 8.0 + +! REAL(r_2), SAVE :: smp_cor = 8.0 ! dependends on FEEDBACK - rk4417 - phase2 + REAL(r_2), PARAMETER :: sucmin = -1.0e8 ! minimum soil pressure head [mm] REAL(r_2), PARAMETER :: volwatmin = 1e-4 !min soil water [mm] REAL(r_2), PARAMETER :: wtd_uncert = 0.1 ! uncertaintiy in wtd calcultations [mm] @@ -71,176 +81,187 @@ MODULE cable_gw_hydro_module REAL(r_2), PARAMETER :: Sy_deep = 0.1 REAL(r_2), PARAMETER :: dz_deep = 50000.0 + REAL(r_2), SAVE :: den_rat=0.921 + INTEGER, PARAMETER :: wtd_iter_max = 20 ! maximum number of iterations to find the water table depth - ! ! This module contains the following subroutines that + REAL(r_2), PARAMETER :: m2mm = 1000.0 + REAL(r_2), PARAMETER :: mm2m = 0.001 + REAL(r_2), PARAMETER :: zero=0.0 + + ! This module contains the following subroutines that !are called from other modules - PUBLIC :: soil_snow_gw,calc_srf_wet_fraction,sli_hydrology,& - pore_space_relative_humidity,set_unsed_gw_vars -CONTAINS + PUBLIC :: soil_snow_gw, calc_srf_wet_fraction, sli_hydrology,& + pore_space_relative_humidity, set_unsed_gw_vars, den_rat + !,set_den_rat ! MMY@Nov2022 comment out set_den_rat +CONTAINS ! ----------------------------------------------------------------------------- - SUBROUTINE GWsoilfreeze(dels, soil, ssnow,tgg_old) - !NOTE: this is only included because gw_model uses parameters XXX_vec - !these are r_2. this breaks bitwise compatibility with trunk - !if acceptable this routine does the same thing but with r_2 soil params - ! if max_ice_frac always set to frozen_limit and tgg_tmp is always CTFRZ -IMPLICIT NONE + SUBROUTINE GWsoilfreeze(dels, soil, ssnow) + + !*## Purpose + ! Snow freezes and melts. + ! this is only included because gw_model uses parameters XXX_vec + ! these are r_2. this breaks bitwise compatibility with trunk + ! if acceptable this routine does the same thing but with r_2 soil params + ! if max_ice_frac always set to frozen_limit and tgg_tmp is always CTFRZ + ! + !## Method + ! + ! The subroutine takes a float as input, prints its value then adds 2/3, + ! prints it and returns it. + ! + ! The equation used is: + ! \[ myarg = myarg + !frac{2.}{3.} \] + ! + !## References + ! + ! The guidelines for documentation can be found in: + ! [CABLE developer guide](https://cable-lsm.github.io/CABLE/developer_guide/doc_guide/science_doc/) REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow TYPE(soil_parameter_type), INTENT(INOUT) :: soil - REAL, INTENT(INOUT), DIMENSION(mp,ms) :: tgg_old - REAL(r_2), DIMENSION(mp) :: sicefreeze - REAL(r_2), DIMENSION(mp) :: sicemelt - REAL(r_2), DIMENSION(mp,ms) :: wbice_delta,avail_por + + REAL , DIMENSION(mp,ms) :: tgg_tmp !calc wb ice at approx tgg + REAL(r_2), DIMENSION(mp,ms) :: wbice_delta,avail_por,delta_ice_vol REAL(r_2), DIMENSION(mp) :: ice_mass,liq_mass,tot_mass + INTEGER :: i,j,k - REAL, DIMENSION(mp,ms) :: tgg_tmp - REAL(r_2),DIMENSION(mp,ms) :: xx,max_ice_frac,iceF,den_css !Decker and Zeng 2009 + REAL(r_2) :: Aconst,Dconst + REAL, DIMENSION(mp,ms) :: gammzz_snow + REAL(r_2),DIMENSION(mp,ms) :: xx,max_ice_frac,den_css !Decker and Zeng 2009 ! MMY iceF is useless - max_ice_frac(:,:) = 0.0 - DO k=1,ms - DO i=1,mp - IF (ssnow%tgg(i,k) .LT. CTFRZ .AND. soil%ssat_vec(i,k) .GT. 1.0e-8) THEN - max_ice_frac(i,k) = (1. - EXP(-2.*(ssnow%wb(i,k)/soil%ssat_vec(i,k))**4.0 *& - (ssnow%tgg(i,k)-CTFRZ)))/EXP(1. - ssnow%wb(i,k)/soil%ssat_vec(i,k)) - max_ice_frac(i,k) = MAX(0.4,max_ice_frac(i,k))*ssnow%wb(i,k) + max_ice_frac(:,:) = 0.0 + delta_ice_vol(:,:) = 0.0 + tgg_tmp(:,:) = ssnow%tgg(:,:) + gammzz_snow(:,:) = 0._r_2 - wbice_delta(i,k) = MAX(0.,max_ice_frac(i,k) - ssnow%wbice(i,k)) + k=1 + do i=1,mp + if (ssnow%isflag(i) .eq. 0 .and. soil%isoilm(i) .ne. 9) then + gammzz_snow(i,k) = real(Ccgsnow,r_2) * real(ssnow%snowd(i),r_2) + end if + end do - avail_por(i,k) = soil%ssat_vec(i,k) - ssnow%wbliq(i,k)+& - Cdensity_ice/Cdensity_liq*wbice_delta(i,k) - & - ssnow%wbice(i,k) + do k=1,ms + do i=1,mp - wbice_delta(i,k) = MIN(wbice_delta(i,k),avail_por(i,k)) + if ( (ssnow%tgg(i,k) .lt. CTFRZ) .and. & + (ssnow%tgg(i,k) .lt. ssnow%otgg(i,k)) ) then ! MMY freezing - max_ice_frac(i,k) = ssnow%wbice(i,k) + wbice_delta(i,k) + ssnow%otgg(i,k) = min(ssnow%otgg(i,k),CTFRZ) - ENDIF - END DO - END DO + tgg_tmp(i,k) = ssnow%tgg(i,k) + Aconst = -2.0*( (0.2+ssnow%wb(i,k)/soil%ssat_vec(i,k))**4.0 ) + Dconst = exp(1. - min(1.0,0.2+ssnow%wb(i,k)/soil%ssat_vec(i,k))) - tgg_tmp(:,:) = tgg_old(:,:) - DO k=1,ms - DO i=1,mp - IF (soil%isoilm(i) .EQ. 9) THEN - tgg_tmp(i,k) = CTFRZ - ELSE - IF (ssnow%tgg(i,k) .LE. CTFRZ) THEN - IF (tgg_old(i,k) .GT. CTFRZ) THEN - tgg_tmp(i,k) = CTFRZ - END IF - ELSE - IF (tgg_old(i,k) .LE. CTFRZ) THEN - tgg_tmp(i,k) = CTFRZ - END IF - END IF - END IF - END DO - END DO + max_ice_frac(i,k) = (1._r_2 - exp(2._r_2*(min(1.,ssnow%wb(i,k)/soil%ssat_vec(i,k))**4.0) *& + real(tgg_tmp(i,k)-CTFRZ,r_2)))/exp(1._r_2- min(1.0,ssnow%wb(i,k)/soil%ssat_vec(i,k))) + if (soil%isoilm(i) .eq. 9) max_ice_frac(i,k) = 0.85_r_2 + ! MMY maximum ice can get under current wb + max_ice_frac(i,k) = min(0.9_r_2,max_ice_frac(i,k)) - !allow more freezing for permenant glacier ice regions - DO i=1,mp - IF (soil%isoilm(i) .EQ. 9) max_ice_frac(i,:) = 0.85_r_2*ssnow%wb(i,:) - END DO + !delta_ice_vol(i,k) = max(0._r_2, ssnow%wb(i,k)*max_ice_frac(i,k) - ssnow%wbice(i,k)) ! MMY + ! MMY ssnow%wb(i,k)*max_ice_frac(i,k) here should be ice volume, thus divided by den_rat + delta_ice_vol(i,k) = max(0._r_2, ssnow%wb(i,k)*max_ice_frac(i,k)/den_rat - ssnow%wbice(i,k)) ! MMY - DO k = 1, ms - DO i=1,mp + ! MMY wbliq can be lower than watr due to condensation, but to avoid large negative water potential + ! set watr as wbliq lower boundary + !check amount of water we have + delta_ice_vol(i,k) = min((ssnow%wbliq(i,k)-soil%watr(i,k))/den_rat, & + max(0._r_2, delta_ice_vol(i,k) ) ) - ice_mass(i) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*REAL(Cdensity_ice,r_2) - liq_mass(i) = ssnow%wbliq(i,k)*soil%zse_vec(i,k)*REAL(Cdensity_liq,r_2) - tot_mass(i) = liq_mass(i) + ice_mass(i) - - IF (ssnow%tgg(i,k) .LE. CTFRZ .AND. & - ssnow%tgg(i,k) .LT. tgg_tmp(i,k) .AND. & - max_ice_frac(i,k) - ssnow%wbice(i,k) > .001) THEN - - sicefreeze(i) = MIN( MAX( 0.0_r_2, ( max_ice_frac(i,k) - & - ssnow%wbice(i,k) ) ) * soil%zse_vec(i,k) * Cdensity_ice, & - ( tgg_tmp(i,k) - ssnow%tgg(i,k) ) * ssnow%gammzz(i,k) / Chlf ) - - ssnow%wbice(i,k) = MIN( ssnow%wbice(i,k) +& - sicefreeze(i)/soil%zse_vec(i,k)/Cdensity_ice,& - max_ice_frac(i,k) ) - ssnow%gammzz(i,k) = MAX(soil%heat_cap_lower_limit(i,k), & - (1.0-soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k) & - + (ssnow%wb(i,k)-ssnow%wbice(i,k)) * REAL(Ccs_rho_wat,r_2) & - + ssnow%wbice(i,k) * REAL(Ccs_rho_ice,r_2)& - )*soil%zse_vec(i,k) - - IF (k .EQ. 1 .AND. ssnow%isflag(i) .EQ. 0) THEN - ssnow%gammzz(i,k) = ssnow%gammzz(i,k) + Ccgsnow * ssnow%snowd(i) - END IF + delta_ice_vol(i,k) = min(delta_ice_vol(i,k),max(0._r_2,& + real(ssnow%otgg(i,k)-ssnow%tgg(i,k),r_2)*ssnow%gammzz(i,k)/& + (soil%zse_vec(i,k)*real(CHLF*Cdensity_ice,r_2)) ) ) - ssnow%tgg(i,k) = ssnow%tgg(i,k) + REAL(sicefreeze(i)) & - * Chlf / REAL(ssnow%gammzz(i,k) ) + elseif ((ssnow%tgg(i,k) .gt. CTFRZ) .and. & + (ssnow%tgg(i,k) .gt. ssnow%otgg(i,k)) .and. ssnow%wbice(i,k) .gt. 0.0) then ! MMY melting - ELSEIF ( ssnow%tgg(i,k) .GT. tgg_tmp(i,k) .AND. & - ssnow%wbice(i,k) .GT. max_ice_frac(i,k) ) THEN + ssnow%otgg(i,k) = CTFRZ - sicemelt(i) = MIN( ssnow%wbice(i,k) * soil%zse_vec(i,k) * Cdensity_ice, & - ( ssnow%tgg(i,k) - tgg_tmp(i,k) ) * ssnow%gammzz(i,k) / Chlf ) + delta_ice_vol(i,k) = ssnow%wbice(i,k) - ssnow%wbice(i,k) = MAX( 0.0_r_2, ssnow%wbice(i,k) - sicemelt(i) & - / (soil%zse_vec(i,k) * Cdensity_ice) ) + delta_ice_vol(i,k) = min(delta_ice_vol(i,k), max(0._r_2,& + real(ssnow%tgg(i,k)-ssnow%otgg(i,k),r_2) * ssnow%gammzz(i,k) /& + (soil%zse_vec(i,k)*real(CHLF*Cdensity_ice,r_2)) ) ) - ssnow%gammzz(i,k) =MAX(soil%heat_cap_lower_limit(i,k), & - (1.0-soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k)& - + (ssnow%wb(i,k)-ssnow%wbice(i,k)) * REAL(Ccs_rho_wat,r_2) & - + ssnow%wbice(i,k) * REAL(Ccs_rho_ice,r_2)& - )*soil%zse_vec(i,k) + endif + end do + end do - IF (k .EQ. 1 .AND. ssnow%isflag(i) .EQ. 0) THEN - ssnow%gammzz(i,k) = ssnow%gammzz(i,k) + Ccgsnow * ssnow%snowd(i) - END IF - ssnow%tgg(i,k) = ssnow%tgg(i,k) - REAL(sicemelt(i)) & - * Chlf / REAL(ssnow%gammzz(i,k)) - !if for some reason end up here - ELSEIF( tgg_tmp(i,k) .GE. CTFRZ .AND. & - ssnow%tgg(i,k) > tgg_tmp(i,k) .AND. ssnow%wbice(i,k) > 0. ) THEN - - sicemelt(i) = MIN( ssnow%wbice(i,k) * soil%zse_vec(i,k) * Cdensity_ice, & - ( ssnow%tgg(i,k) - tgg_tmp(i,k) ) * ssnow%gammzz(i,k) / Chlf ) - - ssnow%wbice(i,k) = MAX( 0.0_r_2, ssnow%wbice(i,k) - sicemelt(i) & - / (soil%zse_vec(i,k) * Cdensity_ice) ) - - ssnow%gammzz(i,k) =MAX(soil%heat_cap_lower_limit(i,k), & - (1.0-soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k)& - + (ssnow%wb(i,k)-ssnow%wbice(i,k)) * REAL(Ccs_rho_wat,r_2) & - + ssnow%wbice(i,k) * REAL(Ccs_rho_ice,r_2)& - )*soil%zse_vec(i,k) - - IF (k .EQ. 1 .AND. ssnow%isflag(i) .EQ. 0) THEN - ssnow%gammzz(i,k) = ssnow%gammzz(i,k) + Ccgsnow * ssnow%snowd(i) - END IF - ssnow%tgg(i,k) = ssnow%tgg(i,k) - REAL(sicemelt(i)) & - * Chlf / REAL(ssnow%gammzz(i,k)) + DO k = 1, ms + DO i=1,mp + + if ( (ssnow%tgg(i,k) .lt. CTFRZ) .and. & + (delta_ice_vol(i,k) .gt. 0.0) ) then + + ssnow%wbice(i,k) = ssnow%wbice(i,k) + delta_ice_vol(i,k) + ssnow%wbliq(i,k) = ssnow%wbliq(i,k) - delta_ice_vol(i,k)*den_rat + + ssnow%gammzz(i,k) = max(soil%heat_cap_lower_limit(i,k), & + (1.0-soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k) & + + ssnow%wbliq(i,k) * REAL(Ccswat*Cdensity_liq,r_2) & + + ssnow%wbice(i,k) * REAL(Ccsice*Cdensity_ice,r_2)& + )*soil%zse_vec(i,k) + gammzz_snow(i,k) + ! MMY freezing releases heat + ssnow%tgg(i,k) = ssnow%tgg(i,k) + real( delta_ice_vol(i,k)*soil%zse_vec(i,k) *& + real(CHLF*Cdensity_ice,r_2) / ssnow%gammzz(i,k) ) + + + ssnow%wb(i,k) = ssnow%wbliq(i,k) + den_rat*ssnow%wbice(i,k) + + ssnow%wmliq(i,k) = ssnow%wbliq(i,k)*soil%zse_vec(i,k)*Cdensity_liq + ssnow%wmice(i,k) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*Cdensity_ice + ssnow%wmtot(i,k) = ssnow%wb(i,k) *soil%zse_vec(i,k)*Cdensity_liq + + elseif ((ssnow%tgg(i,k) .gt. CTFRZ) .and. & + delta_ice_vol(i,k) .gt. 0.0 ) then ! MMY + ! ssnow%wbice(i,k) .gt. 0.0 ) then ! MMY It's not right. + ! MMY When previous tgg > current tgg > 0, wbice won't melt + + ssnow%wbice(i,k) = ssnow%wbice(i,k) - delta_ice_vol(i,k) + ssnow%wbliq(i,k) = ssnow%wbliq(i,k) + delta_ice_vol(i,k)*den_rat + + ssnow%gammzz(i,k) = max(soil%heat_cap_lower_limit(i,k), & + (1.0-soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k) & + + ssnow%wbliq(i,k) * REAL(Ccswat*Cdensity_liq,r_2) & + + ssnow%wbice(i,k) * REAL(Ccsice*Cdensity_ice,r_2)& + )*soil%zse_vec(i,k) + gammzz_snow(i,k) + + ! MMY melting absorbs heat + ssnow%tgg(i,k) = ssnow%tgg(i,k) - real( delta_ice_vol(i,k)*soil%zse_vec(i,k) *& + real(CHLF*Cdensity_ice,r_2) / ssnow%gammzz(i,k) ) + + ssnow%wb(i,k) = ssnow%wbliq(i,k) + den_rat*ssnow%wbice(i,k) + + ssnow%wmliq(i,k) = ssnow%wbliq(i,k)*soil%zse_vec(i,k)*Cdensity_liq + ssnow%wmice(i,k) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*Cdensity_ice + ssnow%wmtot(i,k) = ssnow%wb(i,k) *soil%zse_vec(i,k)*Cdensity_liq END IF - !update the liq and ice volume and mass - ice_mass(i) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*REAL(Cdensity_ice,r_2) - liq_mass(i) = tot_mass(i) - ice_mass(i) - ssnow%wbliq(i,k) = liq_mass(i) / soil%zse_vec(i,k)/REAL(Cdensity_liq,r_2) - ssnow%wbice(i,k) = ice_mass(i) / soil%zse_vec(i,k)/REAL(Cdensity_ice,r_2) - ssnow%wb(i,k) = ssnow%wbliq(i,k) + ssnow%wbice(i,k) - END DO - END DO + + END DO + END DO END SUBROUTINE GWsoilfreeze - ! ---------------------------------------------------------------------------- + ! ----------------------------------------------------------------------------- ! - ! ---------------------------------------------------------------------------- + !! ----------------------------------------------------------------------------- ! SUBROUTINE remove_transGW(dels, soil, ssnow, canopy, veg) + + !*## Purpose + ! + !NOTE: this is only included because gw_model uses parameters XXX_vec !these are r_2. this breaks bitwise compatibility with trunk !if acceptable this routine does the same thing but with r_2 soil params @@ -252,7 +273,7 @@ SUBROUTINE remove_transGW(dels, soil, ssnow, canopy, veg) TYPE(soil_parameter_type), INTENT(INOUT) :: soil TYPE(veg_parameter_type), INTENT(INOUT) :: veg REAL(r_2), DIMENSION(mp,0:ms+1) :: diff - REAL(r_2), DIMENSION(mp) :: xx,xxd,evap_cur + REAL(r_2), DIMENSION(mp) :: xx,xxd REAL(r_2), DIMENSION(mp,ms) :: zse_mp_mm INTEGER :: k,i @@ -274,7 +295,7 @@ SUBROUTINE remove_transGW(dels, soil, ssnow, canopy, veg) IF (canopy%fevc(i) .GT. 0._r_2) THEN - xx(i) = canopy%fevc(i) * dels / Chl * veg%froot(i,k) + diff(i,k-1) + xx(i) = canopy%fevc(i) * dels / CHL * veg%froot(i,k) + diff(i,k-1) diff(i,k) = MAX(0._r_2,ssnow%wbliq(i,k)-soil%swilt_vec(i,k)) & * zse_mp_mm(i,k) xxd(i) = xx(i) - diff(i,k) @@ -300,7 +321,7 @@ SUBROUTINE remove_transGW(dels, soil, ssnow, canopy, veg) canopy%fevc = 0.0_r_2 END WHERE DO k = 1,ms - ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/zse_mp_mm(:,k) + ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/(soil%zse_vec(:,k)*m2mm) ENDDO ENDIF @@ -309,7 +330,7 @@ SUBROUTINE remove_transGW(dels, soil, ssnow, canopy, veg) DO i=1,mp ssnow%wmliq(i,k) = ssnow%wbliq(i,k)*zse_mp_mm(i,k)!mass ssnow%wmtot(i,k) = ssnow%wmliq(i,k) + ssnow%wmice(i,k) !mass - ssnow%wb(i,k) = ssnow%wbliq(i,k) + ssnow%wbice(i,k) !volume + ssnow%wb(i,k) = ssnow%wbliq(i,k) + den_rat * ssnow%wbice(i,k) !volume ! MMY END DO END DO @@ -317,16 +338,19 @@ SUBROUTINE remove_transGW(dels, soil, ssnow, canopy, veg) END SUBROUTINE remove_transGW -!============================================================================== +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! MD GW code from here on ! +!!!!!!!!!!!!!!MD GW code from here on!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !---------------------------------------------------------------------- -!============================================================================== +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------- SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call ) + !* Calculate surface runoff + USE cable_common_module, ONLY : gw_params,cable_user - + !*## Purpose + ! Calculate IMPLICIT NONE REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow ! soil+snow variables @@ -361,17 +385,19 @@ SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call ) END IF !amount of ice in surface layer DO i = 1,mp - efpor(i) = MAX(0.001_r_2, soil%ssat_vec(i,1) - ssnow%wbice(i,1)) - icemass = ssnow%wbice(i,1) * dzmm - liqmass = (ssnow%wb(i,1)-ssnow%wbice(i,1)) * dzmm + efpor(i) = max(0.01_r_2, soil%ssat_vec(i,1) - den_rat*ssnow%wbice(i,1)) + icemass = ssnow%wmice(i,1) + liqmass = ssnow%wmliq(i,1) + totmass = MAX(liqmass+icemass,REAL(1e-2,r_2)) - icef(i) = MAX(0._r_2,MIN(1._r_2,gw_params%IceBeta*icemass / totmass)) + icef(i) = max(0._r_2,min(1._r_2,icemass / totmass)) END DO !sat fraction assuming topo controlled subgrid soil moisture distribution !called from cable_canopy for srf wet fraction alrady - !call saturated_fraction(ssnow,soil,veg) + call saturated_fraction(ssnow,soil,veg) ! give ssnow%satfrac value before it is used + !srf frozen fraction. should be based on topography DO i = 1,mp fice = (EXP(-gw_params%IceAlpha*(1._r_2-icef(i)))-EXP(-gw_params%IceAlpha))/(1._r_2-EXP(-gw_params%IceAlpha)) @@ -382,9 +408,16 @@ SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call ) DO i=1,mp tmpa = ssnow%wbliq(i,1) / efpor(i) tmpb = MAX( (tmpa-satfrac_liqice(i))/MAX(0.01_r_2,(1._r_2-satfrac_liqice(i))), 0._r_2) - tmpa = -2._r_2*soil%bch_vec(i,1)*soil%sucs_vec(i,1)/dzmm + tmpa = -2._r_2*soil%bch_vec(i,1)*soil%sucs_vec(i,1)/soil%zse_vec(i,1)/1000._r_2 qinmax = (1._r_2 + tmpa*(tmpb-1._r_2))*soil%hyds_vec(i,1)*EXP(-gw_params%hkrz*(0.5*dzmm/1000.0_r_2-gw_params%zdepth)) - + ! qinmax = (1._r_2 + tmpa*(tmpb-1._r_2))*soil%hyds_vec(i,1) ! MMY@23Apr2023, keep this commented line and the line above + ! *EXP(-gw_params%hkrz*(0.5*dzmm/1000.0_r_2-gw_params%zdepth)) + ! should be in the eq if considering hk changes with soil depth + ! However, *EXP(-gw_params%hkrz*(0.5*dzmm/1000.0_r_2-gw_params%zdepth)) may only + ! make sense when soil properties are uniform in NO-GW codes. + ! In CABLE-GW soil properties vary with depth, EXP(-gw_params%hkrz*(0.5*dzmm/ + ! 1000.0_r_2-gw_params%zdepth)) complicates things. I think it should be + ! removed everywhere. ssnow%rnof1(i) = satfrac_liqice(i) * ssnow%fwtop(i) + & (1._r_2-satfrac_liqice(i))*MAX((ssnow%fwtop(i)-qinmax) , 0._r_2) @@ -403,7 +436,9 @@ SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call ) !--- glacier formation rnof5= 0. - IF (sli_call .OR. cable_runtime%UM) THEN +!$ IF (sli_call .OR. cable_runtime%UM) THEN +! added by rk4417 to conform with MMY modifications + IF (sli_call .OR. cable_runtime%UM .OR. cable_user%gw_model) THEN ! FEEDBACK (MMY asks: why cable_user%gw_model=True doesn't need to consider snow melting) --rk4417 nglacier = 0 ELSE nglacier = 2 @@ -418,7 +453,7 @@ SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call ) !---- change local tg to account for energy - clearly not best method WHERE( ssnow%isflag == 0 ) smasstot = 0.0 - ssnow%tgg(:,1) = ssnow%tgg(:,1) - rnof5 * Chlf & + ssnow%tgg(:,1) = ssnow%tgg(:,1) - rnof5 * CHLF & / REAL( ssnow%gammzz(:,1) ) ssnow%snowd = ssnow%snowd - rnof5 ELSEWHERE @@ -447,138 +482,107 @@ SUBROUTINE ovrlndflx (dels, ssnow, soil,veg, canopy,sli_call ) END SUBROUTINE ovrlndflx - - - -!============================================================================== - SUBROUTINE simple_wtd(ssnow, soil, veg) - !This was only for testing purposes - IMPLICIT NONE - TYPE (soil_snow_type), INTENT(INOUT) :: ssnow ! soil and snow variables - TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters - TYPE (veg_parameter_type), INTENT(INOUT) :: veg - - REAL(r_2), DIMENSION(mp) :: fz, wmean,ztot - REAL(r_2), DIMENSION(mp,ms) :: stot - INTEGER :: k,i - - DO i=1,mp - wmean(i) = 0._r_2 - fz(i) = 5._r_2 - ztot(i) = 0._r_2 - stot(i,:) = (ssnow%wb(i,:)-soil%watr(i,:)) / (soil%ssat_vec(i,:)-soil%watr(i,:)) - END DO - DO k = 1, ms - DO i=1,mp - wmean(i) = wmean(i) + stot(i,k)*soil%zse(k)*1000._r_2 - ztot(i) = ztot(i) + soil%zse(k)*1000._r_2 - END DO - END DO - - DO i=1,mp - wmean(i) = wmean(i) + ssnow%GWwb(i)/soil%GWssat_vec(i) * soil%GWdz(i)*1000._r_2 - ztot(i) = ztot(i) + soil%GWdz(i)*1000._r_2 - - ssnow%wtd(i) = MIN(200000._r_2, fz(i) * (ztot(i) - wmean(i))) - END DO - - END SUBROUTINE simple_wtd -!============================================================================== - - -!============================================================================== - - !---------------------------------------------------------------------- - ! SUBROUTINE iterative_wtd + + SUBROUTINE iterative_wtd (ssnow, soil, veg, include_aquifer) ! + !*## Purpose ! Iteratively calcs the water table depth by equating the mass of water in the ! soil column to the mass of a hydrostatic column inegrated from the surface to the ! water table depth - ! - SUBROUTINE iterative_wtd (ssnow, soil, veg, include_aquifer) + ! ChatGPT: "The water table depth in soil can be calculated iteratively by equating the mass of water in the soil + ! column to the mass of a hydrostatic column integrated from the surface to the water table depth. This + ! method is based on the assumption that the water in the soil column is in hydrostatic equilibrium, + ! which means that the pressure of the water at any point in the column is equal to the weight of the water + ! above that point. + ! The iterative process involves starting with an initial estimate of the water table depth and then + ! calculating the mass of water in the soil column using that depth. This mass is then compared to the + ! mass of the hydrostatic column integrated from the surface to the estimated water table depth. If the + ! two masses are not equal, the estimated water table depth is adjusted and the process is repeated until + ! convergence is achieved. + ! The calculation of the mass of the hydrostatic column requires knowledge of the soil properties, such as + ! porosity and permeability, as well as the water content of the soil. These parameters can be obtained through + ! laboratory testing or estimated from field observations. + ! Once the water table depth is determined, it can be used to assess the groundwater availability and potential + ! impacts on nearby structures or ecosystems." + ! Method: + ! I cannot find where this subroutine come from. It is not from CLM5 or CLM4.5. It also looks not from SIMGM model in Niu et al 2007 + IMPLICIT NONE TYPE (soil_snow_type), INTENT(INOUT) :: ssnow ! soil and snow variables TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters TYPE (veg_parameter_type), INTENT(INOUT) :: veg LOGICAL, INTENT(IN) :: include_aquifer !use GWwb or only wb to find wtd? - - !Local vars - REAL(r_2), DIMENSION(mp,ms) :: dzmm_mp,tmp_def - REAL(r_2), DIMENSION(0:ms) :: zimm - REAL(r_2), DIMENSION(ms) :: zmm - REAL(r_2), DIMENSION(mp) :: GWzimm,temp REAL(r_2), DIMENSION(mp) :: def,defc,total_depth_column - REAL(r_2) :: deffunc,tempa,tempb,derv,calc,tmpc - REAL(r_2), DIMENSION(mp) :: invB,Nsucs_vec !inverse of C&H B,Nsucs_vec - INTEGER :: k,i,wttd,jlp - - !make code cleaner define these here - invB = 1._r_2/soil%bch_vec(:,ms) !1 over C&H B - Nsucs_vec = soil%sucs_vec(:,ms) !psi_saturated mm - dzmm_mp = REAL(SPREAD((soil%zse(:)) * 1000.0,1,mp),r_2) !layer thickness mm - zimm(0) = 0.0_r_2 !depth of layer interfaces mm - - !total depth of soil column - DO k=1,ms - zimm(k) = zimm(k-1) + soil%zse(k)*1000._r_2 - END DO + REAL(r_2), DIMENSION(mp) :: lam,Nsucs_vec !inverse of C&H B,Nsucs_vec + INTEGER :: k,i,jlp - def(:) = 0._r_2 + lam(:) = 1._r_2/soil%bch_vec(:,ms) !1 over C&H B ! replaced block above as per MMY -- rk4417 + Nsucs_vec(:) = abs(soil%sucs_vec(:,ms)) !psi_saturated mm + ! MMY@23Apr2023, add aquifer thickness to the total depth and calcuate aquifer water content IF (include_aquifer) THEN !do we include the aquifer in the calculation of wtd? DO i=1,mp - total_depth_column(i) = zimm(ms) + soil%GWdz(i)*1000._r_2 - def(i) = def(i) + MAX(0._r_2,soil%GWssat_vec(i)-ssnow%GWwb(i))*soil%GWdz(i)*1000._r_2 + total_depth_column(i) = soil%GWdz(i)*m2mm + def(i) = max(0._r_2,soil%GWssat_vec(i)-ssnow%GWwb(i))*soil%GWdz(i)*m2mm END DO - + + ELSE + def(:) = 0._r_2 + total_depth_column(:) = 0._r_2 END IF + !total depth of soil column + do k=1,ms + do i=1,mp + total_depth_column(i) = total_depth_column(i) + soil%zse_vec(i,k)*m2mm + end do + end do + !comute the total mass away from full saturation DO k=1,ms DO i=1,mp def(i) = def(i) + & - MAX(0._r_2,(soil%ssat_vec(i,k)-(ssnow%wbliq(i,k)+ssnow%wbice(i,k)))*dzmm_mp(i,k)) + max(0._r_2,(soil%ssat_vec(i,k)-(ssnow%wbliq(i,k)+den_rat*ssnow%wbice(i,k)))*soil%zse_vec(i,k)*m2mm) END DO !mp END DO !ms !find the deficit if the water table is at the bottom of the soil column DO i=1,mp - defc(i) = (soil%ssat_vec(i,ms))*(total_depth_column(i)+Nsucs_vec(i)/(1._r_2-invB(i))* & - (1._r_2-((Nsucs_vec(i)+total_depth_column(i))/Nsucs_vec(i))**(1._r_2-invB(i)))) + defc(i) = (soil%ssat_vec(i,ms))*(total_depth_column(i)+Nsucs_vec(i)/(1._r_2-lam(i))* & + (1._r_2-((Nsucs_vec(i)+total_depth_column(i))/Nsucs_vec(i))**(1._r_2-lam(i)))) defc(i) = MAX(0.1_r_2,defc(i)) - - !initial guess at wtd - ssnow%wtd(:) = total_depth_column(:)*def(:)/defc(:) END DO - + + !initial guess at wtd ! taken out of do loop as per MMY -- rk4417 + ssnow%wtd(:) = total_depth_column(:)*def(:)/defc(:) !use newtons method to solve for wtd, note this assumes homogenous column but !that is ok DO i=1,mp IF ((soil%isoilm(i) .NE. 9) .AND. (veg%iveg(i) .NE. 16)) THEN - IF (defc(i) > def(i)) THEN !iterate tfor wtd + IF (defc(i) > def(i)) THEN !iterate tfor wtd ! MMY@23Apr2023 if soil deficit lower than total soil deficit jlp=0 mainloop: DO tempa = 1.0_r_2 - tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(-invB(i)) + tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(-lam(i)) derv = (soil%ssat_vec(i,ms))*(tempa-tempb) + & soil%ssat_vec(i,ms) IF (ABS(derv) .LT. REAL(1e-8,r_2)) derv = SIGN(REAL(1e-8,r_2),derv) tempa = 1.0_r_2 - tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(1._r_2-invB(i)) + tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(1._r_2-lam(i)) deffunc = (soil%ssat_vec(i,ms))*(ssnow%wtd(i) +& - Nsucs_vec(i)/(1-invB(i))* & + Nsucs_vec(i)/(1-lam(i))* & (tempa-tempb)) - def(i) calc = ssnow%wtd(i) - deffunc/derv @@ -600,22 +604,21 @@ SUBROUTINE iterative_wtd (ssnow, soil, veg, include_aquifer) END DO mainloop !defc .gt. def - ELSEIF (defc(i) .LT. def(i)) THEN + ELSEIF (defc(i) .LT. def(i)) THEN ! MMY@23Apr2023 if soil deficit larger than total soil deficit jlp=0 mainloop2: DO tmpc = Nsucs_vec(i)+ssnow%wtd(i)-total_depth_column(i) - tempa = (ABS(tmpc/Nsucs_vec(i)))**(-invB(i)) - tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(-invB(i)) + tempa = (abs(tmpc/Nsucs_vec(i)))**(-lam(i)) + tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(-lam(i)) derv = (soil%ssat_vec(i,ms))*(tempa-tempb) IF (ABS(derv) .LT. REAL(1e-8,r_2)) derv = SIGN(REAL(1e-8,r_2),derv) - - tempa = (ABS((Nsucs_vec(i)+ssnow%wtd(i)-total_depth_column(i))/Nsucs_vec(i)))**(1._r_2-invB(i)) - tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(1._r_2-invB(i)) + tempa = (abs((Nsucs_vec(i)+ssnow%wtd(i)-total_depth_column(i))/Nsucs_vec(i)))**(1._r_2-lam(i)) + tempb = (1._r_2+ssnow%wtd(i)/Nsucs_vec(i))**(1._r_2-lam(i)) deffunc = (soil%ssat_vec(i,ms))*(total_depth_column(i) +& - Nsucs_vec(i)/(1._r_2-invB(i))*(tempa-tempb))-def(i) + Nsucs_vec(i)/(1._r_2-lam(i))*(tempa-tempb))-def(i) calc = ssnow%wtd(i) - deffunc/derv IF ((ABS(calc-ssnow%wtd(i))) .LE. wtd_uncert) THEN @@ -655,16 +658,19 @@ SUBROUTINE iterative_wtd (ssnow, soil, veg, include_aquifer) END SUBROUTINE iterative_wtd - !------------------------------------------------------------------------- - ! SUBROUTINE smoistgw (fwtop,dt,ktau,ssnow,soil,prin) - ! solves the modified richards equation (Zeng and Decker 2009) to find - ! vertical mocement of soil water. Bottom boundary condition is determined - ! using a single layer groundwater module - ! - SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) - USE cable_common_module -USE trimb_mod, ONLY : trimb + SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) + + !* Solve the modified Richards equation + ! [Zeng and Decker 2009](https://doi.org/10.1175/2008JHM1011.1) + ! to find vertical movement of soil water. The modified Richards equation considers + ! soil moisture distribution at hydrostatic equilibrium to avoid the truncation + ! errors due to finite grid spacing used in the partial differential equation. + ! Bottom boundary condition is determined using a single-layer groundwater module. + + USE cable_common_module + USE trimb_mod, ONLY : trimb ! inserted by rk4417 - phase2 + IMPLICIT NONE REAL, INTENT(IN) :: dels ! time step size (s) @@ -693,33 +699,30 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) REAL(r_2), DIMENSION(mp) :: dqodw0 REAL(r_2), DIMENSION(mp) :: dqodw1,dqodw2 REAL(r_2), DIMENSION(mp) :: s1,s2,tmpi,temp0,voleq1,tempi - REAL(r_2), DIMENSION(ms) :: dzmm - REAL(r_2), DIMENSION(mp,ms) :: dzmm_mp - REAL(r_2), DIMENSION(0:ms+1) :: zimm - REAL(r_2), DIMENSION(ms) :: zmm - REAL(r_2), DIMENSION(mp) :: GWzimm,xs,zaq,s_mid,GWdzmm + + REAL(r_2), DIMENSION(mp,0:ms+1) :: zimm + REAL(r_2), DIMENSION(mp,ms) :: zmm + + REAL(r_2), DIMENSION(mp) :: GWzimm,xs,zaq,s_mid REAL(r_2), DIMENSION(mp) :: xs1,GWmsliq!xsi !mass (mm) of liquid over/under saturation, mass of aquifer water REAL(r_2) :: xsi REAL(r_2), DIMENSION(mp,ms+1) :: del_wb !MD DEBUG VARS INTEGER :: imp,ims,k_drain - !make code cleaner define these here - dzmm(:) = 1000.0_r_2 * REAL(soil%zse(:),r_2) - DO i=1,mp - dzmm_mp(i,:) = dzmm(:) - END DO - - zimm(0) = 0._r_2 + zimm(:,0) = 0._r_2 + ! MMY@23Apr2023 calcuate soil layer boundary DO k=1,ms - zimm(k) = zimm(k-1) + dzmm(k) + zimm(:,k) = zimm(:,k-1) + soil%zse_vec(:,k)*m2mm END DO - zmm(1:ms) = zimm(0:(ms-1)) + 0.5_r_2*dzmm(1:ms) + + ! MMY@23Apr2023 calcuate the middle point of each soil layer + zmm(:,1:ms) = zimm(:,1:ms) - 0.5*soil%zse_vec(:,1:ms)*m2mm DO i=1,mp - GWdzmm(i) = REAL(soil%GWdz(i),r_2)*1000._r_2 - GWzimm(i) = zimm(ms)+GWdzmm(i) - zaq(i) = zimm(ms) + 0.5_r_2*GWdzmm(i) + ! MMY@23Apr2023 calcuate the bottom boundary of aquifer and the middle point of the aquifer + GWzimm(i) = zimm(i,ms)+soil%GWdz(i)*m2mm + zaq(i) = zimm(i,ms) + 0.5_r_2*soil%GWdz(i)*m2mm END DO ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -734,17 +737,20 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) END DO END DO - !equilibrium water content - CALL calc_equilibrium_water_content(ssnow,soil) - - CALL calc_soil_hydraulic_props(ssnow,soil,veg) + !> soil hydraulic state/props + CALL calc_soil_hydraulic_props(ssnow,soil,veg) + + !> equilibrium water content + CALL calc_equilibrium_water_content(ssnow,soil) - CALL subsurface_drainage(ssnow,soil,veg,dzmm) + !> calculate the subsurface runoff + CALL subsurface_drainage(ssnow,soil,veg) + !* Calculate soil moisture for soil column and aquifer k = 1 !top soil layer DO i=1,mp qin(i) = ssnow%sinfil(i) - den(i) = (zmm(k+1)-zmm(k)) + den(i) = 0.5*(soil%zse_vec(i,k)+soil%zse_vec(i,k+1) )*m2mm dne(i) = (ssnow%zq(i,k+1)-ssnow%zq(i,k)) num(i) = (ssnow%smp(i,k+1)-ssnow%smp(i,k)) - dne(i) qout(i) = -ssnow%hk(i,k)*num(i)/den(i) @@ -752,18 +758,19 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) dqodw2(i) = -( ssnow%hk(i,k)*ssnow%dsmpdw(i,k+1) + num(i)*ssnow%dhkdw(i,k))/den(i) rt(i,k) = qin(i) - qout(i) at(i,k) = 0._r_2 - bt(i,k) = dzmm(k)/dels + dqodw1(i) + bt(i,k) = m2mm*soil%zse_vec(i,k)/dels + dqodw1(i) ct(i,k) = dqodw2(i) END DO DO k = 2, ms - 1 !middle soil layers DO i=1,mp - den(i) = (zmm(k) - zmm(k-1)) + + den(i) = 0.5*(soil%zse_vec(i,k)+soil%zse_vec(i,k-1) )*m2mm dne(i) = (ssnow%zq(i,k)-ssnow%zq(i,k-1)) num(i) = (ssnow%smp(i,k)-ssnow%smp(i,k-1)) - dne(i) qin(i) = -ssnow%hk(i,k-1)*num(i)/den(i) dqidw0(i) = -(-ssnow%hk(i,k-1)*ssnow%dsmpdw(i,k-1) + num(i)*ssnow%dhkdw(i,k-1))/den(i) dqidw1(i) = -( ssnow%hk(i,k-1)*ssnow%dsmpdw(i,k) + num(i)*ssnow%dhkdw(i,k-1))/den(i) - den(i) = (zmm(k+1)-zmm(k)) + den(i) = 0.5*(soil%zse_vec(i,k)+soil%zse_vec(i,k+1) )*m2mm dne(i) = (ssnow%zq(i,k+1)-ssnow%zq(i,k)) num(i) = (ssnow%smp(i,k+1)-ssnow%smp(i,k)) - dne(i) qout(i) = -ssnow%hk(i,k)*num(i)/den(i) @@ -771,20 +778,22 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) dqodw2(i) = -( ssnow%hk(i,k)*ssnow%dsmpdw(i,k+1) + num(i)*ssnow%dhkdw(i,k))/den(i) rt(i,k) = qin(i) - qout(i) at(i,k) = -dqidw0(i) - bt(i,k) = dzmm(k)/dels - dqidw1(i) + dqodw1(i) + bt(i,k) = m2mm*soil%zse_vec(i,k)/dels - dqidw1(i) + dqodw1(i) + ! MMY ??? why the sign od bt is different to CLM5's ct(i,k) = dqodw2(i) END DO END DO k = ms !Bottom soil layer DO i=1,mp - den(i) = (zmm(k) - zmm(k-1)) + + den(i) = 0.5*(soil%zse_vec(i,k)+soil%zse_vec(i,k-1) )*m2mm dne(i) = (ssnow%zq(i,k)-ssnow%zq(i,k-1)) num(i) = (ssnow%smp(i,k)-ssnow%smp(i,k-1)) - dne(i) qin(i) = -ssnow%hk(i,k-1)*num(i)/den(i) dqidw0(i) = -(-ssnow%hk(i,k-1)*ssnow%dsmpdw(i,k-1) + num(i)*ssnow%dhkdw(i,k-1))/den(i) dqidw1(i) = -( ssnow%hk(i,k-1)*ssnow%dsmpdw(i,k) + num(i)*ssnow%dhkdw(i,k-1))/den(i) - den(i) = zaq(i) - zmm(k) + den(i) = zaq(i) - zmm(i,k) dne(i) = (ssnow%GWzq(i)-ssnow%zq(i,k)) num(i) = (ssnow%GWsmp(i)-ssnow%smp(i,k)) - dne(i) qout(i) = 0._r_2 @@ -792,32 +801,32 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) dqodw2(i) = 0._r_2 rt(i,k) = qin(i) - qout(i) at(i,k) = -dqidw0(i) - bt(i,k) = dzmm(k)/dels - dqidw1(i) + dqodw1(i) + bt(i,k) = m2mm*soil%zse_vec(i,k)/dels - dqidw1(i) + dqodw1(i) ct(i,k) = dqodw2(i) END DO - CALL aquifer_recharge(dels,ssnow,soil,veg,zaq,zmm,dzmm) + CALL aquifer_recharge(dels,ssnow,soil,veg) - CALL trimb(at,bt,ct,rt,ms) !use the defulat cable tridiag solution + CALL trimb(at,bt,ct,rt,ms) !use the default cable tridiag solution DO k=1,ms DO i=1,mp - ssnow%wbliq(i,k) = ssnow%wbliq(i,k) + rt(i,k) - ssnow%qhlev(i,k)*dels/dzmm(k) !volutermic liquid + ssnow%wbliq(i,k) = ssnow%wbliq(i,k) + rt(i,k) - ssnow%qhlev(i,k)*dels/(m2mm*soil%zse_vec(i,k)) !volutermic liquid END DO END DO DO i=1,mp - ssnow%wbliq(i,ms) = ssnow%wbliq(i,ms) - ssnow%Qrecharge(i)*dels/dzmm(ms) + ssnow%wbliq(i,ms) = ssnow%wbliq(i,ms) - ssnow%Qrecharge(i)*dels/(m2mm*soil%zse_vec(i,ms)) END DO DO i=1,mp - ssnow%GWwb(i) = ssnow%GWwb(i) + (ssnow%Qrecharge(i)-ssnow%qhlev(i,ms+1))*dels/GWdzmm(i) + ssnow%GWwb(i) = ssnow%GWwb(i) + (ssnow%Qrecharge(i)-ssnow%qhlev(i,ms+1))*dels/(m2mm*soil%GWdz(i)) END DO !determine the available pore space !volumetric DO k=1,ms DO i=1,mp - eff_por(i,k) = MAX(0._r_2, soil%ssat_vec(i,k) - ssnow%wbice(i,k) ) + eff_por(i,k) = max(0._r_2, soil%ssat_vec(i,k) - den_rat * ssnow%wbice(i,k) ) END DO END DO @@ -825,24 +834,24 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) xsi = 0._r_2 IF (ssnow%GWwb(i) .GT. soil%GWssat_vec(i)) THEN - xsi = (ssnow%GWwb(i) - soil%GWssat_vec(i))*GWdzmm(i) + xsi = (ssnow%GWwb(i) - soil%GWssat_vec(i))*m2mm*soil%GWdz(i) ssnow%GWwb(i) = soil%GWssat_vec(i) END IF DO k=1,ms IF (ssnow%wbliq(i,k) .GT. eff_por(i,k)) THEN - xsi = xsi + (ssnow%wbliq(i,k) - eff_por(i,k))*dzmm(k) + xsi = xsi + (ssnow%wbliq(i,k) - eff_por(i,k))*(m2mm*soil%zse_vec(i,k)) ssnow%wbliq(i,k) = eff_por(i,k) END IF END DO DO k = ms,1,-1 !loop from bottom to top adding extra water to each layer IF (xsi .GT. 0._r_2) THEN - IF (xsi .LT. (eff_por(i,k)-ssnow%wbliq(i,k))*dzmm(k)) THEN - ssnow%wbliq(i,k) = ssnow%wbliq(i,k) + xsi/dzmm(k) + if (xsi .lt. (eff_por(i,k)-ssnow%wbliq(i,k))*(m2mm*soil%zse_vec(i,k))) then + ssnow%wbliq(i,k) = ssnow%wbliq(i,k) + xsi/(m2mm*soil%zse_vec(i,k)) xsi = 0._r_2 ELSE - xsi = xsi - (eff_por(i,k) - ssnow%wbliq(i,k))*dzmm(k) + xsi = xsi - (eff_por(i,k) - ssnow%wbliq(i,k))*(m2mm*soil%zse_vec(i,k)) ssnow%wbliq(i,k) = eff_por(i,k) END IF END IF @@ -853,23 +862,30 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) xsi = 0._r_2 END IF + ! MMY use watr to replace volwatmin to avoid questionable smp or hk (when wbliq Soil is too dry, aquifer extracts water from lateral drainage, found in SUBROUTINE smoistgw" + ! MMY@23Apr2023, it's an alarm for doubtful low aquifer moisture END IF END DO @@ -879,14 +895,15 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) !update mass variables ssnow%wmliq(i,k) = ssnow%wbliq(i,k) * & - soil%zse_vec(i,k)*Cdensity_liq + soil%zse_vec(i,k)*real(Cdensity_liq,r_2) ssnow%wmice(i,k) = ssnow%wbice(i,k) * & - soil%zse_vec(i,k)*Cdensity_ice - ssnow%wb(i,k) = ssnow%wbliq(i,k) + ssnow%wbice(i,k) + soil%zse_vec(i,k)*real(Cdensity_ice,r_2) + ssnow%wb(i,k) = ssnow%wbliq(i,k) + den_rat * ssnow%wbice(i,k) ssnow%wmtot(i,k) = ssnow%wmliq(i,k) + ssnow%wmice(i,k) END DO END DO + !> give value to subsurface runoff DO i=1,mp ssnow%rnof2(i) = ssnow%qhz(i) !rnof2 is default cable deep runoff var END DO !mp loop @@ -894,7 +911,25 @@ SUBROUTINE smoistgw (dels,ktau,ssnow,soil,veg,canopy) END SUBROUTINE smoistgw - + SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) + !*## Purpose + ! + ! This is the main subroutine for Mark Decker's new hydrology scheme. + ! This subroutine should be called in subroutine cbm but not yet. + ! + !## Method + ! + ! The subroutine takes a float as input, prints its value then adds 2/3, + ! prints it and returns it. + ! + ! The equation used is: + ! \[ myarg = myarg + !frac{2.}{3.} \] + ! + !## References + ! + ! The guidelines for documentation can be found in: + ! [CABLE developer guide](https://cable-lsm.github.io/CABLE/developer_guide/doc_guide/science_doc/) + ! ! Inputs: ! dt_in - time step in sec @@ -909,11 +944,10 @@ END SUBROUTINE smoistgw ! ivegt - vegetation type ! Output ! ssnow - SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) - USE cable_IO_vars_module, ONLY: wlogn USE cable_common_module -USE snow_processes_soil_thermal_mod, ONLY : snow_processes_soil_thermal + USE snow_processes_soil_thermal_mod, ONLY : snow_processes_soil_thermal ! inserted by rk4417 - phase2 + REAL , INTENT(IN) :: dels ! integration time step (s) TYPE(soil_parameter_type), INTENT(INOUT) :: soil TYPE(soil_snow_type) , INTENT(INOUT) :: ssnow @@ -923,11 +957,11 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) TYPE (balances_type) , INTENT(INOUT) :: bal INTEGER :: k,i - REAL, DIMENSION(mp) :: snowmlt - REAL, DIMENSION(mp) :: GWwb_ic - REAL, DIMENSION(mp,ms) :: tgg_old - REAL, DIMENSION(mp) :: tggsn_old,wbtot_ic,del_wbtot + ! REAL, DIMENSION(mp) :: snowmlt ! MMY@Nov2022 not really use + ! REAL, DIMENSION(mp) :: GWwb_ic ! MMY@Nov2022 not really use + REAL, DIMENSION(mp) :: tggsn_old,del_wbtot ! wbtot_ic ! MMY@Nov2022 not really use REAL(r_2), DIMENSION(mp) :: xx + REAL(r_2), DIMENSION(mp,ms) :: gammzz_snow REAL :: zsetot INTEGER, SAVE :: ktau =0 REAL(r_2) :: wb_lake_T, rnof2_T @@ -935,21 +969,34 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) LOGICAL, SAVE :: first_gw_hydro_call = .TRUE. use_sli = .FALSE. - ktau = ktau +1 - zsetot = SUM(soil%zse) ssnow%tggav = 0. - DO k = 1, ms + !>## Order of procedure + !> 1. Set the physic constants + !> 2. Calculate average soil temp for soil column + DO k = 1, ms ssnow%tggav = ssnow%tggav + soil%zse(k)*ssnow%tgg(:,k)/zsetot - END DO + !> 3. If it is the first time step, define the ratio of ice density to liquid water density + ! and constrains wb and wbice + if (first_gw_hydro_call) then + den_rat = real(Cdensity_ice/Cdensity_liq,r_2) ! MMY @Nov2022 + ! call set_den_rat() ! MMY @Nov2022 + DO k = 1, ms + ssnow%wb(:,k) = MIN(soil%ssat_vec(:,k), MAX(real(ssnow%wb(:,k)), soil%watr(:,k))) + ssnow%wbice(:,k) = MIN(real(ssnow%wb(:,k))/den_rat, real(ssnow%wbice(:,k))) + END DO + end if + !> 4. set a default snow depth in the top snow layer IF( cable_runtime%offline .OR. cable_runtime%mk3l ) ssnow%t_snwlr = 0.05_r_2 + !> 5. initialize hydrology variables + ssnow%otgg = ssnow%tgg ! initialize of ssnow%otgg DO i=1,mp ssnow%fwtop1(i) = 0.0 ssnow%fwtop2(i) = 0.0 @@ -958,119 +1005,151 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) ssnow%rnof1(i) = 0.0 ! initialise surface runoff ssnow%rnof2(i) = 0.0 ! initialise deep drainage ssnow%smelt(i) = 0.0 ! initialise snowmelt - ssnow%dtmlt(i,:) = 0.0 + ssnow%dtmlt(i,:) = 0.0 ! initialise water flux to the soil ssnow%osnowd(i) = ssnow%snowd(i) ! Scaling runoff to kg/m^2/s (mm/s) to match rest of the model ssnow%sinfil(i) = 0.0 ssnow%qhz(i) = 0.0 END DO + !> 6. set heat cap lower limit IF (cable_user%soil_thermal_fix) THEN - soil%heat_cap_lower_limit(:,:) = 0.01 !never allow /0 +!$ soil%heat_cap_lower_limit(:,:) = 0.01 !never allow /0 ! FEEDBACK (MMY asks: I guess we should delete this line since someone changed it on purpose) --rk4417 + soil%heat_cap_lower_limit(:,:) = 0._r_2 !allow /0 to show bugs ! FEEDBACK (MMY asks: I don't know how it shows bugs, can anyone check it???) --rk4417 ELSE + print *, "MMY testing soil%css_vec(:,:) =",soil%css_vec(:,:), "soil%rhosoil_vec(:,:) = ", soil%rhosoil_vec(:,:) !MMY@18May2023 soil%heat_cap_lower_limit(:,:) = soil%css_vec(:,:) * soil%rhosoil_vec(:,:) END IF + !> 7. MMY??? IF( (.NOT.cable_user%cable_runtime_coupled ) .AND. (first_gw_hydro_call)) THEN - IF (cable_runtime%um) canopy%dgdtg = 0.0 ! RML added um condition - ! after discussion with BP - ! N.B. snmin should exceed sum of layer depths, i.e. .11 m - ssnow%wbtot = 0.0 - ssnow%wb(:,:) = MIN( soil%ssat_vec(:,:), MAX ( ssnow%wb(:,:), 0.5*soil%swilt_vec(:,:) ) ) - - DO k = 1, ms - - WHERE( ssnow%tgg(:,k) <= CTFRZ .AND. ssnow%wbice(:,k) <= 0.001*ssnow%wb(:,k) ) & - ssnow%wbice(:,k) = 0.5 * ssnow%wb(:,k) - - !WHERE( ssnow%tgg(:,k) < CTFRZ) & - ! ssnow%wbice(:,k) = 0.8 * ssnow%wb(:,k) - - END DO - - WHERE ( soil%isoilm .EQ. 9)! .and. ssnow%snowd .le. 0.1*max_glacier_snowd) - - ! permanent ice: fix hard-wired number in next version - ssnow%snowd = max_glacier_snowd - ssnow%osnowd = max_glacier_snowd - ssnow%tgg(:,1) = ssnow%tgg(:,1) - 1.0 - - END WHERE - - WHERE ( SPREAD(soil%isoilm,2,ms) .EQ. 9 ) - - ssnow%wb = 0.95 * soil%ssat_vec - ssnow%wbice = 0.95 * ssnow%wb - - END WHERE - + IF (cable_runtime%um) canopy%dgdtg = 0.0 ! RML added um condition + +! MMY@23Apr2023 delete this block, since without it CABLE-GW works well ! FEEDBACK (Please see MMY questions regarding commented block below) --rk4417 +! MMY@13Mar2023 the block below is an initialization for ssnow%wb, ssnow%wbice & ssnow%snowd (if soil%isoilm=9) +! It is duplicate to the section 3. I don't think we need to initializate ssnow%wb and ssnow%wbice again. +! However, I don't know whether we should keep ssnow%snowd = max_glacier_snowd for soil%isoilm=9. +! Please check with Claire and Jhan if possible. + +!$ ! MMY??? should we keep wb & wbice setting for soil%isoilm==9 (permanent ice) and snow covering region +!$ +!$ ! after discussion with BP ! block commented out by rk4417 as missing from MMY +!$ ! N.B. snmin should exceed sum of layer depths, i.e. .11 m ! need to discuss with MMY +!$ ssnow%wbtot = 0.0 +!$ ssnow%wb(:,:) = MIN( soil%ssat_vec(:,:), MAX ( ssnow%wb(:,:), 0.5*soil%swilt_vec(:,:) ) ) +!$ +!$ DO k = 1, ms +!$ +!$ WHERE( ssnow%tgg(:,k) <= CTFRZ .AND. ssnow%wbice(:,k) <= 0.001*ssnow%wb(:,k) ) & +!$ ssnow%wbice(:,k) = 0.5 * ssnow%wb(:,k) +!$ +!$ !WHERE( ssnow%tgg(:,k) < CTFRZ) & +!$ ! ssnow%wbice(:,k) = 0.8 * ssnow%wb(:,k) +!$ +!$ END DO +!$ WHERE ( soil%isoilm .EQ. 9)! .and. ssnow%snowd .le. 0.1*max_glacier_snowd) +!$ +!$ ! permanent ice: fix hard-wired number in next version +!$ ssnow%snowd = max_glacier_snowd +!$ ssnow%osnowd = max_glacier_snowd +!$ ssnow%tgg(:,1) = ssnow%tgg(:,1) - 1.0 +!$ +!$ END WHERE +!$ +!$ WHERE ( SPREAD(soil%isoilm,2,ms) .EQ. 9 ) +!$ +!$ ssnow%wb = 0.95 * soil%ssat_vec +!$ ssnow%wbice = 0.95 * ssnow%wb +!$ +!$ END WHERE +! END DELETE END IF - - tgg_old = ssnow%tgg - - !Start with wb and wbice. Need wbliq, wmliq,wmice,wmtot - !find the mass of ice and liq from the prognostic volumetric values + !> 8. Calculate snow heat capacity to the top soil + ! ??? MMY@Nov2022 gammzz_snow and ssnow%gammzz calculated here are useless, should be commented out ___ + gammzz_snow(:,:) = 0._r_2 ! MMY@23Apr2023 necessary to have this line + ! since only the top soil layer can have snow + ! and gammzz_snow in the other layers should + ! be 0 + do i=1,mp + ! ssnow%isflag = 0 one snow layer; =1 three snow layers + gammzz_snow(i,1) = real((1. - ssnow%isflag(i)) * Ccgsnow * ssnow%snowd(i),r_2) + end do + + !> 9. Find the mass and saturation fraction of ice and liq from the prognostic volumetric values (wb,wbice) + ! and calculate sucs_hys and smp_hys DO k=1,ms DO i=1,mp - ssnow%wbliq(i,k) = ssnow%wb(i,k) - ssnow%wbice(i,k) !liquid volume - ssnow%wmice(i,k) = ssnow%wbice(i,k)*REAL(Cdensity_ice*soil%zse(k),r_2) !ice mass - ssnow%wmliq(i,k) = ssnow%wbliq(i,k)*REAL(Cdensity_liq*soil%zse(k),r_2) !liquid mass + ssnow%wbliq(i,k) = ssnow%wb(i,k) - den_rat*ssnow%wbice(i,k) ! MMY !liquid volume + ssnow%wmice(i,k) = ssnow%wbice(i,k)*real(Cdensity_ice,r_2)*soil%zse_vec(i,k)!ice mass + ssnow%wmliq(i,k) = ssnow%wbliq(i,k)*real(Cdensity_liq,r_2)*soil%zse_vec(i,k) !liquid mass ssnow%wmtot(i,k) = ssnow%wmice(i,k) + ssnow%wmliq(i,k) !liq+ice mass + ssnow%wblf(i,k) = MAX(ssnow%wbliq(i,k)/soil%ssat_vec(i,k),0.01_r_2) ssnow%wbfice(i,k) = MAX(ssnow%wbice(i,k)/soil%ssat_vec(i,k),0._r_2) - + !do not pass with mpi + ssnow%sucs_hys(i,k) = ssnow%hys_fac(i,k)*soil%sucs_vec(i,k) + ssnow%smp_hys(i,k) = min(ssnow%smp_hys(i,k),-ssnow%sucs_hys(i,k)) END DO END DO - + !> 10. calculate heat capacity for each soil layer IF( first_gw_hydro_call ) THEN - - DO i=1,mp - ssnow%gammzz(i,1) = MAX(soil%heat_cap_lower_limit(i,1),& - (1.0-soil%ssat_vec(i,1))*& - soil%css_vec(i,1) * soil%rhosoil_vec(i,1) & - & + ssnow%wbliq(i,1) * Ccs_rho_wat & - & + ssnow%wbice(i,1) * Ccs_rho_ice ) * soil%zse(1) + & - & (1. - ssnow%isflag(i)) * Ccgsnow * ssnow%snowd(i) - - END DO - + DO k=1,ms + DO i=1,mp + ssnow%gammzz(i,k) = max(soil%heat_cap_lower_limit(i,k),& + (1.0-soil%ssat_vec(i,k))*& + soil%css_vec(i,k) * soil%rhosoil_vec(i,k) & + + ssnow%wbliq(i,k) * real(Ccswat*Cdensity_liq,r_2) & + + ssnow%wbice(i,k) * real(Ccsice*Cdensity_ice,r_2) )* & + soil%zse_vec(i,k) + gammzz_snow(i,k) + END DO + END DO ENDIF ! if(.NOT.cable_runtime_coupled) and first_gw_hydro_call + !> wbliq_old is used for hysteresis + ssnow%wbliq_old = ssnow%wbliq - DO i=1,mp - !initial water in the soil column - wbtot_ic(i) = SUM(ssnow%wbliq(i,:)*Cdensity_liq*soil%zse(:),1) + & - SUM(ssnow%wbice(i,:)*Cdensity_ice*soil%zse(:),1) + & - ssnow%GWwb(i)*soil%GWdz(i)*Cdensity_liq - - GWwb_ic(i) = ssnow%GWwb(i) - - END DO - - !improve hiding, call single soilsnow subroutine to do all the - !snow processes and thermal soil calculations - - CALL snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) - - !leave here for now, could move into soilsnow as well - CALL remove_transGW(dels, soil, ssnow, canopy, veg) !transpiration loss per soil layer - - CALL GWsoilfreeze(dels, soil, ssnow,tgg_old) + !> 10. snow processes and thermal soil + CALL snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) + + !> 11. transpiration loss per soil layer ! leave here for now, could move into soilsnow as well + CALL remove_transGW(dels, soil, ssnow, canopy, veg) + !> 12. Snow freezes and melts. + CALL GWsoilfreeze(dels, soil, ssnow) + ssnow%fwtop = canopy%precis/dels + ssnow%smelt/dels !water from canopy and snowmelt [mm/s] + !> 13. calculate water table depth CALL iterative_wtd (ssnow, soil, veg, .TRUE. ) - CALL ovrlndflx (dels, ssnow, soil, veg, canopy,use_sli ) !surface runoff, incorporate ssnow%pudsto? - - ssnow%sinfil = ssnow%fwtop - canopy%segg !canopy%fes/Chl !remove soil evap from throughfall - - CALL smoistgw (dels,ktau,ssnow,soil,veg,canopy) !vertical soil moisture movement. - - ! correction required for energy balance in online simulations + !> 14. surface runoff, incorporate ssnow%pudsto? + CALL ovrlndflx (dels, ssnow, soil, veg, canopy,use_sli ) + !> remove soil evap from throughfall + ssnow%sinfil = ssnow%fwtop - canopy%segg !canopy%fes/CHL + + !> 15. vertical soil moisture movement. + CALL smoistgw (dels,ktau,ssnow,soil,veg,canopy) + + ! _________ MMY@23Apr2023 need to check these hysteresis-related codes, I think they should show before smoistgw ______ + !* if Brook & Cosby water retention hysteresis is used, calculate smp_hys, wb_hys, + ! ssat_hys, hys_fac, sucs_hys, watr_hys... + IF (gw_params%BC_hysteresis) & + CALL swc_hyst_direction(soil,ssnow,veg) + + !call swc_smp_dsmpdw(soil,ssnow) + !> 16. soil matric potential + if (gw_params%BC_hysteresis) then + call brook_corey_hysteresis_swc_smp(soil,ssnow) + elseif (gw_params%HC_SWC) then + call hutson_cass_swc_smp(soil,ssnow) + else + call brook_corey_swc_smp(soil,ssnow) + end if + + !> 17. correction required for energy balance in online simulations IF( cable_runtime%um ) THEN !cls package - rewritten for flexibility @@ -1094,11 +1173,12 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) canopy%fess = canopy%fess + canopy%fes_cor ENDIF ENDIF - + + !> 18. update variables ! MMY??? where is pudsto calculated DO i=1,mp ssnow%pudsto(i) = 0.0 !no puddle ssnow%smelt(i) = ssnow%smelt(i)/dels !change units to mm/s. cable_driver then reverts back to mm - ssnow%runoff(i) = (ssnow%rnof1(i) + ssnow%rnof2(i))! *dels !cable_driver converts from mm/s to mm + ssnow%runoff(i) = (ssnow%rnof1(i) + ssnow%rnof2(i))!*dels !cable_driver converts from mm/s to mm !rnof1 and rnof2 are already in mm/s ! Set weighted soil/snow surface temperature ssnow%tss(i) = (1-ssnow%isflag(i))*ssnow%tgg(i,1) + ssnow%isflag(i)*ssnow%tggsn(i,1) @@ -1109,16 +1189,20 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) ssnow%GWwb(i)*soil%GWdz(i)*Cdensity_liq !for debug water balance. del_wbtot = fluxes = infiltration [though-evap] - trans - qhorz drainage - del_wbtot(i) = dels * (ssnow%sinfil(i) - ssnow%rnof2(i) - canopy%fevc(i) / Chl) - !set below to keep track of water imbalance within the GW module explicitly. also must change cable_checks - !ssnow%wbtot(i) = ssnow%wbtot(i)-wbtot_ic(i) + del_wbtot(i) = dels * (ssnow%sinfil(i) - ssnow%rnof2(i) - canopy%fevc(i) / CHL) END DO - + ! MMY take from subroutine soil_snow, hydraulic redistribution should be included + ! IF( redistrb ) & + ! CALL hydraulic_redistribution( dels, soil, ssnow, canopy, veg, met ) + first_gw_hydro_call=.false. + END SUBROUTINE soil_snow_gw SUBROUTINE calc_equilibrium_water_content(ssnow,soil) - !find layer mean soil moisture and potential at equilibrium with wtd + !* Calculate soil moisture and potential at hydrostatic equilibrium + ! with water table depth [Zeng and Decker 2009](https://doi.org/10.1175/2008JHM1011.1) + ! MMY unfinished IMPLICIT NONE @@ -1129,150 +1213,203 @@ SUBROUTINE calc_equilibrium_water_content(ssnow,soil) REAL(r_2), DIMENSION(ms) :: dzmm !layer thickness for single tile REAL(r_2), DIMENSION(mp) :: GWdzmm !aquifer thickness at each tile REAL(r_2), DIMENSION(mp) :: GWzimm !aquifer layer interface depth - REAL(r_2), DIMENSION(0:ms) :: zimm !layer interface depth in mm - REAL(r_2), DIMENSION(ms) :: zmm !node depths in mm - REAL(r_2) :: tempi, temp0,voleq1,wbrat - REAL(r_2), DIMENSION(mp,ms+1) :: ice_correction - + REAL(r_2), dimension(mp,0:ms) :: zimm !layer interface depth in mm + REAL(r_2), dimension(mp,ms) :: zmm !node depths in mm + REAL(r_2) :: tempi, temp0,voleq1,wbrat,& + zi_smpc,tmp_const,voleq2 INTEGER :: k,i - IF (gw_params%ssgw_ice_switch) THEN - smp_cor = 8.0 - ELSE - smp_cor = 0.0 - END IF +!$ IF (gw_params%ssgw_ice_switch) THEN ! FEEDBACK (MMY asks: ask Claire or Anna what ssgw_ice_switch is for and whether to keep it) --rk4417 +!$ smp_cor = 8.0 +!$ ELSE +!$ smp_cor = 0.0 +!$ END IF !make code cleaner define these here - dzmm = 1000.0_r_2 * REAL(soil%zse(:),r_2) - zimm(0) = 0._r_2 - + zimm(:,:) = 0._r_2 DO k=1,ms - zimm(k) = zimm(k-1) + dzmm(k) - zmm(k) = zimm(k-1) + 0.5_r_2*dzmm(k) + zimm(:,k) = zimm(:,k-1) + m2mm*soil%zse_vec(:,k) + zmm(:,k) = zimm(:,k-1) + 0.5_r_2*m2mm*soil%zse_vec(:,k) END DO - DO i=1,mp - GWdzmm(i) = REAL(soil%GWdz(i),r_2)*1000._r_2 - GWzimm(i) = zimm(ms)+GWdzmm(i) - zaq(i) = zimm(ms) + 0.5_r_2*GWdzmm(i) - END DO - - IF (.NOT.gw_params%ssgw_ice_switch) THEN - ice_correction(:,:) = 1._r_2 - - ELSE - - DO k=1,ms - DO i=1,mp - ice_correction(i,k) = 1._r_2 + smp_cor * ssnow%wbice(i,k) - ice_correction(i,k) = ice_correction(i,k)**(2.0/soil%bch_vec(i,k)) - END DO - END DO - DO i=1,mp - ice_correction(i,ms+1) = 1._r_2 + smp_cor * ssnow%wbice(i,ms) - ice_correction(i,ms+1) = ice_correction(i,ms+1)**(2.0/soil%GWbch_vec(i)) - END DO - - END IF - - +!$ IF (.NOT.gw_params%ssgw_ice_switch) THEN ! FEEDBACK (MMY asks: ask Claire or Anna what ssgw_ice_switch is for and whether to keep it) --rk4417 +!$ ice_correction(:,:) = 1._r_2 +!$ +!$ ELSE +!$ +!$ DO k=1,ms +!$ DO i=1,mp +!$ ice_correction(i,k) = 1._r_2 + smp_cor * ssnow%wbice(i,k) +!$ ice_correction(i,k) = ice_correction(i,k)**(2.0/soil%bch_vec(i,k)) +!$ END DO +!$ END DO +!$ DO i=1,mp +!$ ice_correction(i,ms+1) = 1._r_2 + smp_cor * ssnow%wbice(i,ms) +!$ ice_correction(i,ms+1) = ice_correction(i,ms+1)**(2.0/soil%GWbch_vec(i)) +!$ END DO +!$ +!$ END IF + + do i=1,mp + GWzimm(i) = zimm(i,ms)+m2mm*soil%GWdz(i) + zaq(i) = zimm(i,ms) + 0.5_r_2*m2mm*soil%GWdz(i) + end do + + soil%wbc_GW(:) = 0.0 + soil%wbc_vec(:,:) = 0.0 + soil%smpc_vec(:,:) = 1.0e+30 + soil%smpc_GW(:) = 1.0e+30 + + if (.not.gw_params%BC_hysteresis) then !not needed but make code clear + ssnow%watr_hys = soil%watr(:,:) + ssnow%ssat_hys = soil%ssat_vec(:,:) + ssnow%sucs_hys = soil%sucs_vec + + if (gw_params%HC_SWC) then + soil%wbc_vec(:,:) = 2.0*soil%bch_vec(:,:)*soil%ssat_vec(:,:)/& + (1.0+2.0*soil%bch_vec(:,:)) + soil%smpc_vec(:,:) = -soil%sucs_vec(:,:) * & + (2.0*soil%bch_vec(:,:)/& + ((1.0+2.0*soil%bch_vec(:,:))))**(-soil%bch_vec(:,:)) + soil%wbc_GW(:) = 2.0*soil%GWbch_vec(:)*soil%GWssat_vec(:)/& + (1.0+2.0*soil%GWbch_vec(:)) + soil%smpc_GW(:) = -soil%GWsucs_vec(:) * & + (2.0*soil%GWbch_vec(:)/& + ((1.0+2.0*soil%GWbch_vec(:))))**(-soil%GWbch_vec(:)) + end if + end if + !equilibrium water content DO k=1,ms DO i=1,mp - IF ((ssnow%wtd(i) .LE. zimm(k-1))) THEN !fully saturated + if ((ssnow%wtd(i) .le. zimm(i,k-1))) then !fully saturated ! replaced block above by below as per MMY -- rk4417 - ssnow%wbeq(i,k) = soil%ssat_vec(i,k) + ssnow%wbeq(i,k) = ssnow%ssat_hys(i,k) - ELSEIF ((ssnow%wtd(i) .LE. zimm(k)) .AND. & - (ssnow%wtd(i) .GT. zimm(k-1))) THEN + elseif ((ssnow%wtd(i) .le. zimm(i,k)) .and. & + (ssnow%wtd(i) .gt. zimm(i,k-1))) then tempi = 1._r_2 temp0 = & - (((soil%sucs_vec(i,k)+ssnow%wtd(i)-zimm(k-1))/& - soil%sucs_vec(i,k)))**(1._r_2-1._r_2/soil%bch_vec(i,k)) - voleq1 = -soil%sucs_vec(i,k)*(soil%ssat_vec(i,k)-soil%watr(i,k))/& - (1._r_2-1._r_2/soil%bch_vec(i,k))/& - (ssnow%wtd(i)-zimm(k-1))*(tempi-temp0) - ssnow%wbeq(i,k) = (voleq1*(ssnow%wtd(i)-zimm(k-1)) +& - (soil%ssat_vec(i,k)-soil%watr(i,k))& - *(zimm(k)-ssnow%wtd(i)))/(zimm(k)-zimm(k-1))& - + soil%watr(i,k) - ssnow%wbeq(i,k) = ssnow%wbeq(i,k)*ice_correction(i,k) - ELSE - - tempi = (((soil%sucs_vec(i,k)+ssnow%wtd(i)-zimm(k))/& - soil%sucs_vec(i,k)))**(1._r_2-1._r_2/soil%bch_vec(i,k)) - temp0 = (((soil%sucs_vec(i,k)+ssnow%wtd(i)-zimm(k-1))/& - soil%sucs_vec(i,k)))**(1._r_2-1._r_2/soil%bch_vec(i,k)) - ssnow%wbeq(i,k) = -soil%sucs_vec(i,k)*(soil%ssat_vec(i,k)-soil%watr(i,k))/& - (1._r_2-1._r_2/soil%bch_vec(i,k))/& - (zimm(k)-zimm(k-1))*(tempi-temp0)+soil%watr(i,k) - ssnow%wbeq(i,k) = ssnow%wbeq(i,k)*ice_correction(i,k) - END IF - - ssnow%wbeq(i,k) = MIN(MAX(ssnow%wbeq(i,k),soil%watr(i,k)),soil%ssat_vec(i,k)) - - ssnow%wbeq(i,k) = ssnow%wbeq(i,k)*ice_correction(i,k) - - wbrat = MIN(MAX((& - ssnow%wbeq(i,k) - soil%watr(i,k))/(soil%ssat_vec(i,k)-soil%watr(i,k)),& - 0.001_r_2),1._r_2) - - ssnow%zq(i,k) = MAX(& - -soil%sucs_vec(i,k)*(wbrat**(-soil%bch_vec(i,k))),sucmin) - - IF (gw_params%ssgw_ice_switch) THEN - ssnow%zq(i,k) = MAX(& - -soil%sucs_vec(i,k)*(1._r_2+smp_cor*ssnow%wbice(i,k))*(wbrat**(-soil%bch_vec(i,k))),sucmin) - END IF - + (((ssnow%sucs_hys(i,k)+ssnow%wtd(i)-zimm(i,k-1))/& + ssnow%sucs_hys(i,k)))**(1._r_2-1._r_2/soil%bch_vec(i,k)) + voleq1 = -ssnow%sucs_hys(i,k)*(ssnow%ssat_hys(i,k)-ssnow%watr_hys(i,k))/& + (1._r_2-1._r_2/soil%bch_vec(i,k))/& + (ssnow%wtd(i)-zimm(i,k-1))*(tempi-temp0) + ssnow%watr_hys(i,k) ! MMY@23Apr2023 compared with trunk version, Mark added ssnow%watr_hys(i,k) to this equation, but why??? + ssnow%wbeq(i,k) = (voleq1*(ssnow%wtd(i)-zimm(i,k-1)) +& ! MMY@23Apr2023 compared with trunk version, Mark remove ssnow%watr_hys(i,k) from this equation but why??? + (ssnow%ssat_hys(i,k))& + *(zimm(i,k)-ssnow%wtd(i)))/(zimm(i,k)-zimm(i,k-1)) + ! MMY@23Apr2023 compared with trunk version, Mark removed ssnow%wbeq(i,k) = ssnow%wbeq(i,k)*ice_correction(i,k), why??? + else + + tempi = (((ssnow%sucs_hys(i,k)+ssnow%wtd(i)-zimm(i,k))/& + ssnow%sucs_hys(i,k)))**(1._r_2-1._r_2/soil%bch_vec(i,k)) + temp0 = (((ssnow%sucs_hys(i,k)+ssnow%wtd(i)-zimm(i,k-1))/& + ssnow%sucs_hys(i,k)))**(1._r_2-1._r_2/soil%bch_vec(i,k)) + ssnow%wbeq(i,k) = -ssnow%sucs_hys(i,k)*(ssnow%ssat_hys(i,k)-ssnow%watr_hys(i,k))/& + (1._r_2-1._r_2/soil%bch_vec(i,k))/& + (zimm(i,k)-zimm(i,k-1))*(tempi-temp0)+ssnow%watr_hys(i,k) + ! MMY@23Apr2023 compared with trunk version, Mark removed ssnow%wbeq(i,k) = ssnow%wbeq(i,k)*ice_correction(i,k), why??? + end if + + if (.not.gw_params%HC_SWC) then + ssnow%wbeq(i,k) = min(max(ssnow%wbeq(i,k),ssnow%watr_hys(i,k)),ssnow%ssat_hys(i,k)) + + wbrat = min(max((& + ssnow%wbeq(i,k) - ssnow%watr_hys(i,k))/& + (ssnow%ssat_hys(i,k)-ssnow%watr_hys(i,k)),& + 0.01_r_2),1._r_2) + ssnow%zq(i,k) = max(& + -ssnow%sucs_hys(i,k)*(wbrat**(-soil%bch_vec(i,k))),sucmin) + + else + if (ssnow%wbeq(i,k) .lt. soil%wbc_vec(i,k)) then + wbrat = min(max((& + ssnow%wbeq(i,k) - soil%watr(i,k))/(soil%ssat_vec(i,k)-soil%watr(i,k)),& + 0.01_r_2),1._r_2) + + ssnow%zq(i,k) = max(& + -soil%sucs_vec(i,k)*(wbrat**(-soil%bch_vec(i,k))),sucmin) + else + ssnow%zq(i,k) = -soil%sucs_vec(i,k)* & + sqrt(1._r_2 - ssnow%wbeq(i,k)/soil%ssat_vec(i,k))/& + (sqrt(1._r_2-soil%wbc_vec(i,k)/soil%ssat_vec(i,k))*& + ( (soil%wbc_vec(i,k)/soil%ssat_vec(i,k))**soil%bch_vec(i,k) )) + + end if + + end if END DO !mp END DO !ms DO i=1,mp - !Aquifer Equilibrium water content - IF (ssnow%wtd(i) .LE. zimm(ms)) THEN !fully saturated - ssnow%GWwbeq(i) = soil%GWssat_vec(i)-soil%GWwatr(i) + if (ssnow%wtd(i) .le. zimm(i,ms)) then !fully saturated !MMY@23Apr2023 I don't know why Mark changed it but as codes from CABLE-GW works, I suggest to delete the codes from trunk + + ssnow%GWwbeq(i) = soil%GWssat_vec(i) ! MMY@23Apr2023 in trunk it has '-soil%GWwatr(i)' - ELSEIF ((ssnow%wtd(i) .GT. GWzimm(i))) THEN !fully unsaturated + elseif ((ssnow%wtd(i) .gt. GWzimm(i))) then !fully unsaturated tempi = & - (((soil%GWsucs_vec(i)+ssnow%wtd(i)-GWzimm(i))/& - soil%GWsucs_vec(i)))**(1._r_2-1._r_2/soil%GWbch_vec(i)) - temp0 = (((soil%GWsucs_vec(i)+ssnow%wtd(i)-zimm(ms))/& - soil%GWsucs_vec(i)))**(1._r_2-1._r_2/soil%GWbch_vec(i)) + (((soil%GWsucs_vec(i)+ssnow%wtd(i)-GWzimm(i))/& + soil%GWsucs_vec(i)))**(1._r_2-1._r_2/soil%GWbch_vec(i)) + temp0 = (((soil%GWsucs_vec(i)+ssnow%wtd(i)-zimm(i,ms))/& + soil%GWsucs_vec(i)))**(1._r_2-1._r_2/soil%GWbch_vec(i)) ssnow%GWwbeq(i) = -soil%GWsucs_vec(i)*soil%GWssat_vec(i)/& - (1._r_2-1._r_2/soil%GWbch_vec(i))/& - (GWzimm(i)-zimm(ms))*(tempi-temp0) + soil%GWwatr(i) + (1._r_2-1._r_2/soil%GWbch_vec(i))/& + (GWzimm(i)-zimm(i,ms))*(tempi-temp0) + soil%GWwatr(i) - ELSE - - tempi = 1._r_2 - temp0 = (((soil%GWsucs_vec(i)+ssnow%wtd(i)-zimm(ms))/& - soil%GWsucs_vec(i)))**(1._r_2-1._r_2/soil%GWbch_vec(i)) - voleq1 = -soil%GWsucs_vec(i)*(soil%GWssat_vec(i)-soil%GWwatr(i))/& - (1._r_2-1._r_2/soil%GWbch_vec(i))/& - (ssnow%wtd(i)-zimm(ms))*(tempi-temp0) + soil%GWwatr(i) - ssnow%GWwbeq(i) = (voleq1*(ssnow%wtd(i)-zimm(ms)) + & - (soil%GWssat_vec(i)-soil%GWwatr(i))*& - (GWzimm(i)-ssnow%wtd(i)))/(GWzimm(i)-zimm(ms)) + soil%GWwatr(i) - - END IF - - ssnow%GWwbeq(i) = MIN(MAX(ssnow%GWwbeq(i),soil%GWwatr(i)),soil%GWssat_vec(i)) + else - ssnow%GWzq(i) = -soil%GWsucs_vec(i)*(MAX((ssnow%GWwbeq(i)-soil%GWwatr(i))/ & - (soil%GWssat_vec(i)-soil%GWwatr(i)),0.001_r_2))**(-soil%GWbch_vec(i)) - ssnow%GWzq(i) = MAX(sucmin, ssnow%GWzq(i)) + tempi = 1._r_2 + temp0 = & + (((soil%GWsucs_vec(i)+ssnow%wtd(i)-zimm(i,ms))/& + soil%GWsucs_vec(i)))**(1._r_2-1._r_2/soil%GWbch_vec(i)) + voleq1 = -soil%GWsucs_vec(i)*(soil%GWssat_vec(i)-soil%GWwatr(i))/& + (1._r_2-1._r_2/soil%GWbch_vec(i))/& + (ssnow%wtd(i)-zimm(i,ms))*(tempi-temp0) + soil%GWwatr(i) + ssnow%GWwbeq(i) = (voleq1*(ssnow%wtd(i)-zimm(i,ms)) +& ! MMY@23Apr2023 this block doesn't consider 'soil%GWwatr(i))' as the trunk why??? + (soil%GWssat_vec(i))& + *(GWzimm(i)-ssnow%wtd(i)))/(GWzimm(i)-zimm(i,ms)) + + end if + + if (.not.gw_params%HC_SWC) then + + ssnow%GWwbeq(i) = min(max(ssnow%GWwbeq(i),soil%GWwatr(i)),soil%GWssat_vec(i)) + + ssnow%GWzq(i) = -soil%GWsucs_vec(i)*(max((ssnow%GWwbeq(i)-soil%GWwatr(i))/ & + (soil%GWssat_vec(i)-soil%GWwatr(i)),0.01_r_2))**(-soil%GWbch_vec(i)) + ssnow%GWzq(i) = max(sucmin, ssnow%GWzq(i)) + + else + if (ssnow%GWwbeq(i) .lt. soil%wbc_GW(i)) then + wbrat = min(max((& + ssnow%wbeq(i,k) - soil%GWwatr(i))/(soil%GWssat_vec(i)-soil%GWwatr(i)),& + 0.01_r_2),1._r_2) + + ssnow%GWzq(i) = max(& + -soil%GWsucs_vec(i)*(wbrat**(-soil%GWbch_vec(i))),sucmin) + else + ssnow%gwzq(i) = -soil%gwsucs_vec(i)* & + sqrt(1._r_2 - ssnow%gwwbeq(i)/soil%gwssat_vec(i))/& + (sqrt(1._r_2-soil%wbc_gw(i)/soil%gwssat_vec(i))*& + ( (soil%wbc_gw(i)/soil%GWssat_vec(i))**soil%GWbch_vec(i) )) + + end if + end if END DO END SUBROUTINE calc_equilibrium_water_content + SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) + !* Calculate the wet fraction of the surafce + ! following [Decker, 2015](http://doi.wiley.com/10.1002/2015MS000507) + !* This subrountine is called in subrountine surf_wetness_fact in cable_canopy.F90 IMPLICIT NONE TYPE(soil_snow_type), INTENT(INOUT) :: ssnow ! soil+snow variables @@ -1283,12 +1420,14 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) !local variables REAL(r_2), DIMENSION(mp) :: icef,satfrac_liqice,S REAL(r_2) :: fice,xx - REAL(r_2) :: dzmm_one,liqmass,icemass,totmass + REAL(r_2) :: liqmass,icemass,totmass INTEGER :: i,j,k REAL(r_2) :: wb_unsat,wb_lin,funcval REAL(r_2) :: derv,slopeSTDmm,func_step REAL(r_2) :: wb_evap_threshold +! REAL(r_2) :: pi_temp=3.1415927 ! no need, use CPI instead - rk4417 - phase2 + !> if true, using Dani Or soil evaporation scheme IF (cable_user%or_evap) THEN CALL saturated_fraction(ssnow,soil,veg) @@ -1306,19 +1445,20 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) END DO ELSEIF (cable_user%gw_model) THEN - + !> if true, using Dani Or soil evaporation scheme CALL saturated_fraction(ssnow,soil,veg) DO i = 1,mp - dzmm_one = 1000._r_2 * REAL(soil%zse_vec(i,1),r_2) - icemass = ssnow%wbice(i,1) * dzmm_one - liqmass = (ssnow%wb(i,1)-ssnow%wbice(i,1)) * dzmm_one + ! MMY@23Apr2023, calc_srf_wet_fraction is called before ssnow%wmice and ssnow%wmliq calculated, + ! here needs to calcuate icemass and liqmass. + + icemass = ssnow%wbice(i,1)*real(Cdensity_ice,r_2)*soil%zse_vec(i,1) !ice mass + liqmass = ssnow%wbliq(i,1)*real(Cdensity_liq,r_2)*soil%zse_vec(i,1) totmass = MAX(liqmass+icemass,REAL(1e-2,r_2)) - icef(i) = MAX(0._r_2,MIN(1._r_2, gw_params%IceBeta*icemass / totmass)) + icef(i) = max(0._r_2,min(1._r_2,icemass / totmass)) ! MMY END DO - !srf frozen fraction. should be based on topography DO i = 1,mp fice = (EXP(-gw_params%IceAlpha*(1._r_2-icef(i)))-& @@ -1326,9 +1466,9 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) (1._r_2-EXP(-gw_params%IceAlpha)) fice = MIN(1._r_2,MAX(0._r_2,fice)) - satfrac_liqice(i) = fice + (1._r_2-fice)*ssnow%satfrac(i) + satfrac_liqice(i) = max(0.,min(0.95,fice + (1._r_2-fice)*ssnow%satfrac(i) ) ) ! MMY - wb_unsat = ((ssnow%wb(i,1)-ssnow%wbice(i,1)) -& + wb_unsat = ((ssnow%wb(i,1)-ssnow%wbice(i,1)*den_rat) -& ! MMY add *den_rat ! MMY@23Apr2023, since this subroutine is called before den_rat calculated in soil_snow_gw, here den_rat's value from definition at the start of module ssnow%satfrac(i)*soil%ssat_vec(i,1))/(1.-ssnow%satfrac(i)) wb_unsat = MIN(soil%ssat_vec(i,1),MAX(0.,wb_unsat)) @@ -1340,7 +1480,9 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) IF (wb_unsat .GE. wb_evap_threshold) THEN xx = 1. ELSE - xx = 0.25 * (1._r_2 - COS(Cpi*wb_unsat/(wb_evap_threshold)))**2.0 +! xx = 0.25 * (1._r_2 - cos(pi_temp*wb_unsat/(wb_evap_threshold)))**2.0 ! MMY@23Apr2023 add pi_temp to this subroutine, C=>phys which doesn't include pi value +! replaced above by below - rk4417 - phase2 + xx = 0.25 * (1._r_2 - cos(CPI*wb_unsat/(wb_evap_threshold)))**2.0 END IF ssnow%wetfac(i) = MAX(0.0,MIN(1.0,satfrac_liqice(i) +& @@ -1363,14 +1505,15 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) IF( ssnow%wbice(i,1) > 0. )& ssnow%wetfac(i) = ssnow%wetfac(i) * & REAL(MAX( 0.5_r_2, 1._r_2 - MIN( 0.2_r_2, & - ( ssnow%wbice(i,1) / ssnow%wb(i,1) )**2 ) ) ) + ( ssnow%wbice(i,1) * den_rat & ! MMY not sure whether needed but add * den_rat + / ssnow%wb(i,1) )**2 ) ) ) IF( ssnow%snowd(i) > 0.1) ssnow%wetfac(i) = 0.9 - IF ( veg%iveg(i) == 16 .AND. met%tk(i) >= Ctfrz + 5. ) & + IF ( veg%iveg(i) == 16 .AND. met%tk(i) >= CTFRZ + 5. ) & ssnow%wetfac(i) = 1.0 ! lakes: hard-wired number to be removed - IF( veg%iveg(i) == 16 .AND. met%tk(i) < Ctfrz + 5. ) & + IF( veg%iveg(i) == 16 .AND. met%tk(i) < CTFRZ + 5. ) & ssnow%wetfac(i) = 0.7 ! lakes: hard-wired number to be removed ENDDO @@ -1383,137 +1526,9 @@ SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) END SUBROUTINE calc_srf_wet_fraction - !SUBROUTINE calc_soil_hydraulic_props(ssnow,soil,veg) - ! USE cable_common_module - ! TYPE(soil_parameter_type), INTENT(INOUT) :: soil - ! TYPE(soil_snow_type) , INTENT(INOUT) :: ssnow - ! TYPE(veg_parameter_type) , INTENT(INOUT) :: veg - ! - ! INTEGER :: i,k,kk - ! - ! REAL(r_2), DIMENSION(mp) :: s1, & !temporary variables for calculating hydraulic properties - ! s2, & - ! s_mid, & - ! liq_ratio, & - ! Dliq_ratio_Dz - ! - ! REAL(r_2), DIMENSION(0:ms) :: zimm !depths at interface between layers - ! REAL(r_2), dimension(mp,ms) ::wb_temp - ! - ! !soil matric potential, hydraulic conductivity, and - ! derivatives of each with respect to water (calculated using total (not liquid)) - ! - ! do k=1,ms - ! do i=1,mp - ! ssnow%icefrac(i,k) = ssnow%wbice(i,k)/(max(ssnow%wb(i,k),0.01_r_2)) - ! ssnow%fracice(i,k) = (exp(-gw_params%IceAlpha*(1._r_2-ssnow%icefrac(i,k)))& - ! -exp(-gw_params%IceAlpha))/(1._r_2-exp(-gw_params%IceAlpha)) - ! end do - ! end do - ! - ! ssnow%fracice(:,:) = max( min( ssnow%fracice, 1._r_2), 0._r_2) - ! - ! if (gw_params%ssgw_ice_switch) then - ! wb_temp = ssnow%wbliq - ! else - ! wb_temp = ssnow%wb - ! end if - ! - ! do k=1,ms-1 - ! kk=k+1 - ! do i=1,mp - ! - ! s1(i) = 0.5_r_2*(max(wb_temp(i,k)-soil%watr(i,k),0.) + & - ! max(wb_temp(i,kk)-soil%watr(i,kk),0.)) / & - ! (0.5_r_2*((soil%ssat_vec(i,k)-soil%watr(i,k)) + & - ! (soil%ssat_vec(i,kk)-soil%watr(i,kk)))) - ! - ! s1(i) = min(max(s1(i),0.01_r_2),1._r_2) - ! s2(i) = soil%hyds_vec(i,k)*s1(i)**(2._r_2*soil%bch_vec(i,k)+2._r_2) - ! - ! ssnow%hk(i,k) = s1(i)*s2(i) - ! ssnow%dhkdw(i,k) = (2._r_2*soil%bch_vec(i,k)+3._r_2)*s2(i)*& - ! 0.5_r_2/(soil%ssat_vec(i,k)-soil%watr(i,k)) - ! if (.not.gw_params%ssgw_ice_switch) then - ! ssnow%hk(i,k) = ssnow%hk(i,k)*& - ! (1.0 - 0.5_r_2*(ssnow%fracice(i,k)+ssnow%fracice(i,kk))) - ! ssnow%dhkdw(i,k) = ssnow%dhkdw(i,k)*& - ! (1.0 - 0.5_r_2*(ssnow%fracice(i,k)+ssnow%fracice(i,kk))) - ! end if - ! end do - ! end do - ! - ! liq_ratio(:) = 1.0 - ! - ! k = ms - ! do i=1,mp - ! - ! if (gw_params%ssgw_ice_switch) then - ! liq_ratio(i) =min(1.,max(0.,wb_temp(i,k)/max(ssnow%wb(i,k),1e-6) ) ) - ! end if - ! - ! s1(i) = 0.5_r_2*(max(wb_temp(i,k)-soil%watr(i,k),0.) + & - ! max(liq_ratio(i)*ssnow%GWwb(i)-soil%GWwatr(i),0.)) / & - ! (0.5_r_2*(soil%ssat_vec(i,k)-soil%watr(i,k) +& - ! soil%GWssat_vec(i)-soil%GWwatr(i))) - ! - ! s1(i) = min(max(s1(i),0.01_r_2),1._r_2) - ! s2(i) = soil%hyds_vec(i,k)*s1(i)**(2._r_2*soil%bch_vec(i,k)+2._r_2) - ! - ! ssnow%hk(i,k) = s1(i)*s2(i) - ! ssnow%dhkdw(i,k) = (2._r_2*soil%bch_vec(i,k)+3._r_2)*& - ! s2(i)*0.5_r_2/(soil%ssat_vec(i,k)-soil%watr(i,k)) - ! if (.not.gw_params%ssgw_ice_switch) then - ! ssnow%hk(i,k) = (1.-ssnow%fracice(i,k))*ssnow%hk(i,k) - ! ssnow%dhkdw(i,k) = (1.-ssnow%fracice(i,k))*ssnow%dhkdw(i,k) - ! endif - ! end do - ! - ! do k=1,ms - ! do i=1,mp - ! s_mid(i) = (ssnow%wb(i,k)-soil%watr(i,k))/& !+dri*ssnow%wbice(:,k) - ! (soil%ssat_vec(i,k)-soil%watr(i,k)) - ! - ! s_mid(i) = min(max(s_mid(i),0.001_r_2),1._r_2) - ! - ! ssnow%smp(i,k) = -soil%sucs_vec(i,k)*s_mid(i)**(-soil%bch_vec(i,k)) - ! - ! ssnow%smp(i,k) = max(min(ssnow%smp(i,k),-soil%sucs_vec(i,k)),sucmin) - ! - ! ssnow%dsmpdw(i,k) = -soil%bch_vec(i,k)*ssnow%smp(i,k)/& - ! (max(s_mid(i)*(soil%ssat_vec(i,k)-soil%watr(i,k)),0.001_r_2)) - ! end do - ! end do - ! - ! do i=1,mp - ! !Aquifer properties - ! s_mid(i) = (ssnow%GWwb(i)*liq_ratio(i)-soil%GWwatr(i))/& - ! (soil%GWssat_vec(i)-soil%GWwatr(i)) - ! s_mid(i) = min(max(s_mid(i),0.001_r_2),1._r_2) - ! s2(i) = soil%GWhyds_vec(i)*s_mid(i)**(2._r_2*soil%GWbch_vec(i)+2._r_2) - ! - ! ssnow%GWhk(i) =s_mid(i)*s2(i) - ! - ! ssnow%GWdhkdw(i) = (2._r_2*soil%GWbch_vec(i)+3._r_2)*& - ! s2(i)*0.5_r_2/(soil%GWssat_vec(i)-soil%GWwatr(i)) - ! - ! if (.not.gw_params%ssgw_ice_switch) then - ! ssnow%GWhk(i) = (1.-ssnow%fracice(i,ms)) * ssnow%GWhk(i) - ! ssnow%GWdhkdw(i) = (1.-ssnow%fracice(i,ms)) * ssnow%GWdhkdw(i) - ! endif - ! - ! s_mid(i) = (ssnow%GWwb(i)-soil%GWwatr(i))/(soil%GWssat_vec(i)-soil%GWwatr(i)) - ! s_mid(i) = min(max(s_mid(i),0.001_r_2),1._r_2) - ! - ! ssnow%GWsmp(i) = -soil%GWsucs_vec(i)*s_mid(i)**(-soil%GWbch_vec(i)) - ! ssnow%GWsmp(i) = max(min(ssnow%GWsmp(i),-soil%GWsucs_vec(i)),sucmin) - ! ssnow%GWdsmpdw(i) = -soil%GWbch_vec(i)*ssnow%GWsmp(i)/& - ! (s_mid(i)*(soil%GWssat_vec(i)-soil%GWwatr(i))) - ! end do - ! - !END SUBROUTINE calc_soil_hydraulic_props - SUBROUTINE calc_soil_hydraulic_props(ssnow,soil,veg) + !* Calculate soil hydraulic conductivity and dhk/dw + USE cable_common_module TYPE(soil_parameter_type), INTENT(INOUT) :: soil TYPE(soil_snow_type) , INTENT(INOUT) :: ssnow @@ -1532,31 +1547,72 @@ SUBROUTINE calc_soil_hydraulic_props(ssnow,soil,veg) REAL(r_2), DIMENSION(mp,ms+1) :: hk_ice_factor !soil matric potential, hydraulic conductivity, and derivatives of each with respect to water (calculated using total (not liquid)) - - DO k=1,ms - DO i=1,mp - ssnow%icefrac(i,k) = ssnow%wbice(i,k)/(MAX(ssnow%wb(i,k),0.01_r_2)) - ssnow%fracice(i,k) = (EXP(-gw_params%IceAlpha*(1._r_2-ssnow%icefrac(i,k)))& - -EXP(-gw_params%IceAlpha))/(1._r_2-EXP(-gw_params%IceAlpha)) - END DO - END DO + ! __ MMY@23Apr2023 need to test, since itself cannot provide all hysteresis parameters and become redundant ____ + if (gw_params%BC_hysteresis) then + !swc_smp_dsmpdw => brook_corey_hysteresis_swc_smp + ssnow%sucs_hys(:,:) = ssnow%hys_fac(:,:)*soil%sucs_vec(:,:) + elseif (gw_params%HC_SWC) then + !swc_smp_dsmpdw => hutson_cass_swc_smp + ssnow%sucs_hys(:,:) = soil%sucs_vec(:,:) + else + !swc_smp_dsmpdw => brook_corey_swc_smp + ssnow%sucs_hys(:,:) = soil%sucs_vec(:,:) + end if + + ! ___ MMY COMMENT: for Hutson Cass SWC potential ___ + do k=1,ms + do i=1,mp + soil%wbc_vec(i,k) = 2.0*soil%bch_vec(i,k)*soil%ssat_vec(i,k)/& + (1.0+2.0*soil%bch_vec(i,k)) + soil%smpc_vec(i,k) = -soil%sucs_vec(i,k) * & + (2.0*soil%bch_vec(i,k)/& + ((1.0+2.0*soil%bch_vec(i,k))))**(-soil%bch_vec(i,k)) + end do + end do + do i=1,mp + soil%wbc_GW(i) = 2.0*soil%GWbch_vec(i)*soil%GWssat_vec(i)/& + (1.0+2.0*soil%GWbch_vec(i)) + soil%smpc_GW(i) = -soil%GWsucs_vec(i) * & + (2.0*soil%GWbch_vec(i)/& + ((1.0+2.0*soil%GWbch_vec(i))))**(-soil%GWbch_vec(i)) + end do + ! ___________________________________________________ + + do k=1,ms + do i=1,mp + ! _____ MMY it doesn't match to icef in SUBROUTINE ovrlndflx ______ + !ssnow%icefrac(i,k) = ssnow%wbice(i,k)/(max(ssnow%wb(i,k),0.01_r_2)) ! MMY + ! ________________ MMY to unify the calculation ____________________ + ssnow%icefrac(i,k) = max(0._r_2, min(1._r_2, & + ssnow%wmice(i,k) / max(ssnow%wmtot(i,k), 0.01_r_2)) ) + ! __________________________________________________________________ + ssnow%fracice(i,k) = (exp(-gw_params%IceAlpha*(1._r_2-ssnow%icefrac(i,k)))& + -exp(-gw_params%IceAlpha))/(1._r_2-exp(-gw_params%IceAlpha)) + end do + end do ssnow%fracice(:,:) = MAX( MIN( ssnow%fracice, 1._r_2), 0._r_2) IF (gw_params%ssgw_ice_switch) THEN wb_temp(:,1:ms) = ssnow%wbliq(:,:) wb_temp(:,ms+1) = ssnow%GWwb(:) - smp_cor = 8.0 DO k=1,ms kk = MIN(k+1,ms) DO i=1,mp IF (soil%isoilm(i) .EQ. 9) THEN - hk_ice_factor(i,k) = 10.0**(-gw_params%ice_impedence) +!$ hk_ice_factor(i,k) = 10.0**(-gw_params%ice_impedence) ! replaced as per MMY -- rk4417 ! MMY@23Apr2023 I don't know why the sign is changed, need to test its impact ... + !hk_ice_factor(i,k) = (1.0-soil%ssat_vec(i,k))**(gw_params%ice_impedence) + hk_ice_factor(i,k) = 10.0**(gw_params%ice_impedence) + ELSE - hk_ice_factor(i,k) = 10.0**(-gw_params%ice_impedence* & - ( 0.5*(ssnow%wbice(i,k)/MAX(1.0e-8,ssnow%wb(i,k)) + & - ssnow%wbice(i,kk)/MAX(1.0e-8,ssnow%wb(i,kk))) ) & - ) +!$ hk_ice_factor(i,k) = 10.0**(-gw_params%ice_impedence* & ! replaced as per MMY -- rk4417 ! MMY@23Apr2023 keep the commented lines as a reminder for future check +!$ ( 0.5*(ssnow%wbice(i,k)/MAX(1.0e-8,ssnow%wb(i,k)) + & +!$ ssnow%wbice(i,kk)/MAX(1.0e-8,ssnow%wb(i,kk))) ) & +!$ ) + !hk_ice_factor(i,k) = sqrt((1.0-ssnow%wbice(i,k))**(gw_params%ice_impedence) *& + ! (1.0-ssnow%wbice(i,kk))**(gw_params%ice_impedence)) + hk_ice_factor(i,k) = min(10.0**(gw_params%ice_impedence*ssnow%wbice(i,k)/(soil%ssat_vec(i,k)-soil%watr(i,k))) ,& ! MMY@23Apr2023 I don't know why the eq is changed, need to test its impact ... + 10.0**(gw_params%ice_impedence*ssnow%wbice(i,kk)/(soil%ssat_vec(i,kk)-soil%watr(i,kk)))) END IF END DO END DO @@ -1564,7 +1620,6 @@ SUBROUTINE calc_soil_hydraulic_props(ssnow,soil,veg) ELSE wb_temp(:,1:ms) = ssnow%wb(:,:) wb_temp(:,ms+1) = ssnow%GWwb(:) - smp_cor = 0.0 DO k=1,ms kk = MIN(k+1,ms) DO i=1,mp @@ -1579,84 +1634,47 @@ SUBROUTINE calc_soil_hydraulic_props(ssnow,soil,veg) DO i=1,mp liq_ratio(i) =MIN(1.,MAX(0.,wb_temp(i,k)/MAX(ssnow%wb(i,k),1e-6) ) ) END DO - !aquifer ice - wb_temp(:,ms+1) = liq_ratio(:) * wb_temp(:,ms+1) - - !potential from soil water rention function - !defined as layer average - DO k=1,ms - DO i=1,mp - s_mid(i) = (wb_temp(i,k)-soil%watr(i,k))/& !+dri*ssnow%wbice(:,k) - (soil%ssat_vec(i,k)-soil%watr(i,k)) - - s_mid(i) = MIN(MAX(s_mid(i),0.001_r_2),1._r_2) - - ssnow%smp(i,k) = -soil%sucs_vec(i,k)*s_mid(i)**(-soil%bch_vec(i,k))*& - ((1._r_2 + smp_cor*ssnow%wbice(i,k))**2.0) - ssnow%smp(i,k) = MAX(MIN(ssnow%smp(i,k),-soil%sucs_vec(i,k)),sucmin) - - ssnow%dsmpdw(i,k) = -soil%bch_vec(i,k)*ssnow%smp(i,k)/& - (MAX(s_mid(i)*(soil%ssat_vec(i,k)-soil%watr(i,k)),0.001_r_2)) *& - ((1._r_2 + smp_cor*ssnow%wbice(i,k))**2.0) - END DO - END DO - - !Aquifer potential - DO i=1,mp - s_mid(i) = (wb_temp(i,ms+1)-soil%GWwatr(i))/& - (soil%GWssat_vec(i)-soil%GWwatr(i)) - s_mid(i) = MIN(MAX(s_mid(i),0.001_r_2),1._r_2) - s2(i) = soil%GWhyds_vec(i)*s_mid(i)**(2._r_2*soil%GWbch_vec(i)+2._r_2) - - !ssnow%GWhk(i) =s_mid(i)*s2(i) * hk_ice_factor(i,ms+1) - !ssnow%GWdhkdw(i) = (2._r_2*soil%GWbch_vec(i)+3._r_2)*& - ! s2(i)*0.5_r_2/(soil%GWssat_vec(i)-soil%GWwatr(i)) *& - ! hk_ice_factor(i,ms+1) - - ssnow%GWhk(i) = soil%GWhyds_vec(i) * hk_ice_factor(i,ms+1)*& - EXP(-ssnow%wtd(i)/1000._r_2/& - (1.0/(120*(soil%drain_dens(i)+1.0e-3)))) - !d(h)*Sy=dW - ssnow%GWdhkdw(i) = ssnow%GWhk(i)/(soil%GWssat_vec(i)-soil%GWwatr(i))*& - (0.001/(120*(soil%drain_dens(i)+1.0e-3))) + !aquifer ice ! MMY should be aquifer liq + wb_temp(:,ms+1) = liq_ratio(:) * wb_temp(:,ms+1) + ! MMY for ssgw_ice_switch true, wb_temp(:,ms+1) = ssnow%GWwb(:) * wb_temp(:,ms)/max(ssnow%wb(:,ms)) + ! MMY for ssgw_ice_switch false, wb_temp(:,ms+1) = ssnow%GWwb(:) - s_mid(i) = (wb_temp(i,ms+1)-soil%GWwatr(i))/(soil%GWssat_vec(i)-soil%GWwatr(i)) - s_mid(i) = MIN(MAX(s_mid(i),0.001_r_2),1._r_2) + !soil potential (head) calculations can use brooks-corey + !or hutson-cass SWC - ssnow%GWsmp(i) = -soil%GWsucs_vec(i)*s_mid(i)**(-soil%GWbch_vec(i)) - ssnow%GWsmp(i) = MAX(MIN(ssnow%GWsmp(i),-soil%GWsucs_vec(i)),sucmin) - ssnow%GWdsmpdw(i) = -soil%GWbch_vec(i)*ssnow%GWsmp(i)/& - (s_mid(i)*(soil%GWssat_vec(i)-soil%GWwatr(i))) - END DO + !call swc_smp_dsmpdw(soil,ssnow) + if (gw_params%BC_hysteresis) then + call brook_corey_hysteresis_swc_smp(soil,ssnow) + elseif (gw_params%HC_SWC) then + call hutson_cass_swc_smp(soil,ssnow) + else + call brook_corey_swc_smp(soil,ssnow) + end if !hydraulic conductivity !Interfacial so uses layer i and i+1 - DO k=1,ms - kk=MIN(ms+1,k+1) - DO i=1,mp - - IF (k .LT. ms) THEN - s1(i) = 0.5_r_2*((wb_temp(i,k)-soil%watr(i,k)) + & - (wb_temp(i,kk)-soil%watr(i,kk))) / & - (0.5_r_2*((soil%ssat_vec(i,k)-soil%watr(i,k)) + & - (soil%ssat_vec(i,kk)-soil%watr(i,kk)))) - ELSE - s1(i) = 0.5_r_2*((wb_temp(i,k)-soil%watr(i,k)) + & - (wb_temp(i,kk)-soil%GWwatr(i))) / & - (0.5_r_2*((soil%ssat_vec(i,k)-soil%watr(i,k)) + & - (soil%GWssat_vec(i)-soil%GWwatr(i)))) - END IF - s1(i) = MIN(MAX(s1(i),0.01_r_2),1._r_2) - s2(i) = soil%hyds_vec(i,k)*s1(i)**(2._r_2*soil%bch_vec(i,k)+2._r_2) - + do k=1,ms-1 + ! kk=min(ms+1,k+1) ! MMY + kk= k+1 ! MMY + do i=1,mp + ! ____________________ MMY ______________________ + s1(i) = 0.5_r_2*((wb_temp(i,k)-soil%watr(i,k)) + & + (wb_temp(i,kk)-soil%watr(i,kk))) / & + (0.5_r_2*((soil%ssat_vec(i,k)-soil%watr(i,k)) + & + (soil%ssat_vec(i,kk)-soil%watr(i,kk)))) + s1(i) = min(max(s1(i),0.01_r_2),1._r_2) + s2(i) = soil%hyds_vec(i,k)*s1(i)**(2._r_2*soil%bch_vec(i,k)+2._r_2) + ! _________________________________________________ ssnow%hk(i,k) = s1(i)*s2(i)*hk_ice_factor(i,k) ssnow%dhkdw(i,k) = (2._r_2*soil%bch_vec(i,k)+3._r_2)*s2(i)*& - 0.5_r_2/(soil%ssat_vec(i,k)-soil%watr(i,k))*& - hk_ice_factor(i,k) - END DO - END DO + 0.5_r_2/(soil%ssat_vec(i,k)-soil%watr(i,k))*& + hk_ice_factor(i,k) + ! MMY Note that the dhkdw equation doesn't exactly follow the finite difference + end do + end do + !conductivity between soil column and aquifer k = ms kk = ms+1 @@ -1674,13 +1692,23 @@ SUBROUTINE calc_soil_hydraulic_props(ssnow,soil,veg) ssnow%dhkdw(i,k) = (2._r_2*soil%bch_vec(i,k)+3._r_2)*& hk_ice_factor(i,k)*& s2(i)*0.5_r_2/(soil%ssat_vec(i,k)-soil%watr(i,k)) + !Aquifer + + s2(i) = soil%GWhyds_vec(i)*s1(i)**(2._r_2*soil%GWbch_vec(i)+2._r_2) + ssnow%GWhk(i) =s1(i)*s2(i) * hk_ice_factor(i,ms+1) + ssnow%GWdhkdw(i) = (2._r_2*soil%GWbch_vec(i)+3._r_2)*& + s2(i)*0.5_r_2/(soil%GWssat_vec(i)-soil%GWwatr(i)) *& + hk_ice_factor(i,ms+1) END DO END SUBROUTINE calc_soil_hydraulic_props + SUBROUTINE aquifer_recharge(dt,ssnow,soil,veg) + !* Calculate the water movement between the bottom soil layer and + ! groundwater aquifer, following the equation (21) in + ! [Zeng and Decker 2009](https://doi.org/10.1175/2008JHM1011.1) - SUBROUTINE aquifer_recharge(dt,ssnow,soil,veg,zaq,zmm,dzmm) USE cable_common_module IMPLICIT NONE @@ -1688,31 +1716,65 @@ SUBROUTINE aquifer_recharge(dt,ssnow,soil,veg,zaq,zmm,dzmm) TYPE (soil_snow_type), INTENT(INOUT) :: ssnow ! soil and snow variables TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters TYPE (veg_parameter_type), INTENT(INOUT) :: veg - REAL(r_2), DIMENSION(:), INTENT(in) :: zaq - REAL(r_2), DIMENSION(:), INTENT(in) :: zmm,dzmm + REAL(r_2), dimension(mp) :: zaq + REAL(r_2), dimension(mp,ms) :: zmm,csum_dzmm + integer :: i,k + + !> Doing the recharge outside of the soln of Richards Equation makes it easier to track total recharge amount. + + csum_dzmm(:,1) = m2mm*soil%zse_vec(:,1) + do k=2,ms + csum_dzmm(:,k) = csum_dzmm(:,k-1) + m2mm*soil%zse_vec(:,k) + end do - INTEGER :: i !Doing the recharge outside of the soln of Richards Equation makes it easier to track total recharge amount. !Add to ssnow at some point - DO i=1,mp - IF ((ssnow%wtd(i) .LE. SUM(dzmm,dim=1)) .OR. & - (veg%iveg(i) .GE. 16) .OR. & - (soil%isoilm(i) .EQ. 9)) THEN - - ssnow%Qrecharge(i) = 0._r_2 - ELSE - ssnow%Qrecharge(i) = -0.5*(ssnow%hk(i,ms)*ssnow%GWhk(i))*& - ((-ssnow%smp(i,ms)) -& - (-ssnow%zq(i,ms))) / & - (ssnow%wtd(i) - & - (zmm(ms)-0.5*soil%zse_vec(i,ms)*1000.0)) - END IF - END DO + select case (gw_params%aquifer_recharge_function) + !* if water table depth is above the bottom soil coloumn, land cover type is lake + ! or soil type is permanent ice, there is no water movement between bottom soil + ! layer and groundwater aquifer. Otherwise, the water movement will be calcuated + ! as the equation (21) in [Zeng and Decker 2009](https://doi.org/10.1175/2008JHM1011.1) + case(0) + ssnow%Qrecharge(:) = 0._r_2 + case(1) + do i=1,mp + if ((ssnow%wtd(i) .le. csum_dzmm(i,ms-1)) .or. & + (veg%iveg(i) .ge. 16) .or. & + (soil%isoilm(i) .eq. 9)) then + + ssnow%Qrecharge(i) = 0._r_2 + else + ssnow%Qrecharge(i) = -0.5*(ssnow%hk(i,ms)+ssnow%GWhk(i))*& + ((-ssnow%smp(i,ms)) -& + (-ssnow%zq(i,ms))) / & + (ssnow%wtd(i) - & + (csum_dzmm(i,ms)-0.5*soil%zse_vec(i,ms)*m2mm)) + end if + end do + + case default + do i=1,mp + if ((ssnow%wtd(i) .le. csum_dzmm(i,ms)) .or. & + (veg%iveg(i) .ge. 16) .or. & + (soil%isoilm(i) .eq. 9)) then + + ssnow%Qrecharge(i) = 0._r_2 + else + ssnow%Qrecharge(i) = -(ssnow%hk(i,ms)+ssnow%GWhk(i))*& + ((ssnow%GWsmp(i)-ssnow%smp(i,ms)) -& + (ssnow%GWzq(i)-ssnow%zq(i,ms))) / & + (m2mm*(soil%GWdz(i)+soil%zse_vec(i,ms))) + end if + end do + end select END SUBROUTINE aquifer_recharge - SUBROUTINE subsurface_drainage(ssnow,soil,veg,dzmm) + SUBROUTINE subsurface_drainage(ssnow,soil,veg) + !* Calculate horizontal drainage based on the equation (18) in + ! [Decker, 2015](http://doi.wiley.com/10.1002/2015MS000507) + USE cable_common_module IMPLICIT NONE @@ -1720,11 +1782,9 @@ SUBROUTINE subsurface_drainage(ssnow,soil,veg,dzmm) TYPE (soil_snow_type), INTENT(INOUT) :: ssnow ! soil and snow variables TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters TYPE (veg_parameter_type), INTENT(INOUT) :: veg - REAL(r_2), DIMENSION(:), INTENT(in) :: dzmm REAL(r_2), DIMENSION(mp) :: sm_tot,&!sum var - ice_factor_tot !avg ice factor - REAL(r_2), DIMENSION(mp,ms) :: ice_factor !ice limitation on - !subsurface drainage + ice_factor_tot !avg ice factor + REAL(r_2), dimension(mp,ms) :: dzmm, ice_factor !ice limitation on subsurface drainage INTEGER, DIMENSION(mp) :: k_drain INTEGER :: i,k @@ -1733,126 +1793,98 @@ SUBROUTINE subsurface_drainage(ssnow,soil,veg,dzmm) Efold_mod(:) = 1.0 !Efold_mod(1:4) = (/0.2,0.2,0.2,0.2/) !Efold_mod(9) = 0.25 - DO i=1,mp - - !Note: future revision will have interaction with river here. nned to - !work on router and add river type cells - ssnow%qhz(i) = MIN(MAX(soil%slope(i),0.000001),0.9)*& - gw_params%MaxHorzDrainRate* & - EXP(-ssnow%wtd(i)/1000._r_2/& - (1.0/(60.0*gw_params%EfoldHorzDrainRate*& - (soil%drain_dens(i)+1.0e-3))*Efold_mod(veg%iveg(i)))) - - - IF (gw_params%subsurface_sat_drainage) THEN - !drain from sat layers - k_drain(i) = ms+1 - DO k=ms,2,-1 - IF (ssnow%wtd(i) .LE. SUM(dzmm(1:k),dim=1)) THEN - k_drain(i) = k + 1 - END IF - END DO - k_drain(i) = MAX(k_drain(i),3) - ELSE - k_drain(i) = 2 - END IF - - END DO - - IF (gw_params%ssgw_ice_switch) THEN - DO k=1,ms - DO i=1,mp - ice_factor(i,k) = (10.0**(-gw_params%ice_impedence*ssnow%wbice(i,k)/& - (ssnow%wb(i,k)+1.0e-12))) - END DO - END DO - ELSE - DO k=1,ms - DO i=1,mp - ice_factor(i,k) = (1._r_2-ssnow%fracice(i,k)) - END DO - END DO - - END IF - DO i=1,mp - - ice_factor(i,:) = 0._r_2 - ice_factor_tot(i)= 0._r_2 - ssnow%qhlev(i,:) = 0._r_2 - sm_tot(i) = 0._r_2 - - END DO - - DO i=1,mp - IF (gw_params%subsurface_sat_drainage) THEN - sm_tot(i) = MAX((ssnow%GWwb(i) - soil%watr(i,ms)),0._r_2)*& - ice_factor(i,ms) + dzmm(:,:) = m2mm*soil%zse_vec(:,:) - ice_factor_tot(i) = (SUM(ice_factor(i,k_drain(i):ms)*& - soil%zse_vec(i,k_drain(i):ms),dim=1)+& - ice_factor(i,ms)*soil%GWdz(i))/& - (SUM(soil%zse_vec(i,:),dim=1)+& - soil%GWdz(i)) + do i=1,mp - DO k=k_drain(i),ms - sm_tot(i) = sm_tot(i) + MAX(ssnow%wbliq(i,k)-soil%watr(i,k),0._r_2)*& - ice_factor(i,k) - END DO - - ssnow%qhz(i) = ssnow%qhz(i)*ice_factor_tot(i) !reduced due to ice - - IF (sm_tot(i) .GE. 1.0e-12) THEN - DO k=k_drain(i),ms - ssnow%qhlev(i,k) = ssnow%qhz(i)*ice_factor(i,k)*& - MAX(0.0,ssnow%wbliq(i,k)-soil%watr(i,k))/sm_tot(i) - END DO - ssnow%qhlev(i,ms+1) = MAX((ssnow%GWwb(i) - soil%watr(i,ms)),0.0)*& - ice_factor(i,ms)*ssnow%qhz(i)/sm_tot(i) - ENDIF - - ELSE !second option - IF (k_drain(i) .LE. ms) THEN - DO k=k_drain(i),ms - sm_tot(i) = sm_tot(i) + MAX(ssnow%wbliq(i,k)-soil%watr(i,k),0._r_2)*& - ice_factor(i,k) - END DO - ice_factor_tot(i) = (SUM(ice_factor(i,k_drain(i):ms)*& - soil%zse_vec(i,k_drain(i):ms),dim=1))/& - (SUM(soil%zse_vec(i,:),dim=1)) - - ssnow%qhz(i) = ssnow%qhz(i)*ice_factor_tot(i) !reduced due to ice - IF (sm_tot(i) .GE. 1.0e-12) THEN - DO k=k_drain(i),ms - ssnow%qhlev(i,k) = (ssnow%qhz(i)*ice_factor(i,k)/sm_tot(i))*& - MAX(0.0,ssnow%wbliq(i,k)-soil%watr(i,k)) - END DO - ENDIF - - ELSE - ice_factor_tot(i) = ice_factor(i,ms) - ssnow%qhz(i) = ssnow%qhz(i)*ice_factor_tot(i) - ssnow%qhlev(i,ms+1) = ssnow%qhz(i)*MAX(ssnow%GWwb(i)-soil%watr(i,ms),0.0) + !Note: future revision will have interaction with river here. nned to + !work on router and add river type cells - END IF + !* Calulate the total subsurface runoff by following the equation (18) in + ! [Decker, 2015](http://doi.wiley.com/10.1002/2015MS000507) + ssnow%qhz(i) = soil%slope(i)*soil%qhz_max(i)*& + exp( -mm2m*ssnow%wtd(i)/ soil%qhz_efold(i) ) + ! MMY@23Apr2023, defined in derived_parameters: + ! soil%qhz_max(:) = real(gw_params%MaxHorzDrainRate,r_2) + ! soil%qhz_efold(:) = real(gw_params%EfoldHorzDrainScale,r_2)*soil%drain_dens(:)+real(gw_params%EfoldHorzDrainRate,r_2) + ! Note that, this equation is different to the one in the trunk, I have no idea about why it is changed + + !* calculate how many layers are below water table depth and can generate subsurface runoff + !drain from sat layers + k_drain(i) = ms+1 + do k=ms,2,-1 + if (ssnow%wtd(i) .le. sum(dzmm(i,1:k),dim=1)) then + k_drain(i) = k + 1 + end if + end do + k_drain(i) = max(k_drain(i),3) + + end do + + if (gw_params%ssgw_ice_switch) then + do k=1,ms + do i=1,mp + ice_factor(i,k) = (1.0-ssnow%wbice(i,k))**gw_params%ice_impedence + ! MMY@23Apr2023 in trunk: + ! ice_factor(i,k) = (10.0**(-gw_params%ice_impedence*ssnow%wbice(i,k)/(ssnow%wb(i,k)+1.0e-12))) + end do + end do + else + do k=1,ms + do i=1,mp + ice_factor(i,k) = (1._r_2-ssnow%fracice(i,k)) + end do + end do + + end if + + ssnow%qhlev(:,:) = 0._r_2 + sm_tot(:) = 0._r_2 + + do i=1,mp + + if (gw_params%subsurface_sat_drainage) then ! MMY@23Apr2023 I think subsurface_sat_drainage should be the default setting + ! MMY it doesn't make sense since sm_tot == 0. + ! sm_tot(i) = sum(ssnow%qhlev(i,k_drain(i):ms+1),dim=1) ! MMY + ! ______ MMY find this calculation from version : cable-2.2.3-pore-scale-model _______________ + sm_tot(i) = max(ssnow%GWwb(i)-soil%GWwatr(i),0._r_2) + do k=k_drain(i),ms + sm_tot(i) = sm_tot(i) + max(ssnow%wbliq(i,k)-soil%watr(i,k),0._r_2)!*dzmm(k) + end do + ! ____________________________________________________________________________________________ + if (sm_tot(i) .lt. 1.0e-8) then + sm_tot(i) = 1._r_2 + end if + end if + end do + + + do i=1,mp + + ssnow%qhlev(i,ms+1) = max((ssnow%GWwb(i) - soil%GWwatr(i)),0._r_2)*& + ice_factor(i,ms)*ssnow%qhz(i)/sm_tot(i) + + do k=k_drain(i),ms + ssnow%qhlev(i,k) = max(ssnow%wbliq(i,k)-ssnow%watr_hys(i,k),0._r_2)*& + ice_factor(i,k)*ssnow%qhz(i)/sm_tot(i) + end do - END IF !incase every layer is frozen very dry - ssnow%qhz(i) = SUM(ssnow%qhlev(i,:),dim=1) + ssnow%qhz(i) = sum(ssnow%qhlev(i,:),dim=1) !Keep "lakes" saturated forcing qhz = 0. runoff only from lakes !overflowing - IF (soil%isoilm(i) .EQ. 9 .OR. veg%iveg(i) .GE. 16) THEN + if (soil%isoilm(i) .eq. 9 .or. veg%iveg(i) .ge. 16) then ssnow%qhz(i) = 0._r_2 ssnow%qhlev(i,:) = 0._r_2 - END IF + end if - END DO + end do END SUBROUTINE subsurface_drainage - SUBROUTINE saturated_fraction(ssnow,soil,veg) TYPE (soil_snow_type), INTENT(INOUT) :: ssnow ! soil and snow variables TYPE (soil_parameter_type), INTENT(IN) :: soil ! soil parameters @@ -1875,8 +1907,8 @@ SUBROUTINE saturated_fraction(ssnow,soil,veg) S(:) = 0._r_2 DO k=1,gw_params%level_for_satfrac S(:) = S(:) + MAX(0.01,MIN(1.0, & - (ssnow%wb(:,k)-ssnow%wbice(:,k)-soil%watr(:,k))/& - MAX(0.001,soil%ssat_vec(:,k)-soil%watr(:,k)) ) )*soil%zse_vec(:,k) + (ssnow%wb(:,k)-den_rat*ssnow%wbice(:,k)-ssnow%watr_hys(:,k))/& + max(0.001,ssnow%ssat_hys(:,k)-soil%watr(:,k)-den_rat*ssnow%wbice(:,k)) ) )*soil%zse_vec(:,k) END DO S(:) = S(:)/SUM(soil%zse(1:gw_params%level_for_satfrac),dim=1) !srf frozen fraction. should be based on topography @@ -1888,7 +1920,7 @@ SUBROUTINE saturated_fraction(ssnow,soil,veg) 1e-5),10000._r_2)) ! ensure some variability ssnow%satfrac(i) = MAX(0._r_2,MIN(0.99_r_2,& !note UM wants std03, and erf is not included then - 1._r_2 - my_erf( slopeSTDmm / SQRT(2.0* S(i)) ) ) ) + 1._r_2 - my_erf( slopeSTDmm * sqrt(2.0* S(i)) ) ) ) ! MMY change / to *@24 Sep 2021 ELSEIF (veg%iveg(i) .LT. 16) THEN ssnow%satfrac(i) = 0._r_2 ELSE @@ -1906,9 +1938,12 @@ SUBROUTINE pore_space_relative_humidity(ssnow,soil,veg) REAL(r_2), DIMENSION(mp) :: unsat_wb,unsat_smp - ! Need a matching array of NEGATIVE ones for the intrinsic sign func below - REAL(r_2), DIMENSION(mp,ms) :: minus_ones - + ! Need a matching array of ones to use in Mark's call to the intrinsic + ! sign func below + ! mgk, 24/07/2018 +!$ REAL(r_2), DIMENSION(mp,ms) :: ones + REAL(r_2), DIMENSION(mp,ms) :: minus_ones ! rk4417 as per MMY + INTEGER :: i !if gw_model = true @@ -1918,24 +1953,26 @@ SUBROUTINE pore_space_relative_humidity(ssnow,soil,veg) !cable_um_init_subrs.F90 or cable_parameters: !ssat_vec(i,:) = ssat(i) !so ssat_vec can be used although soilsnow uses ssat - - minus_ones = -1.0_r_2 + + ! mgk, 24/07/2018 - fix to compile + minus_ones = -1. DO i=1,mp IF (veg%iveg(i) .LT. 16 .AND. soil%isoilm(i) .NE. 9 .AND. & ssnow%snowd(i) .LE. 1e-8 ) THEN - unsat_wb(i) = (ssnow%wb(i,1) - soil%ssat_vec(i,1)*& + unsat_wb(i) = (ssnow%wb(i,1) - ssnow%ssat_hys(i,1)*& MIN(0.95,MAX(0.0,ssnow%satfrac(i))))/(1.0 - MIN(0.95,MAX(0.0,ssnow%satfrac(i)))) - - unsat_wb(i) = MAX(soil%watr(i,1)+1e-2, MIN(soil%ssat_vec(i,1), unsat_wb(i) ) ) - + ! MMY ??? I'm not sure whether we should take wbice*den_rat into consideration here, + ! ??? so I didn't modify. + unsat_wb(i) = max(ssnow%watr_hys(i,1)+1e-2, min(ssnow%ssat_hys(i,1), unsat_wb(i) ) ) + !unsat_smp(i) = sign(soil%sucs_vec(i,1),-1.0) * min(1.0, & ! (max(0.001, (unsat_wb(i)-soil%watr(i,1))/(soil%ssat_vec(i,1)-& ! soil%watr(i,1)) ) )** (-soil%bch_vec(i,1)) ) ! mgk, 24/07/2018 - fix to compile - unsat_smp(i) = SIGN(soil%sucs_vec(i,1),minus_ones(i,1)) * MIN(1.0, & + unsat_smp(i) = SIGN(soil%sucs_vec(i,1),minus_ones(i,1)) * MIN(1.0, & (MAX(0.001, (unsat_wb(i)-soil%watr(i,1))/(soil%ssat_vec(i,1)-& soil%watr(i,1)) ) )** (-soil%bch_vec(i,1)) ) @@ -1963,9 +2000,9 @@ SUBROUTINE sli_hydrology(dels,ssnow,soil,veg,canopy) LOGICAL, SAVE :: sli_call = .TRUE. - REAL(r_2), DIMENSION(ms) :: dzmm - REAL(r_2), DIMENSION(mp) :: zmm - REAL(r_2), DIMENSION(mp) :: zaq + ssnow%sucs_hys = soil%sucs_vec + ssnow%ssat_hys = soil%ssat_vec + ssnow%watr_hys = soil%watr CALL iterative_wtd (ssnow, soil, veg, cable_user%test_new_gw) @@ -1973,21 +2010,18 @@ SUBROUTINE sli_hydrology(dels,ssnow,soil,veg,canopy) CALL ovrlndflx (dels, ssnow, soil, veg,canopy,sli_call ) - dzmm = REAL(soil%zse(:),r_2)*1000._r_2 - - CALL subsurface_drainage(ssnow,soil,veg,dzmm) - - zmm(:) = 1000._r_2*(SUM(REAL(soil%zse,r_2),dim=1)) - zaq(:) = zmm(:) + 0.5_r_2*soil%GWdz(:)*1000._r_2 - - CALL aquifer_recharge(dels,ssnow,soil,veg,zaq,zmm,zmm) - - + CALL subsurface_drainage(ssnow,soil,veg) + call aquifer_recharge(dels,ssnow,soil,veg) END SUBROUTINE sli_hydrology - +! MMY@23Apr2023 Please keep SUBROUTINE set_unsed_gw_vars to see whether can call it in future codes to normalize +! CABLE writing. Open to any suggestion. +!----------------------------------------------------------------------------------------- +! subroutine set_unsed_gw_vars checked by rk4417. This subroutine is not called anywhere. +!----------------------------------------------------------------------------------------- + SUBROUTINE set_unsed_gw_vars(ssnow,soil,canopy) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow ! soil+snow variables TYPE(soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters @@ -2031,4 +2065,499 @@ REAL(r_2) FUNCTION my_erf(x) END FUNCTION my_erf + subroutine brook_corey_swc_smp(soil,ssnow) + !* Using Brook_Corey equation to calculate soil matric potential + + type(soil_parameter_type), intent(inout) :: soil + type(soil_snow_type), intent(inout) :: ssnow + + real(r_2), dimension(mp,ms+1) :: wb_temp + real(r_2), dimension(mp) :: s_mid + integer :: i,k + + if (gw_params%ssgw_ice_switch) then + wb_temp(:,1:ms) = ssnow%wbliq(:,:) + wb_temp(:,ms+1) = ssnow%GWwb(:) + do i=1,mp + if (ssnow%wb(i,ms) .gt. 1.0e-8) then + wb_temp(i,ms+1) = ssnow%GWwb(i) * ssnow%wbliq(i,ms)/ssnow%wb(i,ms) + end if + end do + else + wb_temp(:,1:ms) = ssnow%wb(:,:) + wb_temp(:,ms+1) = ssnow%GWwb(:) + end if + !potential from soil water rention function + !defined as layer average + do k=1,ms + do i=1,mp + s_mid(i) = (wb_temp(i,k)-soil%watr(i,k))/& !+dri*ssnow%wbice(:,k) + (soil%ssat_vec(i,k)-soil%watr(i,k)) + + s_mid(i) = min(max(s_mid(i),0.01_r_2),1._r_2) + + ssnow%smp(i,k) = -soil%sucs_vec(i,k)*s_mid(i)**(-soil%bch_vec(i,k)) + + ssnow%smp(i,k) = max(min(ssnow%smp(i,k),-soil%sucs_vec(i,k)),sucmin) + + ssnow%dsmpdw(i,k) = -soil%bch_vec(i,k)*ssnow%smp(i,k)/& + (max(s_mid(i)*(wb_temp(i,k)-soil%watr(i,k)),0.01_r_2)) !MMY + ! MMY BUG in (max(s_mid(i)*(soil%ssat_vec(i,k)-soil%watr(i,k)),0.01_r_2)) + end do + end do + ! _____________________ MMY BUG it forgot to calculate the aquifer __________________ + !Aquifer potential + do i=1,mp + s_mid(i) = (wb_temp(i,ms+1)-soil%GWwatr(i))/(soil%GWssat_vec(i)-soil%GWwatr(i)) + s_mid(i) = min(max(s_mid(i),0.001_r_2),1._r_2) + + ssnow%GWsmp(i) = -soil%GWsucs_vec(i)*s_mid(i)**(-soil%GWbch_vec(i)) + ssnow%GWsmp(i) = max(min(ssnow%GWsmp(i),-soil%GWsucs_vec(i)),sucmin) + ssnow%GWdsmpdw(i) = -soil%GWbch_vec(i)*ssnow%GWsmp(i)/& + (max(s_mid(i)*(wb_temp(i,ms+1)-soil%GWwatr(i)),0.01_r_2)) !MMY + end do + ! ____________________________________________________________________________________ + end subroutine brook_corey_swc_smp + + ! MMY@13Mar2023 Please keep subroutine hutson_cass_swc_smp, test this subroutine when there is no problem + ! in the merged version. + subroutine hutson_cass_swc_smp(soil,ssnow) + + !* Using Hutson-CASS equation to calculate soil matric potential + + type(soil_parameter_type), intent(inout) :: soil + type(soil_snow_type), intent(inout) :: ssnow + + real(r_2), dimension(mp,ms+1) :: wb_temp + real(r_2), dimension(mp) :: s_mid + integer :: i,k + + do k=1,ms + do i=1,mp + soil%wbc_vec(i,k) = 2.0*soil%bch_vec(i,k)*soil%ssat_vec(i,k)/& + (1.0+2.0*soil%bch_vec(i,k)) + soil%smpc_vec(i,k) = -soil%sucs_vec(i,k) * & + (2.0*soil%bch_vec(i,k)/& + ((1.0+2.0*soil%bch_vec(i,k))))**(-soil%bch_vec(i,k)) + end do + end do + + do i=1,mp + soil%wbc_GW(i) = 2.0*soil%GWbch_vec(i)*soil%GWssat_vec(i)/& + (1.0+2.0*soil%GWbch_vec(i)) + soil%smpc_GW(i) = -soil%GWsucs_vec(i) * & + (2.0*soil%GWbch_vec(i)/& + ((1.0+2.0*soil%GWbch_vec(i))))**(-soil%GWbch_vec(i)) + end do + !potential from soil water rention function + !defined as layer average + do k=1,ms + do i=1,mp + s_mid(i) = (wb_temp(i,k)-soil%watr(i,k))/& !+dri*ssnow%wbice(:,k) + (soil%ssat_vec(i,k)-soil%watr(i,k)) + + s_mid(i) = min(max(s_mid(i),0.01_r_2),0.9999999_r_2) + + if (wb_temp(i,k) .lt. soil%wbc_vec(i,k)) then + ssnow%smp(i,k) = -soil%sucs_vec(i,k)*s_mid(i)**(-soil%bch_vec(i,k)) + + ssnow%smp(i,k) = max(min(ssnow%smp(i,k),-soil%sucs_vec(i,k)),sucmin) + + ssnow%dsmpdw(i,k) = -soil%bch_vec(i,k)*ssnow%smp(i,k)/& + (max(s_mid(i)*(wb_temp(i,k)-soil%watr(i,k)),0.001_r_2)) ! MMY + ! MMY BUG in (max(s_mid(i)*(soil%ssat_vec(i,k)-soil%watr(i,k)),0.001_r_2)) + else + !!! MMY BUG: in Hutson & Cass 1987 using 'wb/ssat' as independent factor + !!! MMY but here Mark uses effective saturation (wb-watr)/(ssat-watr) + !!! MMY in wet soil and wb/ssat in wet soil. This leads to a iscontinuity + !!! MMY in the water retention curve (smp-wb relationship), and causes + !!! MMY abs(smp) in wet soil is larger than in dry soil. I suggest to + !!! MMY use wb/ssat in HC_SWC + ssnow%smp(i,k) = -soil%sucs_vec(i,k)* sqrt(1._r_2 - wb_temp(i,k)/soil%ssat_vec(i,k))/& + (sqrt(1._r_2-soil%wbc_vec(i,k)/soil%ssat_vec(i,k))*& + ( (soil%wbc_vec(i,k)/soil%ssat_vec(i,k))**soil%bch_vec(i,k) )) + + ssnow%dsmpdw(i,k) = -ssnow%smp(i,k)/(soil%ssat_vec(i,k)*& + (1._r_2 - wb_temp(i,k)/soil%ssat_vec(i,k)) ) + end if + end do + end do + + !Aquifer potential + do i=1,mp + s_mid(i) = (wb_temp(i,ms+1)-soil%GWwatr(i))/(soil%GWssat_vec(i)-soil%GWwatr(i)) + s_mid(i) = min(max(s_mid(i),0.01_r_2),1._r_2) + if (wb_temp(i,ms+1) .lt. soil%wbc_GW(i)) then + ssnow%GWsmp(i) = -soil%GWsucs_vec(i)*s_mid(i)**(-soil%GWbch_vec(i)) + + ssnow%GWsmp(i) = max(min(ssnow%GWsmp(i),-soil%GWsucs_vec(i)),sucmin) + + ssnow%GWdsmpdw(i) = -soil%GWbch_vec(i)*ssnow%GWsmp(i)/& + (max(s_mid(i)*(soil%GWssat_vec(i)-soil%GWwatr(i)),0.001_r_2)) + else + ssnow%GWsmp(i) = -soil%GWsucs_vec(i)* sqrt(1._r_2 - wb_temp(i,ms+1)/soil%GWssat_vec(i))/& + (sqrt(1._r_2-soil%wbc_GW(i)/soil%GWssat_vec(i))*& + ( (soil%wbc_GW(i)/soil%GWssat_vec(i))**soil%GWbch_vec(i) )) + + ssnow%GWsmp(i) = max(min(ssnow%GWsmp(i),-soil%GWsucs_vec(i)),sucmin) + + ssnow%GWdsmpdw(i) = -ssnow%GWsmp(i)/(soil%GWssat_vec(i)*& + (1._r_2 - wb_temp(i,ms+1)/soil%GWssat_vec(i)) ) + end if + end do + end subroutine hutson_cass_swc_smp + + ! MMY@13Mar2023 Please keep subroutine brook_corey_hysteresis_swc_smp temporarily. + ! Hysteresis in water retention curve isn't commonly considered in + ! land surface models. I doubt whether this subrountine is necessary + ! for CABLE. But keep it at this moment, may delete it if this subroutine + ! is actually broken. + + subroutine brook_corey_hysteresis_swc_smp(soil,ssnow) + + !* Using Brook_Corey equation to calculate soil matric potential + ! but considering the hysteresis of water retention curve + + type(soil_parameter_type), intent(inout) :: soil + type(soil_snow_type), intent(inout) :: ssnow + + real(r_2), dimension(mp,ms) :: delta_wbliq + real(r_2), dimension(mp) :: s_mid + integer :: i,k,klev + real(r_2), dimension(mp,ms+1) :: wb_temp + real(r_2) :: tmp_smp + + if (gw_params%ssgw_ice_switch) then + wb_temp(:,1:ms) = ssnow%wbliq(:,:) + wb_temp(:,ms+1) = ssnow%GWwb(:) * ssnow%wbliq(:,ms)/soil%ssat_vec(:,ms) + else + wb_temp(:,1:ms) = ssnow%wb(:,:) + wb_temp(:,ms+1) = ssnow%GWwb(:) + end if + + !ensure sucs_hys is set + ssnow%sucs_hys = ssnow%hys_fac * soil%sucs_vec +! do k=1,ms +! do i=1,mp +! !on dry or wet? +! IF (ssnow%hys_fac(i,k) .gt. 0.75) then +! !on drying +! ssnow%watr_hys(i,k) = soil%watr(i,k) +! ssnow%hys_fac(i,k) = 1.0 !drying +! ssnow%sucs_hys(i,k) = ssnow%hys_fac(i,k)*soil%sucs_vec(i,k) +! if (ssnow%smp_hys(i,k) .lt. -ssnow%sucs_hys(i,k)) then +! tmp_smp = (min(0.99,-ssnow%smp_hys(i,k)/ssnow%sucs_hys(i,k)))**(-soil%bch_vec(i,k)) +! ssnow%ssat_hys(i,k) = (ssnow%wb_hys(i,k) - ssnow%watr_hys(i,k)*(1.0-tmp_smp))/& +! tmp_smp +! else +! ssnow%ssat_hys(i,k) = ssnow%wb_hys(i,k) +! end if +! ssnow%ssat_hys(i,k) = max(0.001+ssnow%watr_hys(i,k),ssnow%ssat_hys(i,k)) +! +! ELSE !wetting +! ssnow%ssat_hys(i,k) = gw_params%ssat_wet_factor*soil%ssat_vec(i,k) +! ssnow%hys_fac(i,k) = 0.5 !wetting +! ssnow%sucs_hys(i,k) = ssnow%hys_fac(i,k)*soil%sucs_vec(i,k) +! +! if (ssnow%smp_hys(i,k) .lt. -ssnow%sucs_hys(i,k) ) then +! tmp_smp = (min(0.99,-ssnow%smp_hys(i,k)/ssnow%sucs_hys(i,k)))**(-soil%bch_vec(i,k)) +! ssnow%watr_hys(i,k) = (ssnow%wb_hys(i,k) - ssnow%ssat_hys(i,k) *tmp_smp)/& +! (1._r_2 - tmp_smp) +! else +! ssnow%watr_hys(i,k) = soil%watr(i,k) +! end if +! +! ssnow%watr_hys(i,k) = max(0._r_2,min(ssnow%ssat_hys(i,k)-0.001, ssnow%watr_hys(i,k) ) ) +! END IF +! END DO +! END DO + + do k=1,ms + do i=1,mp + s_mid(i) = (wb_temp(i,k) - ssnow%watr_hys(i,k)) / & + (ssnow%ssat_hys(i,k) - ssnow%watr_hys(i,k)) + + s_mid(i) = min(max(s_mid(i),0.01_r_2),1._r_2) + ! MMY ??? why 0.01_r_2 here but 0.001_r_2 in aquifer + ssnow%smp(i,k) = -ssnow%sucs_hys(i,k)*s_mid(i)**(-soil%bch_vec(i,k)) + + ssnow%smp(i,k) = max(min(ssnow%smp(i,k),-ssnow%sucs_hys(i,k)),sucmin) + + ssnow%dsmpdw(i,k) = -soil%bch_vec(i,k)*ssnow%smp(i,k)/s_mid(i) + end do + end do + + !Aquifer potential + do i=1,mp + s_mid(i) = (wb_temp(i,ms+1)-soil%GWwatr(i))/(soil%GWssat_vec(i)-soil%GWwatr(i)) + s_mid(i) = min(max(s_mid(i),0.001_r_2),1._r_2) + ! MMY ??? + ssnow%GWsmp(i) = -soil%GWsucs_vec(i)*s_mid(i)**(-soil%GWbch_vec(i)) + ssnow%GWsmp(i) = max(min(ssnow%GWsmp(i),-soil%GWsucs_vec(i)),sucmin) + ssnow%GWdsmpdw(i) = -soil%GWbch_vec(i)*ssnow%GWsmp(i)/& + (max(s_mid(i)*(wb_temp(i,ms+1)-soil%GWwatr(i)),0.01_r_2)) !MMY + end do + + end subroutine brook_corey_hysteresis_swc_smp + + ! MMY@13Mar2023 Please keep subroutine swc_hyst_direction temporarily. Subroutine + ! swc_hyst_direction should be kept or deleted as subroutine + ! brook_corey_hysteresis_swc_smp + +subroutine swc_hyst_direction(soil,ssnow,veg) + !* Calculate soil water content and soil hydraulic parameters + ! with the consideration of hysteresis effect + + type(soil_parameter_type), intent(inout) :: soil + type(soil_snow_type), intent(inout) :: ssnow + TYPE(veg_parameter_type) , INTENT(INOUT) :: veg ! veg parameters + + real(r_2), dimension(mp,ms) :: delta_wbliq,psi_tmp + integer :: i,k,klev + real(r_2) :: tmp_smp + + delta_wbliq = ssnow%wbliq - ssnow%wbliq_old + !soil hydraulic state/props so smp out matches the wb out + !call swc_smp_dsmpdw(soil,ssnow) + + if (gw_params%BC_hysteresis) then ! MMY??? redundant,calling swc_hyst_direction requires gw_params%BC_hysteresis, + call brook_corey_hysteresis_swc_smp(soil,ssnow) + elseif (gw_params%HC_SWC) then + call hutson_cass_swc_smp(soil,ssnow) + else + call brook_corey_swc_smp(soil,ssnow) + end if + + !switch drying/wetting curve + do k=1,ms + do i=1,mp + if (delta_wbliq(i,k) .gt. 0.0 .and. & + ssnow%hys_fac(i,k) .gt. 0.75.and.& + ssnow%wb(i,k).lt.gw_params%ssat_wet_factor * soil%ssat_vec(i,k)) then !avoid testing .eq. 1.0 + !layer was drying now it is wetting! + ssnow%smp_hys(i,k) = ssnow%smp(i,k) + ssnow%wb_hys(i,k) = ssnow%wb(i,k) + ssnow%ssat_hys(i,k)= gw_params%ssat_wet_factor * soil%ssat_vec(i,k) + ssnow%hys_fac(i,k) = 0.5 + ssnow%sucs_hys(i,k) = ssnow%hys_fac(i,k)*soil%sucs_vec(i,k) + elseif (delta_wbliq(i,k) .le. 0.0 .and. & + ssnow%hys_fac(i,k) .lt. 0.75 .and.& + ssnow%wb(i,k) .gt. soil%watr(i,k)) then + !swtiched wetting to drying + ssnow%smp_hys(i,k) = ssnow%smp(i,k) + ssnow%wb_hys(i,k) = ssnow%wb(i,k) + ssnow%hys_fac(i,k) = 1.0 + ssnow%sucs_hys(i,k) = ssnow%hys_fac(i,k)*soil%sucs_vec(i,k) + end if + + if (delta_wbliq(i,k) .gt. 0.0) then + + ssnow%ssat_hys(i,k)= gw_params%ssat_wet_factor * soil%ssat_vec(i,k) + if (ssnow%smp_hys(i,k) .lt. -ssnow%sucs_hys(i,k) ) then + + tmp_smp = (min(0.99,-ssnow%smp_hys(i,k)/ssnow%sucs_hys(i,k)))**(-soil%bch_vec(i,k)) + ssnow%watr_hys(i,k) = (ssnow%wb_hys(i,k) - soil%ssat_vec(i,k) * tmp_smp)/& + (1.0 - tmp_smp) + else + ssnow%watr_hys(i,k) = ssnow%wb_hys(i,k) + end if + + ssnow%watr_hys(i,k) = max(0._r_2,min(ssnow%ssat_hys(i,k)-0.001, ssnow%watr_hys(i,k) ) ) + + else !drying ! MMY add a space + ssnow%watr_hys(i,k)= soil%watr(i,k) !this is how we tell if drying + !watr_hys .eq. watr + if (ssnow%smp_hys(i,k) .lt. -ssnow%sucs_hys(i,k)) then + tmp_smp = (min(0.99,-ssnow%smp_hys(i,k)/ssnow%sucs_hys(i,k)))**(-soil%bch_vec(i,k)) + ssnow%ssat_hys(i,k) = (ssnow%wb_hys(i,k) - soil%watr(i,k)*(1.0-tmp_smp))/& + tmp_smp + else + ssnow%ssat_hys(i,k) = ssnow%wb_hys(i,k) + end if + + ssnow%ssat_hys(i,k) = max(ssnow%watr_hys(i,k)+0.001,ssnow%ssat_hys(i,k)) + end if + + end do + end do + + if (cable_user%gw_model .and. gw_params%bc_hysteresis) then + do klev=1,ms + do i=1,mp + if (soil%isoilm(i) .ne. 9 .and. veg%iveg(i) .le. 16) then + + psi_tmp(i,klev) = abs(psi_c(veg%iveg(i))) + + soil%swilt_vec(i,klev) = (ssnow%ssat_hys(i,klev)-ssnow%watr_hys(i,klev)) * & + (psi_tmp(i,klev)/soil%sucs_vec(i,klev))& + **(-1.0/soil%bch_vec(i,klev))+& + ssnow%watr_hys(i,klev) + soil%sfc_vec(i,klev) = (gw_params%sfc_vec_hk/soil%hyds_vec(i,klev))& + **(1.0/(2.0*soil%bch_vec(i,klev)+3.0)) *& + (ssnow%ssat_hys(i,klev)-ssnow%watr_hys(i,klev)) + ssnow%watr_hys(i,klev) + end if + end do + end do + + end if + + end subroutine swc_hyst_direction + +!`!taken from +!`!/g/data1/w35/mrd561/CABLE/CMIP6-GM2_dev_testing_tiles/core/biogeophys +!`! to test +!`SUBROUTINE GWsoilfreeze(dels, soil, ssnow) +!` !NOTE: this is only included because gw_model uses parameters XXX_vec +!` !these are r_2. this breaks bitwise compatibility with trunk +!` !if acceptable this routine does the same thing but with r_2 soil params +!` ! if max_ice_frac always set to frozen_limit and tgg_tmp is always CTFRZ +!` +!` REAL, INTENT(IN) :: dels ! integration time step (s) +!` TYPE(soil_snow_type), INTENT(INOUT) :: ssnow +!` TYPE(soil_parameter_type), INTENT(INOUT) :: soil +!` REAL , DIMENSION(mp,ms) :: tgg_old,tgg_new,tgg_tmp !tgg_old is previous point of when crosses freezing +!` REAL(r_2), DIMENSION(mp,ms) :: sicefreeze +!` REAL(r_2), DIMENSION(mp,ms) :: sicemelt +!` REAL(r_2), DIMENSION(mp,ms) :: wbice_delta,avail_por,delta_ice_vol +!` REAL(r_2), DIMENSION(mp) :: ice_mass,liq_mass,tot_mass +!` INTEGER :: i,j,k +!` REAL(r_2) :: func,funcderv,Aconst,Dconst,t_zero,t_one,dtmp +!` REAL, DIMENSION(mp,ms) :: gammzz_snow +!` REAL(r_2),DIMENSION(mp,ms) :: xx,max_ice_frac,iceF,den_css !Decker and Zeng 20.9 +!` REAL(r_2) :: delta_wbliq,tmp_var +!` +!` max_ice_frac(:,:) = zero +!` delta_ice_vol(:,:) = zero +!` tgg_old(:,:) = ssnow%otgg(:,:) +!` tgg_new(:,:) = ssnow%tgg(:,:) +!` tgg_tmp(:,:) = tgg_old(:,:) +!` +!` gammzz_snow(:,:) = zero +!` k=1 +!` do i=1,mp +!` if (ssnow%isflag(i) .eq. zero.and. soil%isoilm(i) .ne. 9) then +!` gammzz_snow(i,k) = real(Ccgsnow,r_2) * real(ssnow%snowd(i),r_2) +!` end if +!` end do +!` +!` sicefreeze(:,:) = -1._r_2 +!` sicemelt(:,:) = -1._r_2 +!` +!` do k=1,ms +!` do i=1,mp +!` +!` ssnow%wmice(i,k) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*real(Cdensity_ice,r_2) +!` ssnow%wmliq(i,k) = ssnow%wbliq(i,k)*soil%zse_vec(i,k)*real(Cdensity_liq,r_2) +!` ssnow%wmtot(i,k) = ssnow%wmice(i,k) + ssnow%wmliq(i,k) +!` +!` if ((ssnow%tgg(i,k) .lt. CTFRZ) .and. & +!` (ssnow%tgg(i,k) .lt. ssnow%otgg(i,k))) then +!` +!` ssnow%otgg(i,k) = min(ssnow%otgg(i,k),CTFRZ) +!` +!` if (ssnow%wb(i,k) .gt. (soil%watr(i,k)+1.0e-6) ) then +!` iceF(i,k) = max(zero,min(0.95_r_2,max(zero,ssnow%wbliq(i,k)-soil%watr(i,k))/& +!` (ssnow%wb(i,k)-soil%watr(i,k)) ) ) +!` else +!` iceF(i,k) = 0._r_2 +!` end if +!` tgg_tmp(i,k) = iceF(i,k)*ssnow%otgg(i,k) + & +!` (1._r_2 - iceF(i,k))*ssnow%tgg(i,k) +!` +!` Aconst =max(0.1,min(0.99_r_2,(ssnow%wb(i,k)-soil%watr(i,k))/(soil%ssat_vec(i,k)-soil%watr(i,k)) ) ) +!` Dconst = exp(1. - Aconst ) +!` +!` if (tgg_tmp(i,k) .lt. CTFRZ) then +!` max_ice_frac(i,k) = max(zero,ssnow%wb(i,k)*& +!` (1._r_2 - exp(2._r_2*(tgg_tmp(i,k)-CTFRZ)*Aconst*Aconst)) / Dconst) +!` +!` else +!` max_ice_frac(i,k) = zero +!` end if +!` +!` max_ice_frac(i,k) = min(ssnow%wb(i,k)*0.9, max_ice_frac(i,k) ) +!` +!` delta_ice_vol(i,k) = max(zero, max_ice_frac(i,k) - ssnow%wbice(i,k) ) +!` +!` !check amount of water we have +!` delta_ice_vol(i,k) = min(ssnow%wbliq(i,k)*real(Cdensity_liq/Cdensity_ice,r_2), delta_ice_vol(i,k) ) +!` +!` sicefreeze(i,k) = min(delta_ice_vol(i,k)*soil%zse_vec(i,k)*Cdensity_ice, & +!` max(zero,(ssnow%otgg(i,k)-ssnow%tgg(i,k))*ssnow%gammzz(i,k)/CHLF) ) +!` +!` elseif ((ssnow%tgg(i,k) .gt. CTFRZ) .and. & +!` (ssnow%tgg(i,k) .gt. ssnow%otgg(i,k)) .and. & +!` ssnow%wbice(i,k) .gt. zero ) then +!` +!` ssnow%otgg(i,k) = min(ssnow%otgg(i,k),CTFRZ) +!` +!` tgg_tmp(i,k) = CTFRZ +!` +!` delta_ice_vol(i,k) = ssnow%wbice(i,k) +!` +!` sicemelt(i,k) = min(delta_ice_vol(i,k)*soil%zse_vec(i,k)*Cdensity_ice, & +!` max(zero,(ssnow%tgg(i,k)-ssnow%otgg(i,k))*ssnow%gammzz(i,k)/CHLF)) +!` +!` endif +!` end do +!` end do +!` +!` DO k = 1, ms +!` DO i=1,mp +!` +!` +!` if (sicefreeze(i,k) .gt. zero .and. ssnow%tgg(i,k).lt.CTFRZ .and. & +!` ssnow%tgg(i,k).lt.ssnow%otgg(i,k) ) then +!` +!` +!` ssnow%wbice(i,k) = ssnow%wbice(i,k) +& +!` sicefreeze(i,k)/soil%zse_vec(i,k)/real(Cdensity_ice,r_2) +!` +!` delta_wbliq = max(zero,min(ssnow%wbliq(i,k),delta_ice_vol(i,k)*real(Cdensity_ice/Cdensity_liq,r_2))) +!` +!` ssnow%gammzz(i,k) = max(soil%heat_cap_lower_limit(i,k), & +!` (1.0- soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k) & +!` + (ssnow%wbliq(i,k) - delta_wbliq) * REAL(Ccswat*Cdensity_liq,r_2) & +!` + ssnow%wbice(i,k) * REAL(Ccsice*Cdensity_ice,r_2)& +!` )*soil%zse_vec(i,k) + gammzz_snow(i,k) +!` +!` ssnow%tgg(i,k) = ssnow%tgg(i,k) + real(sicefreeze(i,k) )& +!` * CHLF / real(ssnow%gammzz(i,k) ) +!` +!` ssnow%wmice(i,k) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*real(Cdensity_ice,r_2) +!` ssnow%wmliq(i,k) = ssnow%wmtot(i,k) - ssnow%wmice(i,k) +!` +!` ssnow%wbliq(i,k) = ssnow%wmliq(i,k) / (soil%zse_vec(i,k)*m2mm) +!` ssnow%wb(i,k) = ssnow%wbliq(i,k) + den_rat * ssnow%wbice(i,k) +!` +!` elseif (sicemelt(i,k) .gt. zero .and. ssnow%tgg(i,k).gt.CTFRZ & +!` .and. ssnow%wbice(i,k) .gt. 0._r_2) then +!` +!` ssnow%wbice(i,k) = ssnow%wbice(i,k) - delta_ice_vol(i,k) +!` +!` delta_wbliq = delta_ice_vol(i,k)*real(Cdensity_ice/Cdensity_liq,r_2) +!` +!` ssnow%gammzz(i,k) =max(soil%heat_cap_lower_limit(i,k), & +!` (1.0 - soil%ssat_vec(i,k))*soil%css_vec(i,k)*soil%rhosoil_vec(i,k)& +!` + (ssnow%wbliq(i,k)+delta_wbliq) * REAL(Ccswat*Cdensity_liq,r_2) & +!` + ssnow%wbice(i,k) * REAL(Ccsice*Cdensity_ice,r_2)& +!` )*soil%zse_vec(i,k) + gammzz_snow(i,k) +!` +!` ssnow%tgg(i,k) = ssnow%tgg(i,k) - real(sicemelt(i,k) )& +!` * CHLF / real(ssnow%gammzz(i,k) ) +!` +!` ssnow%wmice(i,k) = ssnow%wbice(i,k)*soil%zse_vec(i,k)*real(Cdensity_ice,r_2) +!` ssnow%wmliq(i,k) = ssnow%wmtot(i,k) - ssnow%wmice(i,k) +!` +!` ssnow%wbliq(i,k) = ssnow%wmliq(i,k) / (soil%zse_vec(i,k)*m2mm) +!` ssnow%wb(i,k) = ssnow%wbliq(i,k) + den_rat * ssnow%wbice(i,k) +!` +!` END IF +!` +!` +!` END DO +!` END DO +!` +!`END SUBROUTINE GWsoilfreeze + END MODULE cable_gw_hydro_module diff --git a/src/science/soilsnow/cbl_thermal.F90 b/src/science/soilsnow/cbl_thermal.F90 index 163b9266c2..cb682df1f5 100644 --- a/src/science/soilsnow/cbl_thermal.F90 +++ b/src/science/soilsnow/cbl_thermal.F90 @@ -6,7 +6,11 @@ MODULE snow_processes_soil_thermal_mod CONTAINS -SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) +!SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowmlt) ! replaced by rk4417 - phase2 +! snowmlt is not used in cable_gw_hydro_module and as far as I can tell serves no purpose in the code - rk4417 + +SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) +!* calculate snow processes and thermal soil USE snowl_adjust_mod, ONLY: snowl_adjust USE snowCheck_mod, ONLY: snowCheck @@ -24,10 +28,12 @@ SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal,snowml TYPE(veg_parameter_type), INTENT(INOUT) :: veg TYPE(met_type), INTENT(INOUT) :: met ! all met forcing TYPE (balances_type), INTENT(INOUT) :: bal - REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt - + REAL, DIMENSION(mp) :: snowmlt !track snow melt +! REAL, DIMENSION(:), INTENT(INOUT) :: snowmlt ! replaced by rk4417 - phase2 INTEGER :: k,i + snowmlt = 0.0 ! inserted by rk4417 - phase2 + CALL snowcheck (dels, ssnow, soil, met ) CALL snowdensity (dels, ssnow, soil) diff --git a/src/util/cable_common.F90 b/src/util/cable_common.F90 index ce9550d460..f3d2a99c97 100644 --- a/src/util/cable_common.F90 +++ b/src/util/cable_common.F90 @@ -100,36 +100,91 @@ MODULE cable_common_module ssat_vec_organic = 0.9, & watr_organic = 0.1, & sfc_vec_hk = 1.157407e-06, & - swilt_vec_hk = 2.31481481e-8 - + swilt_vec_hk = 2.31481481e-8, & + +! should the ones below be included also or is it a renaming? - rk4417 - phase2 +! they are included to fix compilation error - rk4417 - phase2 + hyds_organic = 1.0e-4, & + sucs_organic = 10.3, & + bch_organic = 2.91, & + ssat_organic = 0.9, & + css_organic = 4000.0, & + cnsd_organic = 0.1 END TYPE organic_soil_params TYPE gw_parameters_type +! REAL :: & +! MaxHorzDrainRate=2e-4, & !anisintropy * q_max [qsub] +! EfoldHorzDrainRate=2.0, & !e fold rate of q_horz +! MaxSatFraction=2500.0, & !parameter controll max sat fraction +! hkrz=0.5, & !hyds_vec variation with z +! zdepth=1.5, & !level where hyds_vec(z) = hyds_vec(no z) +! frozen_frac=0.05, & !ice fraction to determine first non-frozen layer for qsub +! SoilEvapAlpha = 1.0, & !modify field capacity dependence of soil evap limit +! IceAlpha=3.0, & +! IceBeta=1.0 +! +! REAL :: ice_impedence=5.0 +! +! TYPE(organic_soil_params) :: org +! +! INTEGER :: level_for_satfrac = 6 +! LOGICAL :: ssgw_ice_switch = .FALSE. +! +! LOGICAL :: subsurface_sat_drainage = .TRUE. + +! replaced above block by below - rk4417 - phase2 + REAL :: & MaxHorzDrainRate=2e-4, & !anisintropy * q_max [qsub] EfoldHorzDrainRate=2.0, & !e fold rate of q_horz + EfoldHorzDrainScale=1.0, & !e fold rate of q_horz MaxSatFraction=2500.0, & !parameter controll max sat fraction hkrz=0.5, & !hyds_vec variation with z zdepth=1.5, & !level where hyds_vec(z) = hyds_vec(no z) frozen_frac=0.05, & !ice fraction to determine first non-frozen layer for qsub SoilEvapAlpha = 1.0, & !modify field capacity dependence of soil evap limit IceAlpha=3.0, & - IceBeta=1.0 + IceBeta=1.0, & + sfc_vec_hk = 1.157407e-06, & + swilt_vec_hk = 2.31481481e-8 REAL :: ice_impedence=5.0 - + REAL :: ssat_wet_factor=0.85 + !hysteresis reduces ssat due to air entrapment + TYPE(organic_soil_params) :: org - + INTEGER :: aquifer_recharge_function=-1 !0=>no flux,1=>assume gw at hydrostat eq !inserted line as per MMY -- rk4417 INTEGER :: level_for_satfrac = 6 LOGICAL :: ssgw_ice_switch = .FALSE. + !LOGICAL :: derive_soil_param = .FALSE. ! MMY TRUE: derive soil parameters by cosby or HC-SWC equations hard-coded in CABLE-GW + ! MMY however, sand/silt/clay/org/rhosoil_vec are read from gridinfo + ! MMY FALSE: read soil parameters from land gridinfo file LOGICAL :: subsurface_sat_drainage = .TRUE. - + LOGICAL :: cosby_univariate=.false. + LOGICAL :: cosby_multivariate=.false. + LOGICAL :: HC_SWC=.false. !use Hutson Cass modified brooks corey + !seperates wet/dry to remove need for watr and + !gives better numerical behavoir for soln + LOGICAL :: BC_hysteresis=.false. + END TYPE gw_parameters_type TYPE(gw_parameters_type), SAVE :: gw_params +! psi_c and psi_o below inserted by rk4417 - phase2 + REAL, DIMENSION(17),SAVE :: psi_c = (/-2550000.0,-2550000.0,-2550000.0, & + -2240000.0,-4280000.0,-2750000.0,-2750000.0,& + -2750000.0,-2750000.0,-2750000.0,-2750000.0,-2750000.0,& + -2750000.0,-2750000.0,-2750000.0,-2750000.0,-2750000.0/) + + REAL, DIMENSION(17),SAVE :: psi_o = (/-66000.0,-66000.0,-66000.0,& + -35000.0,-83000.0,-74000.0,-74000.0,& + -74000.0,-74000.0,-74000.0,-74000.0,-74000.0,& + -74000.0,-74000.0,-74000.0,-74000.0,-74000.0/) + REAL, SAVE :: &!should be able to change parameters! max_glacier_snowd=1100.0,& snow_ccnsw = 2.0, & From 315b899be1ad36a03e408cd2efed4226ce1f688c Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 14:37:17 +1000 Subject: [PATCH 2/8] Update src/science/gw_hydro/cable_gw_hydro.F90 Co-authored-by: Claire Carouge --- src/science/gw_hydro/cable_gw_hydro.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index ead9055936..0b63c9c157 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -1,6 +1,6 @@ !============================================================================== ! This source code is part of the -! Austrian Community Atmosphere Biosphere Land Exchange (CABLE) model. +! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model. ! This work is licensed under the CABLE Academic User Licence Agreement ! (the "Licence"). ! You may not use this file except in compliance with the Licence. From f430650ece9e5f767104d63e0a386b4610f1a346 Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 14:41:29 +1000 Subject: [PATCH 3/8] Update src/science/gw_hydro/cable_gw_hydro.F90 Co-authored-by: Claire Carouge --- src/science/gw_hydro/cable_gw_hydro.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index 0b63c9c157..777093b1c3 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -109,18 +109,6 @@ SUBROUTINE GWsoilfreeze(dels, soil, ssnow) ! if acceptable this routine does the same thing but with r_2 soil params ! if max_ice_frac always set to frozen_limit and tgg_tmp is always CTFRZ ! - !## Method - ! - ! The subroutine takes a float as input, prints its value then adds 2/3, - ! prints it and returns it. - ! - ! The equation used is: - ! \[ myarg = myarg + !frac{2.}{3.} \] - ! - !## References - ! - ! The guidelines for documentation can be found in: - ! [CABLE developer guide](https://cable-lsm.github.io/CABLE/developer_guide/doc_guide/science_doc/) REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow From 1c445811e2cb74232cef6f36e0895c2dda83eb13 Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 14:42:43 +1000 Subject: [PATCH 4/8] Update src/science/gw_hydro/cable_gw_hydro.F90 Co-authored-by: Claire Carouge --- src/science/gw_hydro/cable_gw_hydro.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index 777093b1c3..29a4b774eb 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -905,18 +905,6 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) ! This is the main subroutine for Mark Decker's new hydrology scheme. ! This subroutine should be called in subroutine cbm but not yet. ! - !## Method - ! - ! The subroutine takes a float as input, prints its value then adds 2/3, - ! prints it and returns it. - ! - ! The equation used is: - ! \[ myarg = myarg + !frac{2.}{3.} \] - ! - !## References - ! - ! The guidelines for documentation can be found in: - ! [CABLE developer guide](https://cable-lsm.github.io/CABLE/developer_guide/doc_guide/science_doc/) ! ! Inputs: From 6fac982362a111210d807955eb45e25fdd66d27e Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 15:35:56 +1000 Subject: [PATCH 5/8] Update src/science/gw_hydro/cable_gw_hydro.F90 --- src/science/gw_hydro/cable_gw_hydro.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index 29a4b774eb..35d377f692 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -1385,7 +1385,7 @@ END SUBROUTINE calc_equilibrium_water_content SUBROUTINE calc_srf_wet_fraction(ssnow,soil,met,veg) !* Calculate the wet fraction of the surafce ! following [Decker, 2015](http://doi.wiley.com/10.1002/2015MS000507) - !* This subrountine is called in subrountine surf_wetness_fact in cable_canopy.F90 + IMPLICIT NONE TYPE(soil_snow_type), INTENT(INOUT) :: ssnow ! soil+snow variables From 45251dd5e4b13c74f6e7e73e6302006ade6bcf11 Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 15:43:37 +1000 Subject: [PATCH 6/8] Update src/offline/cable_parameters.F90 --- src/offline/cable_parameters.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index e69ead88c9..be0f14f184 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -1798,7 +1798,7 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: ssat_bounded,rho_soil_bulk ! added by rk4417 - phase2 ! line below inserted by rk4417 - phase2 - where(veg%iveg .eq. 17) soil%isoilm = 9 ! MMY@13April it says where iveg = ice, isoilm = permanent ice + ! soil_depth(1) = REAL(soil%zse(1),r_2) ! DO klev=2,ms From 7fdcd1166c4678ee5c4f11d908e10b96eeaf64bc Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 15:45:29 +1000 Subject: [PATCH 7/8] Update src/offline/cable_parameters.F90 --- src/offline/cable_parameters.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index be0f14f184..6f6792d6ff 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -1797,7 +1797,6 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) REAL(r_2), DIMENSION(:,:), ALLOCATABLE :: ssat_bounded,rho_soil_bulk ! added by rk4417 - phase2 -! line below inserted by rk4417 - phase2 ! soil_depth(1) = REAL(soil%zse(1),r_2) From fb4cc3b5d0e5f81591ae39cc185e34699b1a47e3 Mon Sep 17 00:00:00 2001 From: Ramzi Kutteh <98803952+rkutteh@users.noreply.github.com> Date: Thu, 25 Jul 2024 18:22:02 +1000 Subject: [PATCH 8/8] Update src/offline/cable_parameters.F90 Co-authored-by: Claire Carouge --- src/offline/cable_parameters.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 6f6792d6ff..2780d8baf5 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -2288,7 +2288,7 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) END DO bal%osnowd0 = ssnow%osnowd - !! vh_js !! comment out hide% condition + ! vh_js ! comment out hide% condition ! IF (hide%Ticket49Bug6) THEN IF(cable_user%SOIL_STRUC=='sli') THEN