diff --git a/src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual b/src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual deleted file mode 100644 index b1548b0d8..000000000 --- a/src/coupled/AM3/control/cable/util/init/cbl_um_init_veg.F90.manual +++ /dev/null @@ -1,77 +0,0 @@ -MODULE cbl_um_init_veg_mod - -IMPLICIT NONE - -PUBLIC initialize_veg - -CONTAINS - -SUBROUTINE initialize_veg( SurfaceType, SoilType, mp,& - ms, nrb, npft, nsurft, land_pts, l_tile_pts, & - ICE_SurfaceType, ICE_SoilType, LAI_thresh, & - soil_zse, surft_pts, surft_index, tile_frac, & - veg, soil, pars ) - -USE init_cable_parms_mod, ONLY: init_cable_parms_expl -USE params_io_mod_cbl, ONLY: params_io_data_type -USE cable_def_types_mod, ONLY: veg_parameter_type -USE cable_def_types_mod, ONLY: soil_parameter_type, r_2 - -IMPLICIT NONE - -INTEGER, INTENT(IN) :: mp -INTEGER, INTENT(IN) :: ms -INTEGER, INTENT(IN) :: nrb -INTEGER, INTENT(IN) :: nsurft -INTEGER, INTENT(IN) :: npft -INTEGER, INTENT(IN) :: land_pts -INTEGER, INTENT(IN) :: ICE_SurfaceType ! index ICE surface type -INTEGER, INTENT(IN) :: ICE_SoilType ! index soil type -REAL, INTENT(IN) :: lai_thresh -REAL, INTENT(IN) :: soil_zse(ms) ! soil depth per layer -INTEGER, INTENT(IN) :: surft_pts(nsurft) ! # land points per tile -INTEGER, INTENT(IN) :: surft_index(land_pts, nsurft) ! land_pt index of point -LOGICAL, INTENT(IN) :: L_tile_pts(land_pts, nsurft) ! TRUE if active tile -REAL, INTENT(IN) :: tile_frac(land_pts,nsurft) - -INTEGER, INTENT(OUT) :: SurfaceType(mp) ! surface tile PFT/nveg -INTEGER, INTENT(OUT) :: SoilType(mp) ! soil type per tile - -TYPE(veg_parameter_type), INTENT(OUT) :: veg ! vegetation parameters -TYPE(soil_parameter_type), INTENT(OUT) :: soil ! soil parameters -TYPE(params_io_data_type), INTENT(IN) :: pars - -! local vars -INTEGER :: i -REAL :: totdepth - -CALL init_cable_parms_expl( SurfaceType, SoilType, mp, ms, nrb, land_pts, & - nsurft, l_tile_pts, ICE_SurfaceType, & - ICE_SoilType, soil_zse, veg, soil, pars, tile_frac ) - -veg%meth = 1 ! review !should be namelist config if even used -veg%ejmax = 2.*veg%vcmax -veg%gamma = 3.e-2 ! for Haverd2013 switch ! offline set init _parameters -veg%F10 = 0.85 ! offline set init _parameters - -! calculate veg%froot from using rootbeta and soil depth -! (Jackson et al. 1996, Oceologica, 108:389-411) -totdepth = 0.0 -DO i = 1, ms-1 - totdepth = totdepth + soil_zse(i) * 100.0 ! unit in centimetres - veg%froot(:, i) = MIN( 1.0_r_2, 1.0-veg%rootbeta(:)**totdepth ) -END DO -veg%froot(:, ms) = 1.0 - veg%froot(:, ms-1) -DO i = ms-1, 2, -1 - veg%froot(:, i) = veg%froot(:, i)-veg%froot(:,i-1) -END DO - - -RETURN -END SUBROUTINE initialize_veg - -END MODULE cbl_um_init_veg_mod - - - - diff --git a/src/offline/cable_checks.F90 b/src/offline/cable_checks.F90 index 5d1b73a5e..5ad32517c 100644 --- a/src/offline/cable_checks.F90 +++ b/src/offline/cable_checks.F90 @@ -474,48 +474,46 @@ END SUBROUTINE constant_check_range SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met, & air,bal) + !* Calculates the energy and water balances at each time step as well as + ! the cumulative balance over the simulation. ! Input arguments - REAL,INTENT(IN) :: dels ! time step size - INTEGER, INTENT(IN) :: ktau ! timestep number - TYPE (soil_snow_type),INTENT(IN) :: ssnow ! soil data - TYPE (soil_parameter_type),INTENT(IN) :: soil ! soil data - TYPE (canopy_type),INTENT(IN) :: canopy ! canopy variable data - TYPE(met_type),INTENT(IN) :: met ! met data + REAL,INTENT(IN) :: dels + !! Time step length \(s\) + INTEGER, INTENT(IN) :: ktau + !! Number of time steps since the start of the simulation \(-\) + TYPE (soil_snow_type),INTENT(IN) :: ssnow + !! Soil and snow variables data + TYPE (soil_parameter_type),INTENT(IN) :: soil + !! Soil parameters + TYPE (canopy_type),INTENT(IN) :: canopy + !! Canopy variables data + TYPE(met_type),INTENT(IN) :: met + !! Met forcing data TYPE (air_type),INTENT(IN) :: air + !! Air variables data + TYPE (balances_type),INTENT(INOUT) :: bal + !! Water balance variables data ! Local variables - REAL(r_2), DIMENSION(:,:,:),POINTER, SAVE :: bwb ! volumetric soil moisture + REAL(r_2), DIMENSION(:),POINTER, SAVE :: owb ! volumetric soil moisture + ! at previous time step REAL(r_2), DIMENSION(mp) :: delwb ! change in soilmoisture ! b/w tsteps - REAL, DIMENSION(mp) :: canopy_wbal !canopy water balance - TYPE (balances_type),INTENT(INOUT) :: bal - INTEGER :: j, k ! do loop counter + REAL, DIMENSION(mp) :: canopy_wbal !canopy water balance + INTEGER :: j, k ! do loop counter IF(ktau==1) THEN - ALLOCATE( bwb(mp,ms,2) ) + ALLOCATE( owb(mp) ) ! initial vlaue of soil moisture - bwb(:,:,1)=ssnow%wb - bwb(:,:,2)=ssnow%wb - delwb(:) = 0. - ELSE - ! Calculate change in soil moisture b/w timesteps: - IF(MOD(REAL(ktau),2.0)==1.0) THEN ! if odd timestep - bwb(:,:,1)=ssnow%wb - DO k=1,mp ! current smoist - prev tstep smoist - delwb(k) = SUM((bwb(k,:,1) & - - (bwb(k,:,2)))*soil%zse)*1000.0 - END DO - ELSE IF(MOD(REAL(ktau),2.0)==0.0) THEN ! if even timestep - bwb(:,:,2)=ssnow%wb - DO k=1,mp ! current smoist - prev tstep smoist - delwb(k) = SUM((bwb(k,:,2) & - - (bwb(k,:,1)))*soil%zse)*1000.0 - END DO - END IF + owb=ssnow%wbtot END IF + ! Calculate change in soil moisture b/w timesteps: + delwb = ssnow%wbtot - owb + ! Update old soil moisture to this timestep value. + owb = ssnow%wbtot ! net water into soil (precip-(change in canopy water storage) ! - (change in snow depth) - (surface runoff) - (deep drainage) diff --git a/src/params/cable_params_mod.F90 b/src/params/cable_params_mod.F90 deleted file mode 100644 index 8bd5f73ec..000000000 --- a/src/params/cable_params_mod.F90 +++ /dev/null @@ -1,255 +0,0 @@ -MODULE cable_params_mod - -!H! Elevate these to namelist definable -USE cable_other_constants_mod, ONLY: nsl -USE cable_other_constants_mod, ONLY: nrb ! # radiation "bANDS" - !dir/dif components in bands VIS/NIR -USE cable_other_constants_mod, ONLY: nscs ! number of soil carbon stores -USE cable_other_constants_mod, ONLY: nvcs ! number of vegetation carbon stores -USE grid_constants_mod_cbl, ONLY : mstype => nsoil_max ! # of soil types [9] -USE grid_constants_mod_cbl, ONLY : ntype_max ! # of PFTs [17] - -IMPLICIT NONE - -PUBLIC :: veg_parameter_type -PUBLIC :: vegin_type -PUBLIC :: soil_parameter_type -PUBLIC :: soilin_type - -PUBLIC :: vegin -PUBLIC :: soilin - -!----------------------------------------------------------------------------- -! Description: -! Defines variable types and variables for CABLE standalone runs. -! Based on cable_def_types_mod.F90 from the CABLE trunk. -! -! Code Owner: Please refer to ModuleLeaders.txt -! This file belongs in CABLE SCIENCE -!----------------------------------------------------------------------------- - -! Vegetation parameters used: -TYPE veg_parameter_type -!jhan: surely %meth doesn't need to be 'mp' here - INTEGER, DIMENSION(:), POINTER :: & - meth, & ! method for calculation of canopy fluxes and temp. - iveg , & ! vegetation type - iLU ! land use type - REAL, DIMENSION(:), POINTER :: & - canst1, & ! max intercepted water by canopy (mm/LAI) - dleaf, & ! chararacteristc legnth of leaf (m) - ejmax, & ! max pot. electron transp rate top leaf(mol/m2/s) - frac4, & ! fraction of c4 plants - hc, & ! roughness height of canopy (veg - snow) - vlai, & ! leaf area index - xalbnir, & - rp20, & ! plant respiration coefficient at 20 C - rpcoef, & ! temperature coef nonleaf plant respiration (1/C) - rs20, & ! soil respiration at 20 C [mol m-2 s-1] - shelrb, & ! sheltering factor (dimensionless) - vegcf, & ! kdcorbin, 08/10 - tminvj, & ! min temperature of the start of photosynthesis - toptvj, & ! opt temperature of the start of photosynthesis - tmaxvj, & ! max temperature of the start of photosynthesis - vbeta, & ! - vcmax, & ! max RuBP carboxylation rate top leaf (mol/m2/s) - xfang, & ! leaf angle PARAMETER - extkn, & ! extinction coef for vertical - vlaimax, & ! extinction coef for vertical - wai, & ! wood area index (stem+branches+twigs) - a1gs, & ! a1 parameter in stomatal conductance model - d0gs, & ! d0 in stomatal conductance model - alpha, & ! initial slope of J-Q response curve - convex, & ! convexity of J-Q response curve - cfrd, & ! ratio of day respiration to vcmax - gswmin, & ! minimal stomatal conductance - conkc0, & ! Michaelis-menton constant for carboxylase - conko0, & ! Michaelis-menton constant for oxygenase - ekc, & ! activation energy for caroxylagse - eko, & ! acvtivation enegery for oxygenase - g0, & ! Belinda's stomatal model intercept, Ticket #56. - g1 ! Belinda's stomatal model slope, Ticket #56. - - LOGICAL, DIMENSION(:), POINTER :: & - deciduous ! flag used for phenology fix - - REAL, DIMENSION(:,:), POINTER :: & - refl, & - taul, & - froot ! fraction of root in each soil layer - - ! Additional veg parameters: - REAL, DIMENSION(:), POINTER :: rootbeta ! parameter for estimating vertical root mass distribution (froot) - REAL, DIMENSION(:), POINTER :: gamma ! parameter in root efficiency function (Lai and Katul 2000) - REAL, DIMENSION(:), POINTER :: ZR ! maximum rooting depth (cm) - REAL, DIMENSION(:), POINTER :: F10 ! fraction of roots in top 10 cm - - REAL, DIMENSION(:), POINTER :: clitt ! - - ! Additional POP veg param - INTEGER, DIMENSION(:,:), POINTER :: disturbance_interval - REAL, DIMENSION(:,:), POINTER :: disturbance_intensity - - REAL, POINTER :: & - soil(:,:), & - csoil(:,:), & - cplant(:,:), & - ratecs(:,:), & - ratecp(:,:) - - END TYPE veg_parameter_type - - -! Vegetation parameters I/O: -TYPE vegin_type - REAL :: & - canst1(ntype_max), & - dleaf(ntype_max), & - length(ntype_max), & - width(ntype_max), & - vcmax(ntype_max), & - ejmax(ntype_max), & - hc(ntype_max), & - xfang(ntype_max), & - rp20(ntype_max), & - rpcoef(ntype_max), & - rs20(ntype_max), & - wai(ntype_max), & - rootbeta(ntype_max), & - shelrb(ntype_max), & - vegcf(ntype_max), & - frac4(ntype_max), & - xalbnir(ntype_max), & - extkn(ntype_max), & - tminvj(ntype_max), & - tmaxvj(ntype_max), & - vbeta(ntype_max), & - a1gs(ntype_max), & - d0gs(ntype_max), & - alpha(ntype_max), & - convex(ntype_max), & - cfrd(ntype_max), & - gswmin(ntype_max), & - conkc0(ntype_max), & - conko0(ntype_max), & - ekc(ntype_max), & - eko(ntype_max), & - g0(ntype_max), & - g1(ntype_max), & - zr(ntype_max), & - clitt(ntype_max), & - froot(nsl,ntype_max), & - csoil(nscs,ntype_max), & - ratecs(nscs,ntype_max), & - cplant(nvcs,ntype_max), & - ratecp(nvcs,ntype_max), & - refl(nrb,ntype_max), & - taul(nrb,ntype_max) -END TYPE vegin_type - -! Soil parameters used: - TYPE soil_parameter_type - - INTEGER, DIMENSION(:), POINTER :: & - isoilm ! integer soil type - - REAL, DIMENSION(:), POINTER :: & - bch, & ! parameter b in Campbell equation - c3, & ! c3 drainage coeff (fraction) - clay, & ! fraction of soil which is clay - css, & ! soil specific heat capacity [kJ/kg/K] - hsbh, & ! difsat * etasat (=hyds*abs(sucs)*bch) - hyds, & ! hydraulic conductivity @ saturation [m/s], Ksat - i2bp3, & ! par. one in K vis suction (=nint(bch)+2) - ibp2, & ! par. two in K vis suction (fn of pbch) - rhosoil, & ! soil density [kg/m3] - sand, & ! fraction of soil which is sand - sfc, & ! vol H2O @ field capacity - silt, & ! fraction of soil which is silt - ssat, & ! vol H2O @ saturation - sucs, & ! suction at saturation (m) - swilt, & ! vol H2O @ wilting - zse, & ! thickness of each soil layer (1=top) in m - zshh, & ! distance between consecutive layer midpoints (m) - ! vars intro for Ticket #27 - soilcol, & ! keep color for all patches/tiles - albsoilf ! soil reflectance - - REAL, DIMENSION(:,:), POINTER :: & - heat_cap_lower_limit - - REAL, DIMENSION(:,:), POINTER :: & - zse_vec,css_vec,cnsd_vec - - REAL, DIMENSION(:), POINTER :: & - cnsd, & ! thermal conductivity of dry soil [W/m/K] - pwb_min ! working variable (swilt/ssat)**ibp2 - - REAL, DIMENSION(:,:), POINTER :: & - albsoil ! soil reflectance (2nd dim. BP 21Oct2009) - !mrd561 - !MD parameters for GW module that vary with soil layer - REAL, DIMENSION(:,:), POINTER :: & - sucs_vec, & !psi at saturation in [mm] - hyds_vec, & !saturated hydraulic conductivity [mm/s] - bch_vec, & !C and H B [none] - clay_vec, & !fraction of soil that is clay [frac] - sand_vec, & !fraction of soil that is sand [frac] - silt_vec, & !fraction of soil that is silt [frac] - org_vec, & !fration of soil made of organic soils [frac] - rhosoil_vec,& !soil density [kg/m3] - ssat_vec, & !volumetric water content at saturation [mm3/mm3] - watr, & !residual water content of the soil [mm3/mm3] - sfc_vec, & !field capcacity (hk = 1 mm/day) - swilt_vec ! wilting point (hk = 0.02 mm/day) - - REAL, DIMENSION(:), POINTER :: & - drain_dens,&! drainage density ( mean dist to rivers/streams ) - elev, & !elevation above sea level - elev_std, & !elevation above sea level - slope, & !mean slope of grid cell - slope_std !stddev of grid cell slope - - !MD parameters for GW module for the aquifer - REAL, DIMENSION(:), POINTER :: & - GWsucs_vec, & !head in the aquifer [mm] - GWhyds_vec, & !saturated hydraulic conductivity of the aquifer [mm/s] - GWbch_vec, & !clapp and horn b of the aquifer [none] - 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] - GWdz, & !thickness of the aquifer [m] - GWrhosoil_vec !density of the aquifer substrate [kg/m3] - - ! Additional SLI parameters - INTEGER, DIMENSION(:), POINTER :: nhorizons ! number of soil horizons - INTEGER, DIMENSION(:,:), POINTER :: ishorizon ! horizon number 1:nhorizons - REAL, DIMENSION(:), POINTER :: clitt ! litter (tC/ha) - REAL, DIMENSION(:), POINTER :: zeta ! macropore parameter - REAL, DIMENSION(:), POINTER :: fsatmax ! variably saturated area parameter - - END TYPE soil_parameter_type - - - -! Soil parameters I/O: -TYPE soilin_type - REAL :: & - silt(mstype), & - clay(mstype), & - sand(mstype), & - swilt(mstype), & - sfc(mstype), & - ssat(mstype), & - bch(mstype), & - hyds(mstype), & - sucs(mstype), & - rhosoil(mstype), & - css(mstype) -END TYPE soilin_type - -!Instantiate types -TYPE(vegin_type) :: vegin !read from namelist -TYPE(soilin_type) :: soilin !read from namelist - -END MODULE cable_params_mod diff --git a/src/science/soilsnow/cbl_hyd_redistrib.F90 b/src/science/soilsnow/cbl_hyd_redistrib.F90 index 0c63f61e3..ad63b8ef3 100644 --- a/src/science/soilsnow/cbl_hyd_redistrib.F90 +++ b/src/science/soilsnow/cbl_hyd_redistrib.F90 @@ -13,6 +13,8 @@ MODULE hydraulic_redistribution_mod SUBROUTINE hydraulic_redistribution(dels, soil, ssnow, canopy, veg, met) USE cable_common_module, ONLY : wiltParam, satuParam +!data +USE cable_surface_types_mod, ONLY: evergreen_broadleaf, c4_grassland IMPLICIT NONE REAL, INTENT(IN) :: dels ! integration time step (s) @@ -97,10 +99,9 @@ SUBROUTINE hydraulic_redistribution(dels, soil, ssnow, canopy, veg, met) hr_perTime(:,k,j) = hr_perTime(:,k,j)/soil%zse(k) hr_perTime(:,j,k) = hr_perTime(:,j,k)/soil%zse(j) - ! Overwrite to give zero redistribution for all types except - ! evergreen broadleaf (2) and c4 grass (7) - ! NB: Hard-wired numbers should be removed in future version - WHERE( .NOT.(veg%iveg == 2 .OR. veg%iveg == 7 ) ) + ! Zero redistribution: all types except evergreen broadleaf, c4 grass + WHERE( .NOT. ( veg%iveg == evergreen_broadleaf .OR. & + veg%iveg == c4_grassland ) ) hr_perTime(:,k,j) = 0.0 hr_perTime(:,j,k) = 0.0 ENDWHERE @@ -175,10 +176,9 @@ SUBROUTINE hydraulic_redistribution(dels, soil, ssnow, canopy, veg, met) hr_perTime(:,k,j) = hr_perTime(:,k,j)/soil%zse(k) hr_perTime(:,j,k) = hr_perTime(:,j,k)/soil%zse(j) - ! Overwrite to give zero redistribution for all types except - ! evergreen broadleaf (2) and c4 grass (7) - ! NB: Hard-wired numbers should be removed in future version - WHERE( .NOT.( veg%iveg == 2 .OR. veg%iveg == 7 ) ) + ! Zero redistribution: all types except evergreen broadleaf, c4 grass + WHERE( .NOT. ( veg%iveg == evergreen_broadleaf .OR. & + veg%iveg == c4_grassland ) ) hr_perTime(:,k,j) = 0.0 hr_perTime(:,j,k) = 0.0 ENDWHERE diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 6274cedbc..50a438ad8 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -31,14 +31,14 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) ! which can be removed: xx = canopy%fevc * dels / CHL * veg%froot(:,k) + diff(:,k-1) ! kg/m2 diff(:,k) = MAX( 0.0_r_2, ssnow%wb(:,k) - soil%swilt) & ! m3/m3 - * soil%zse(k)*1000.0 + * soil%zse(k)*Cdensity_liq xxd = xx - diff(:,k) WHERE ( xxd .GT. 0.0 ) - ssnow%wb(:,k) = ssnow%wb(:,k) - diff(:,k) / (soil%zse(k)*1000.0) + ssnow%wb(:,k) = ssnow%wb(:,k) - diff(:,k) / (soil%zse(k)*Cdensity_liq) diff(:,k) = xxd ELSEWHERE - ssnow%wb(:,k) = ssnow%wb(:,k) - xx / (soil%zse(k)*1000.0) + ssnow%wb(:,k) = ssnow%wb(:,k) - xx / (soil%zse(k)*Cdensity_liq) diff(:,k) = 0.0 ENDWHERE @@ -52,10 +52,7 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) canopy%fevc = 0.0_r_2 END WHERE DO k = 1,ms - ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*1000.0) - - ! write(59,*) k, ssnow%wb(:,k), ssnow%evapfbl(:,k)/(soil%zse(k)*1000.0) - ! write(59,*) + ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq) ENDDO diff --git a/src/science/soilsnow/cbl_smoisturev.F90 b/src/science/soilsnow/cbl_smoisturev.F90 index ac38c482a..635e4f06b 100644 --- a/src/science/soilsnow/cbl_smoisturev.F90 +++ b/src/science/soilsnow/cbl_smoisturev.F90 @@ -60,9 +60,9 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg) wbl_kp, & ! wh, & ! z3_k, & ! - pwb_wbh ! - + pwb_wbh + REAL :: dfactor, sicemelt REAL, DIMENSION(mp,ms+1) :: z1mult ! change dimension of at,bt,ct from 3*ms to ms (BP Jun2010) @@ -83,10 +83,7 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg) delt, & ! dtt ! - LOGICAL :: is_open ! Is file open? - INTEGER :: & - u, & ! I/O unit - k + INTEGER :: i, k at = 0.0 bt = 1.0 @@ -225,7 +222,7 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg) END DO - ssnow%rnof2 = dels * REAL( fluxh(:,ms) ) * 1000.0 + ssnow%rnof2 = dels * REAL( fluxh(:,ms) ) * Cdensity_liq ! wbh_k represents wblf(k-.5) DO k = 2, ms @@ -425,9 +422,23 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg) CALL trimb(at, bt, ct, ssnow%wblf, ms) DO k = 1, ms - ssatcurr(:,k) = soil%ssat - ssnow%wbice(:,k) - ssnow%wb(:,k) = ssnow%wblf(:,k) * ssatcurr(:,k) + ssnow%wbice(:,k) - ssnow%wbice(:,k) = MIN( ssnow%wbice(:,k), frozen_limit * ssnow%wb(:,k) ) + ssatcurr(:,k) = soil%ssat - ssnow%wbice(:,k) + ssnow%wb(:,k) = ssnow%wblf(:,k) * ssatcurr(:,k) + ssnow%wbice(:,k) + END DO + ! If redistribution of liquid has caused a layer to exceed frozen limit + ! melt the excess and account for latent heat. + dfactor = 1.0 - Cdensity_ice / Cdensity_liq + DO k = 1, ms + DO i = 1, mp + IF (ssnow%wbice(i,k) > frozen_limit * ssnow%wb(i,k)) THEN + sicemelt = (ssnow%wbice(i,k) - frozen_limit * ssnow%wb(i,k)) / & + (1.0 - frozen_limit*dfactor) + ssnow%wbice(i,k) = ssnow%wbice(i,k) - sicemelt + ssnow%wb(i,k) = ssnow%wb(i,k) - dfactor*sicemelt + ssnow%tgg(i,k) = ssnow%tgg(i,k) - & + sicemelt*soil%zse(k)*Cdensity_ice*CHLF / REAL(ssnow%gammzz(i,k) ) + END IF + END DO END DO END SUBROUTINE smoisturev diff --git a/src/science/soilsnow/cbl_snowAccum.F90 b/src/science/soilsnow/cbl_snowAccum.F90 index 28b6a3dd6..27ed28d03 100644 --- a/src/science/soilsnow/cbl_snowAccum.F90 +++ b/src/science/soilsnow/cbl_snowAccum.F90 @@ -47,6 +47,8 @@ SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil ) ssnow%tgg(i,1) = ssnow%tgg(i,1) + canopy%precis(i) * CHLF & / ( REAL( ssnow%gammzz(i,1) ) + Ccswat *canopy%precis(i) ) + ssnow%dtmlt(i,1) = ssnow%dtmlt(i,1) + canopy%precis(i) * CHLF & + / ( REAL( ssnow%gammzz(i,1) ) + Ccswat *canopy%precis(i) ) ! change density due to water being added ssnow%ssdn(i,1) = MIN( max_ssdn, MAX( 120.0, ssnow%ssdn(i,1) & * ssnow%osnowd(i) / MAX( 0.01, ssnow%snowd(i) ) + Cdensity_liq & @@ -87,6 +89,8 @@ SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil ) ssnow%tggsn(i,1) = ssnow%tggsn(i,1) + canopy%precis(i) * CHLF & * osm(i) / (sgamm(i) * ssnow%osnowd(i) ) + ssnow%dtmlt(i,1) = ssnow%dtmlt(i,1) + canopy%precis(i) * CHLF & + * osm(i) / (sgamm(i) * ssnow%osnowd(i) ) ssnow%smass(i,1) = ssnow%smass(i,1) + canopy%precis(i) & * osm(i)/ssnow%osnowd(i) @@ -173,11 +177,8 @@ SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil ) ssnow%sdepth(i,1) = MAX( 0.02, ssnow%smass(i,1) / ssnow%ssdn(i,1) ) ENDIF - canopy%segg(i) = 0.0 - - !INH: cls package - !we still need to conserve moisture/energy when evapsn is limited - !this is a key point of moisture non-conservation + ! Use any remaining energy for evaporation of soil water + canopy%segg(i) = (CHL + CHLF) * (xxx(i) - ssnow%evapsn(i)) / CHL / dels ENDIF diff --git a/src/science/soilsnow/cbl_soilfreeze.F90 b/src/science/soilsnow/cbl_soilfreeze.F90 index ec068b2b0..a1b1a3cde 100644 --- a/src/science/soilsnow/cbl_soilfreeze.F90 +++ b/src/science/soilsnow/cbl_soilfreeze.F90 @@ -26,16 +26,22 @@ SUBROUTINE soilfreeze(dels, soil, ssnow,heat_cap_lower_limit) IF(ssnow%tgg(i,k) < CTFRZ & & .AND. frozen_limit*ssnow%wb(i,k) - ssnow%wbice(i,k) > .001) THEN - sicefreeze(i) = MIN( MAX( 0.0_r_2, ( frozen_limit * ssnow%wb(i,k) - & - ssnow%wbice(i,k) ) ) * soil%zse(k) * 1000.0, & - ( CTFRZ - ssnow%tgg(i,k) ) * ssnow%gammzz(i,k) / CHLF ) + ! Limits on water mass that can be converted to ice due frozen_limit and saturation + sicefreeze(i) = MIN( frozen_limit * ssnow%wb(i,k) - ssnow%wbice(i,k) , & + (soil%ssat(i) - ssnow%wb(i,k)) / MAX((1.0-Cdensity_ice/Cdensity_liq),1.0E-3) ) + ! Apply limits on watermass to freeze by zero and the energy limit + sicefreeze(i) = MIN( MAX( 0.0_r_2, sicefreeze(i) ) * soil%zse(k) * Cdensity_ice, & + ( CTFRZ - ssnow%tgg(i,k) ) * ssnow%gammzz(i,k) / CHLF ) + ssnow%wbice(i,k) = MIN( ssnow%wbice(i,k) + sicefreeze(i) / (soil%zse(k) & - * 1000.0), frozen_limit * ssnow%wb(i,k) ) + * Cdensity_ice), frozen_limit * ssnow%wb(i,k) ) + ssnow%wb(i,k) = ssnow%wb(i,k) + sicefreeze(i) / (soil%zse(k) * Cdensity_ice) & + - sicefreeze(i) / (soil%zse(k) * Cdensity_liq) xx(i) = soil%css(i) * soil%rhosoil(i) max_arg1 = heat_cap_lower_limit(i,k) max_arg2 = REAL((1.0 - soil%ssat(i)) * soil%css(i) * soil%rhosoil(i) ,r_2) & + (ssnow%wb(i,k) - ssnow%wbice(i,k)) * REAL(Ccswat * Cdensity_liq,r_2) & - + ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_liq * 0.9,r_2) + + ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_ice,r_2) ssnow%gammzz(i,k) = MAX( max_arg1, max_arg2 ) * REAL( soil%zse(k),r_2 ) IF (k == 1 .AND. ssnow%isflag(i) == 0) & @@ -46,16 +52,18 @@ SUBROUTINE soilfreeze(dels, soil, ssnow,heat_cap_lower_limit) ELSEIF( ssnow%tgg(i,k) > CTFRZ .AND. ssnow%wbice(i,k) > 0. ) THEN - sicemelt(i) = MIN( ssnow%wbice(i,k) * soil%zse(k) * 1000.0, & + sicemelt(i) = MIN( ssnow%wbice(i,k) * soil%zse(k) * Cdensity_ice, & ( ssnow%tgg(i,k) - CTFRZ ) * ssnow%gammzz(i,k) / CHLF ) ssnow%wbice(i,k) = MAX( 0.0_r_2, ssnow%wbice(i,k) - sicemelt(i) & - / (soil%zse(k) * 1000.0) ) + / (soil%zse(k) * Cdensity_ice) ) + ssnow%wb(i,k) = ssnow%wb(i,k) - sicemelt(i) / (soil%zse(k) * Cdensity_ice) & + + sicemelt(i) / (soil%zse(k) * Cdensity_liq) xx(i) = soil%css(i) * soil%rhosoil(i) max_arg1 = heat_cap_lower_limit(i,k) max_arg2 = REAL((1.0 - soil%ssat(i)) * soil%css(i) * soil%rhosoil(i) ,r_2) & + (ssnow%wb(i,k) - ssnow%wbice(i,k)) * REAL(Ccswat * Cdensity_liq,r_2) & - + ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_liq * 0.9,r_2) + + ssnow%wbice(i,k) * REAL(Ccsice * Cdensity_ice,r_2) ssnow%gammzz(i,k) = MAX( max_arg1, max_arg2 ) * REAL( soil%zse(k),r_2 ) diff --git a/src/science/soilsnow/cbl_soilsnow_data.F90 b/src/science/soilsnow/cbl_soilsnow_data.F90 index a3a34060c..88ac8adc5 100644 --- a/src/science/soilsnow/cbl_soilsnow_data.F90 +++ b/src/science/soilsnow/cbl_soilsnow_data.F90 @@ -8,6 +8,7 @@ MODULE cbl_ssnow_data_mod USE cable_phys_constants_mod, ONLY : CHL => HL USE cable_phys_constants_mod, ONLY : CHLF => HLF 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 : Ccsice => csice diff --git a/src/science/soilsnow/cbl_soilsnow_init_special.F90 b/src/science/soilsnow/cbl_soilsnow_init_special.F90 index 816ca797b..69e25503e 100644 --- a/src/science/soilsnow/cbl_soilsnow_init_special.F90 +++ b/src/science/soilsnow/cbl_soilsnow_init_special.F90 @@ -27,15 +27,12 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap REAL, DIMENSION(mp) :: xx, tgg_old, tggsn_old REAL(r_2), DIMENSION(mp) :: xxx,deltat,sinfil1,sinfil2,sinfil3 REAL :: zsetot -INTEGER, SAVE :: ktau =0 REAL :: heat_cap_lower_limit(mp,ms) -ktau = ktau +1 +! since CMIP5 only ever (potentially) TRUE offline if special initialization +! requested and then only on first timestep +IF( .NOT.cable_user%cable_runtime_coupled ) THEN -IF( .NOT.cable_user%cable_runtime_coupled ) THEN - - IF( ktau_gl <= 1 ) 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 @@ -64,27 +61,29 @@ SUBROUTINE spec_init_soil_snow(dels, soil, ssnow, canopy, met, bal, veg,heat_cap ssnow%wb(:,4) = 0.95 * soil%ssat ssnow%wb(:,5) = 0.95 * soil%ssat ssnow%wb(:,6) = 0.95 * soil%ssat - ssnow%wbice(:,1) = 0.90 * ssnow%wb(:,1) - ssnow%wbice(:,2) = 0.90 * ssnow%wb(:,2) - ssnow%wbice(:,3) = 0.90 * ssnow%wb(:,3) - ssnow%wbice(:,4) = 0.90 * ssnow%wb(:,4) - ssnow%wbice(:,5) = 0.90 * ssnow%wb(:,5) - ssnow%wbice(:,6) = 0.90 * ssnow%wb(:,6) + ssnow%wbice(:,1) = frozen_limit * ssnow%wb(:,1) + ssnow%wbice(:,2) = frozen_limit * ssnow%wb(:,2) + ssnow%wbice(:,3) = frozen_limit * ssnow%wb(:,3) + ssnow%wbice(:,4) = frozen_limit * ssnow%wb(:,4) + ssnow%wbice(:,5) = frozen_limit * ssnow%wb(:,5) + ssnow%wbice(:,6) = frozen_limit * ssnow%wb(:,6) ENDWHERE xx=REAL(heat_cap_lower_limit(:,1)) ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil & & + (ssnow%wb(:,1) - ssnow%wbice(:,1) ) * Ccswat * Cdensity_liq & - & + ssnow%wbice(:,1) * Ccsice * Cdensity_liq * .9, xx ) * soil%zse(1) - END IF + & + ssnow%wbice(:,1) * Ccsice * Cdensity_ice, xx ) * soil%zse(1) + ENDIF ! if(.NOT.cable_runtime_coupled) -IF (ktau <= 1) THEN - xx=heat_cap_lower_limit(:,1) - ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil & - & + (ssnow%wb(:,1) - ssnow%wbice(:,1) ) * Ccswat * Cdensity_liq & - & + ssnow%wbice(:,1) * Ccsice * Cdensity_liq * .9, xx ) * soil%zse(1) + & - & (1. - ssnow%isflag) * Ccgsnow * ssnow%snowd -END IF +! this is done in all cases? overwrites gammzz from above and mediates thru +! snowd% (but only for single layer snow) + +! Evaluated on first timestep offline/ ESM1.5 +xx=heat_cap_lower_limit(:,1) +ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil & + & + (ssnow%wb(:,1) - ssnow%wbice(:,1) ) * Ccswat * Cdensity_liq & + & + ssnow%wbice(:,1) * Ccsice * Cdensity_ice, xx ) * soil%zse(1) + & + & (1. - ssnow%isflag) * Ccgsnow * ssnow%snowd END SUBROUTINE spec_init_soil_snow diff --git a/src/science/soilsnow/cbl_soilsnow_main.F90 b/src/science/soilsnow/cbl_soilsnow_main.F90 index b7c100480..c7fd4d7db 100644 --- a/src/science/soilsnow/cbl_soilsnow_main.F90 +++ b/src/science/soilsnow/cbl_soilsnow_main.F90 @@ -40,6 +40,8 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) USE snow_accum_mod, ONLY: snow_accum USE snowdensity_mod, ONLY: snowDensity +USE cable_phys_constants_mod, ONLY: density_liq, density_ice + REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_parameter_type), INTENT(INOUT) :: soil TYPE(soil_snow_type), INTENT(INOUT) :: ssnow @@ -51,21 +53,20 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) REAL, DIMENSION(mp) :: snowmlt REAL, DIMENSION(mp) :: totwet REAL, DIMENSION(mp) :: weting - REAL, DIMENSION(mp) :: xx + REAL(r_2), DIMENSION(mp) :: xx REAL(r_2), DIMENSION(mp) :: xxx REAL(r_2), DIMENSION(mp) :: deltat,sinfil1,sinfil2,sinfil3 REAL :: zsetot INTEGER, SAVE :: ktau =0 -REAL :: wbliq(mp,ms) ktau = ktau +1 !this is the value it is initialized with in cable_common anyway max_glacier_snowd = 1100.0 ! for ACCESS1.3 onwards. = 50000.0 for ACCESS1.0 zsetot = SUM(soil%zse) - ssnow%tggav = 0. + ssnow%tggav(:) = 0. DO k = 1, ms - ssnow%tggav = ssnow%tggav + soil%zse(k)*ssnow%tgg(:,k)/zsetot + ssnow%tggav(:) = ssnow%tggav(:) + ( (soil%zse(k)/zsetot) * ssnow%tgg(:,k) ) soil%heat_cap_lower_limit(:,k) = MAX( 0.01, soil%css(:) * soil%rhosoil(:) ) END DO @@ -83,8 +84,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ssnow%dtmlt = 0.0 ssnow%osnowd = ssnow%snowd - - wbliq = ssnow%wb - ssnow%wbice + ssnow%wbliq = ssnow%wb - ssnow%wbice !%cable_runtime_coupled special initalizations in um_init NA for ESM1.5 @@ -92,7 +92,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) IF (ktau <= 1) & ssnow%gammzz(:,1) = MAX( (1.0 - soil%ssat) * soil%css * soil%rhosoil & & + (ssnow%wb(:,1) - ssnow%wbice(:,1) ) * Ccswat * Cdensity_liq & - & + ssnow%wbice(:,1) * Ccsice * Cdensity_liq * .9, xx ) * soil%zse(1) + & + & + ssnow%wbice(:,1) * Ccsice * Cdensity_ice, xx ) * soil%zse(1) + & & (1. - ssnow%isflag) * Ccgsnow * ssnow%snowd @@ -133,8 +133,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) CALL remove_trans(dels, soil, ssnow, canopy, veg) - CALL soilfreeze(dels, soil, ssnow, soil%heat_cap_lower_limit) - + CALL soilfreeze(dels, soil, ssnow, soil%heat_cap_lower_limit ) totwet = canopy%precis + ssnow%smelt @@ -191,11 +190,13 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ! Set weighted soil/snow surface temperature ssnow%tss=(1-ssnow%isflag)*ssnow%tgg(:,1) + ssnow%isflag*ssnow%tggsn(:,1) - wbliq = ssnow%wb - ssnow%wbice + ssnow%wbliq = ssnow%wb - ssnow%wbice ssnow%wbtot = 0.0 DO k = 1, ms - ssnow%wbtot = ssnow%wbtot + REAL(ssnow%wb(:,k)*1000.0*soil%zse(k),r_2) + ! tot moisture this timestep + ssnow%wbtot(:) = ssnow%wbtot(:) + & + (ssnow%wbliq(:,k)*density_liq + ssnow%wbice(:,k)*density_ice) * soil%zse(k) END DO diff --git a/src/science/soilsnow/cbl_stempv.F90 b/src/science/soilsnow/cbl_stempv.F90 index 2111c0d9a..cd5a280b2 100644 --- a/src/science/soilsnow/cbl_stempv.F90 +++ b/src/science/soilsnow/cbl_stempv.F90 @@ -90,7 +90,7 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) ssnow%gammzz(:,k) = MAX( heat_cap_lower_limit(:,k), & ( 1.0 - soil%ssat ) * soil%css * soil%rhosoil & + soil%ssat * ( wblfsp * Ccswat * Cdensity_liq + & - ssnow%wbfice(:,k) * Ccsice * Cdensity_liq * 0.9 ) ) & + ssnow%wbfice(:,k) * Ccsice * Cdensity_ice ) ) & * soil%zse(k) ssnow%gammzz(:,k) = ssnow%gammzz(:,k) + Ccgsnow * ssnow%snowd @@ -113,7 +113,7 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) ssnow%gammzz(:,k) = MAX( REAL(heat_cap_lower_limit(:,k)), & ( 1.0 - soil%ssat ) * soil%css * soil%rhosoil & + soil%ssat * ( wblfsp * Ccswat * Cdensity_liq + & - ssnow%wbfice(:,k) * Ccsice * Cdensity_liq * 0.9 ) ) & + ssnow%wbfice(:,k) * Ccsice * Cdensity_ice ) ) & * soil%zse(k) dtg = dels / ssnow%gammzz(:,k) @@ -183,8 +183,7 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) ssnow%gammzz(:,k) = MAX( ( 1.0 - soil%ssat ) * soil%css * & soil%rhosoil + soil%ssat * ( wblfsp * Ccswat * & - Cdensity_liq + ssnow%wbfice(:,k) * Ccsice * Cdensity_liq * & - 0.9) , & + Cdensity_liq + ssnow%wbfice(:,k) * Ccsice * Cdensity_ice) , & heat_cap_lower_limit(:,k) ) * soil%zse(k) dtg = dels / ssnow%gammzz(:,k) diff --git a/src/science/soilsnow/cbl_surfbv.F90 b/src/science/soilsnow/cbl_surfbv.F90 index f4182f3d7..b6117bb70 100644 --- a/src/science/soilsnow/cbl_surfbv.F90 +++ b/src/science/soilsnow/cbl_surfbv.F90 @@ -7,12 +7,12 @@ MODULE surfbv_mod CONTAINS SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) -! subrs -USE smoisturev_mod, ONLY: smoisturev +USE smoisturev_mod, ONLY: smoisturev + USE cable_common_module ! data -USE cable_surface_types_mod, ONLY: lakes_cable -USE cable_common_module +USE cable_surface_types_mod, ONLY: lakes_cable +USE cable_other_constants_mod, ONLY : wilt_limitfactor IMPLICIT NONE @@ -51,15 +51,15 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) DO k = 1, ms xxx = REAL( soil%ssat,r_2 ) ssnow%rnof1 = ssnow%rnof1 + REAL( MAX( ssnow%wb(:,k) - xxx, 0.0_r_2 ) & - * 1000.0 ) * soil%zse(k) - ssnow%wb(:,k) = MAX( 0.01_r_2, MIN( ssnow%wb(:,k), xxx ) ) + * Cdensity_liq ) * soil%zse(k) + ssnow%wb(:,k) = MAX( soil%swilt(:)/(2.*wilt_limitfactor), MIN( ssnow%wb(:,k), xxx ) ) END DO ! for deep runoff use wb-sfc, but this value not to exceed .99*wb-wbice ! account for soil/ice cracking ! fracm = MIN(0.2, 1. - MIN(1., ssnow%wb(:,ms) / soil%sfc ) ) ! ssnow%wb(:,ms) = ssnow%wb(:,ms) & - ! + fracm*ssnow%rnof1/(1000.0*soil%zse(ms)) + ! + fracm*ssnow%rnof1/(Cdensity_liq*soil%zse(ms)) ! ssnow%rnof1 = (1. - fracm) * ssnow%rnof1 ! Scaling runoff to kg/m^2/s to match rest of the model @@ -112,13 +112,13 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) ssnow%sinfil = MIN( ssnow%rnof2, ssnow%wb_lake ) ! water that can be extracted from the rnof2 ssnow%rnof2 = MAX( 0.0, ssnow%rnof2 - ssnow%sinfil ) ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil) - xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - REAL(soil%sfc(:),r_2))*soil%zse(ms)*1000.0) + xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - REAL(soil%sfc(:),r_2))*soil%zse(ms)*Cdensity_liq) ssnow%sinfil = MIN( REAL(xxx), ssnow%wb_lake ) - ssnow%wb(:,ms) = ssnow%wb(:,ms) - REAL(ssnow%sinfil / (soil%zse(ms)*1000.0), r_2) + ssnow%wb(:,ms) = ssnow%wb(:,ms) - REAL(ssnow%sinfil / (soil%zse(ms)*Cdensity_liq), r_2) ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil) - xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - 0.5*(soil%sfc + soil%swilt))*soil%zse(ms)*1000.0) + xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - 0.5*(soil%sfc + soil%swilt))*soil%zse(ms)*Cdensity_liq) ssnow%sinfil = MIN( REAL(xxx), ssnow%wb_lake ) - ssnow%wb(:,ms) = ssnow%wb(:,ms) - ssnow%sinfil / (soil%zse(ms)*1000.0) + ssnow%wb(:,ms) = ssnow%wb(:,ms) - ssnow%sinfil / (soil%zse(ms)*Cdensity_liq) ssnow%wb_lake = MAX( 0.0, ssnow%wb_lake - ssnow%sinfil) ENDWHERE