From 6a04a4448ce5bbe7c3e461d76b090945a200099d Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 24 Jun 2025 16:55:15 +1000 Subject: [PATCH 01/28] New calculation for accounting of root distribution on the transpiration Fixes #560 --- src/science/canopy/cbl_dryLeaf.F90 | 30 ++++++++++++++++++++--- src/science/soilsnow/cbl_remove_trans.F90 | 2 +- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index f71f299fe..116348cef 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -145,6 +145,11 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & REAL, DIMENSION(mp,2) :: gsw_term, lower_limit2 ! local temp var INTEGER :: i, j, k, kk ! iteration count + + ! For the calculation of the amount of transpired water + REAL(r_2) :: xxd, xx + REAL(r_2), DIMENSION(0:ms) :: diff + REAL :: vpd, g1 ! Ticket #56 REAL, DIMENSION(mp,mf) :: & xleuning ! leuning stomatal coeff @@ -513,12 +518,29 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & evapfb(i) = ( 1.0 - canopy%fwet(i)) * REAL( ecx(i) ) *dels & / air%rlam(i) + xx = 0.; xxd = 0.; diff(:) = 0. DO kk = 1,ms - ssnow%evapfbl(i,kk) = MIN( evapfb(i) * veg%froot(i,kk), & - MAX( 0.0, REAL( ssnow%wb(i,kk) ) - & - 1.1 * soil%swilt(i) ) * & - soil%zse(kk) * density_liq ) + ! Root water extraction demand + xx = evapfb(i) * veg%froot(i,kk) + diff(kk-1) + ! Maximum water available at this soil layer + diff(kk) = MAX( 0.0_r_2, ssnow%wb(i,kk) - 1.1 * soil%swilt(i)) & ! m3/m3 + * soil%zse(kk)*density_liq + xxd = xx - diff(kk) + + ! ssnow%evapfbl is the water extracted from this layer + ! diff is the excess water demand that is transferred to the next layer + IF (xxd > 0.0) THEN + ssnow%evapfbl(i,kk) = diff(kk) + diff(kk) = xxd + ELSE + ssnow%evapfbl(i,kk) = xx + diff(kk) = 0.0 + END IF + ! ssnow%evapfbl(i,kk) = MIN( evapfb(i) * veg%froot(i,kk), & + ! MAX( 0.0, REAL( ssnow%wb(i,kk) ) - & + ! 1.1 * soil%swilt(i) ) * & + ! soil%zse(kk) * density_liq ) ENDDO IF (cable_user%soil_struc=='default') THEN diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 50a438ad8..f4a49dc98 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -30,7 +30,7 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) ! Calculate the amount (perhaps moisture/ice limited) ! 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 + diff(:,k) = MAX( 0.0_r_2, ssnow%wb(:,k) - 1.1 * soil%swilt) & ! m3/m3 * soil%zse(k)*Cdensity_liq xxd = xx - diff(:,k) From 3334097813876353eab85cec7f7e5761fcbb01aa Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Mon, 23 Jun 2025 16:41:32 +1000 Subject: [PATCH 02/28] (#560) - Remove unused, duplicate canopy%evapfbl. ssnow%evapfbl is the same variable. --- src/offline/cable_define_types.F90 | 3 --- src/offline/cable_mpicommon.F90 | 6 +++--- src/offline/cable_mpimaster.F90 | 21 --------------------- src/offline/cable_mpiworker.F90 | 22 ---------------------- src/science/canopy/cbl_dryLeaf.F90 | 1 - 5 files changed, 3 insertions(+), 50 deletions(-) diff --git a/src/offline/cable_define_types.F90 b/src/offline/cable_define_types.F90 index 9a2699ab7..f6a92a935 100644 --- a/src/offline/cable_define_types.F90 +++ b/src/offline/cable_define_types.F90 @@ -493,7 +493,6 @@ MODULE cable_def_types_mod ga_cor ! correction to ground heat flux (W/m2) REAL, DIMENSION(:,:), POINTER :: & - evapfbl, & gswx, & ! stom cond for water zetar, & ! stability parameter (ref height) ! vh_js ! @@ -1149,7 +1148,6 @@ SUBROUTINE alloc_canopy_type(var, mp) ALLOCATE( var% fwet(mp) ) ALLOCATE( var% fns_cor(mp) ) !REV_CORR variable ALLOCATE( var% ga_cor(mp) ) !REV_CORR variable - ALLOCATE ( var % evapfbl(mp,ms) ) ALLOCATE( var% epot(mp) ) ALLOCATE( var% fnpp(mp) ) ALLOCATE( var% fevw_pot(mp) ) @@ -1783,7 +1781,6 @@ SUBROUTINE dealloc_canopy_type(var) DEALLOCATE( var% fwet ) DEALLOCATE( var% fns_cor ) !REV_CORR variable DEALLOCATE( var% ga_cor ) !REV_CORR variable - DEALLOCATE ( var % evapfbl ) DEALLOCATE( var% epot ) DEALLOCATE( var% fnpp ) DEALLOCATE( var% fevw_pot ) diff --git a/src/offline/cable_mpicommon.F90 b/src/offline/cable_mpicommon.F90 index 15e1c70d2..5070ca57d 100644 --- a/src/offline/cable_mpicommon.F90 +++ b/src/offline/cable_mpicommon.F90 @@ -29,7 +29,7 @@ MODULE cable_mpicommon ! base number of input fields: must correspond to CALLS to ! MPI_address (field ) in *_mpimaster/ *_mpiworker - INTEGER, PARAMETER :: nparam = 341 + INTEGER, PARAMETER :: nparam = 340 ! MPI: extra params sent only if nsoilparmnew is true INTEGER, PARAMETER :: nsoilnew = 1 @@ -77,7 +77,7 @@ MODULE cable_mpicommon !INTEGER, PARAMETER :: nmat = 29 ! MPI: CABLE_r491, after following up with Bernard on the new variables ! vh sli nmat + 4 36 -> 40 - INTEGER, PARAMETER :: nmat = 40 + INTEGER, PARAMETER :: nmat = 39 ! MPI: number of contig vector parts / worker (results) !INTEGER, PARAMETER :: nvec = 149 @@ -103,7 +103,7 @@ MODULE cable_mpicommon ! MPI: number of fields included in restart_t type for data ! that is returned only for creating a restart file at the end of the run ! MPI: gol124: canopy%rwater removed when Bernard ported to CABLE_r491 - INTEGER, PARAMETER :: nrestart = 17 + INTEGER, PARAMETER :: nrestart = 16 INTEGER, PARAMETER :: nsumcasaflux = 62 INTEGER, PARAMETER :: nsumcasapool = 40 INTEGER, PARAMETER :: nclimate = 30 diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 591caba01..74aefa282 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -2508,14 +2508,6 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& ! blen(bidx) = 1 ! !blen(bidx) = ms * r1len - bidx = bidx + 1 - CALL MPI_Get_address (canopy%evapfbl(off,1), displs(bidx), ierr) - ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491 - CALL MPI_Type_create_hvector (ms, r1len, r1stride, MPI_BYTE, & - & types(bidx), ierr) - blen(bidx) = 1 - !blen(bidx) = ms * r2len - bidx = bidx + 1 CALL MPI_Get_address (canopy%epot(off), displs(bidx), ierr) blen(bidx) = r1len @@ -4788,12 +4780,6 @@ SUBROUTINE master_outtypes (comm,met,canopy,ssnow,rad,bal,air,soil,veg) ! REAL(r_2) ! MPI: gol124: backport to r1134 changes r_2 to r_1 ! MPI: gol124: in newest CABLE-cnp it's r_2 again - midx = midx + 1 - CALL MPI_Get_address (canopy%evapfbl(off,1), maddr(midx), ierr) ! 2 - ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491 - CALL MPI_Type_create_hvector (ms, r1len, r1stride, MPI_BYTE, & - & mat_t(midx, rank), ierr) - CALL MPI_Type_commit (mat_t(midx, rank), ierr) midx = midx + 1 CALL MPI_Get_address (canopy%gswx(off,1), maddr(midx), ierr) ! 2 @@ -7222,13 +7208,6 @@ SUBROUTINE master_restart_types (comm, canopy, air, bgc) ! & types(bidx), ierr) ! blocks(bidx) = 1 - bidx = bidx + 1 - CALL MPI_Get_address (canopy%evapfbl(off,1), displs(bidx), ierr) ! 2 - ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491 - CALL MPI_Type_create_hvector (ms, r1len, r1stride, MPI_BYTE, & - & types(bidx), ierr) - blocks(bidx) = 1 - bidx = bidx + 1 CALL MPI_Get_address (bgc%cplant(off,1), displs(bidx), ierr) CALL MPI_Type_create_hvector (ncp, r1len, r1stride, MPI_BYTE, & diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index 030f62af1..f2066beb1 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -1742,11 +1742,6 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& ! CALL MPI_Get_address (canopy%rwater, displs(bidx), ierr) ! blen(bidx) = ms * r1len - bidx = bidx + 1 - CALL MPI_Get_address (canopy%evapfbl, displs(bidx), ierr) - ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491 - blen(bidx) = ms * r1len - bidx = bidx + 1 CALL MPI_Get_address (canopy%epot, displs(bidx), ierr) blen(bidx) = r1len @@ -3679,18 +3674,6 @@ SUBROUTINE worker_outtype (comm,met,canopy,ssnow,rad,bal,air,soil,veg) ! CALL MPI_Get_address (canopy%rwater(off,1), displs(bidx), ierr) ! blocks(bidx) = r1len * ms - ! midx = midx + 1 - ! REAL(r_2) - ! CALL MPI_Get_address (canopy%evapfbl(off,1), maddr(midx), ierr) ! 2 - !CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & - ! & mat_t(midx, rank), ierr) - - ! TODO: skip, used for restart but not output - bidx = bidx + 1 - CALL MPI_Get_address (canopy%evapfbl(off,1), displs(bidx), ierr) - ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491 - blocks(bidx) = r1len * ms - bidx = bidx + 1 CALL MPI_Get_address (canopy%gswx(off,1), displs(bidx), ierr) blocks(bidx) = r1len * mf @@ -6404,11 +6387,6 @@ SUBROUTINE worker_restart_type (comm, canopy, air, bgc) ! CALL MPI_Get_address (canopy%rwater(off,1), displs(bidx), ierr) ! blocks(bidx) = r1len * ms - bidx = bidx + 1 - CALL MPI_Get_address (canopy%evapfbl(off,1), displs(bidx), ierr) - ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491 - blocks(bidx) = r1len * ms - bidx = bidx + 1 CALL MPI_Get_address (bgc%cplant(off,1), displs(bidx), ierr) blocks(bidx) = r1len * ncp diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 116348cef..75cb50a02 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -678,7 +678,6 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & canopy%frday = 12.0 * SUM(rdy, 2) ! vh ! inserted min to avoid -ve values of GPP canopy%fpn = MIN(-12.0 * SUM(an_y, 2), canopy%frday) - canopy%evapfbl = ssnow%evapfbl DEALLOCATE( gswmin ) From 4eacb1cf5c46c15588f38bab21f41720f29443a1 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 24 Jun 2025 14:54:33 +1000 Subject: [PATCH 03/28] (#560) - change ssnow%evapfbl to double precision --- src/offline/cable_define_types.F90 | 2 +- src/science/canopy/cbl_latent_heat.F90 | 50 +++++++++++++++----------- 2 files changed, 31 insertions(+), 21 deletions(-) diff --git a/src/offline/cable_define_types.F90 b/src/offline/cable_define_types.F90 index f6a92a935..5358b6fbe 100644 --- a/src/offline/cable_define_types.F90 +++ b/src/offline/cable_define_types.F90 @@ -278,7 +278,6 @@ MODULE cable_def_types_mod tggsn, & ! snow temperature in K dtmlt, & ! water flux to the soil albsoilsn, & ! soil + snow reflectance - evapfbl, & ! tilefrac ! factor for latent heat @@ -286,6 +285,7 @@ MODULE cable_def_types_mod wbtot ! total soil water (mm) REAL(r_2), DIMENSION(:,:), POINTER :: & + evapfbl, & ! gammzz, & ! heat capacity for each soil layer wb, & ! volumetric soil moisture (solid+liq) wbice, & ! soil ice diff --git a/src/science/canopy/cbl_latent_heat.F90 b/src/science/canopy/cbl_latent_heat.F90 index 124dea5c5..7329c2edc 100644 --- a/src/science/canopy/cbl_latent_heat.F90 +++ b/src/science/canopy/cbl_latent_heat.F90 @@ -60,41 +60,51 @@ SUBROUTINE Latent_heat_flux( mp, CTFRZ, dels, soil_zse, soil_swilt, & INTEGER :: mp REAL(KIND=r_2), INTENT(OUT) :: canopy_fes(mp) -!! latent heat flux from the ground (Wm\(^{-2}\)) + !! latent heat flux from the ground (Wm\(^{-2}\)) REAL(KIND=r_2), INTENT(OUT) :: canopy_fess(mp) -!! latent heat flux from the soil (Wm\(^{-2}\)) + !! latent heat flux from the soil (Wm\(^{-2}\)) REAL(KIND=r_2), INTENT(OUT) :: canopy_fesp(mp) -!! latent heat flux from any puddles (Wm\(^{-2}\)) + !! latent heat flux from any puddles (Wm\(^{-2}\)) REAL, INTENT(OUT) :: ssnow_cls(mp) -!! factor denoting phase of water flux (=1 if liquid, =1.1335 if ice) + !! factor denoting phase of water flux (=1 if liquid, =1.1335 if ice) REAL, INTENT(IN OUT) :: pwet(mp) -!! factor to reduce soil evaporation due to presence of a puddle (-) + !! factor to reduce soil evaporation due to presence of a puddle (-) REAL, INTENT(IN OUT) :: ssnow_wetfac(mp) -!! wetness factor for soil (between 0 and 1) + !! wetness factor for soil (between 0 and 1) -REAL, INTENT(IN) :: CTFRZ !! temperature at freezing point (K) -REAL, INTENT(IN) :: dels !! time step of CABLE (s) -REAL, INTENT(IN) :: soil_zse !! thickness of topmost soil layer (m) +REAL, INTENT(IN) :: CTFRZ + !! temperature at freezing point (K) +REAL, INTENT(IN) :: dels + !! time step of CABLE (s) +REAL, INTENT(IN) :: soil_zse + !! thickness of topmost soil layer (m) REAL, INTENT(IN) :: soil_swilt(mp) -!! soil moisture content at wilting point (m\(^3\) water m\(^{-3}\) volume of soil) + !! soil moisture content at wilting point + !! (m\(^3\) water m\(^{-3}\) volume of soil) LOGICAL , INTENT(IN) :: cable_user_l_new_reduce_soilevp -!! NAMELIST switch to use alternate soil evaporation scheme + !! NAMELIST switch to use alternate soil evaporation scheme -REAL, INTENT(IN) :: air_rlam(mp) !! density of air (kgm\(^{-3}\)) +REAL, INTENT(IN) :: air_rlam(mp) + !! density of air (kgm\(^{-3}\)) REAL, INTENT(IN) :: ssnow_snowd(mp) -!! depth of snow in liquid water equivalent (mm m\(^{-2}\)) -REAL, INTENT(IN) :: ssnow_pudsto(mp) !! amount of water in puddles (kgm\(^{-2}\)) + !! depth of snow in liquid water equivalent (mm m\(^{-2}\)) +REAL, INTENT(IN) :: ssnow_pudsto(mp) + !! amount of water in puddles (kgm\(^{-2}\)) REAL, INTENT(IN) :: ssnow_pudsmx(mp) -!! maximum amount of water possible in puddles (kgm\(^{-2}\)) + !! maximum amount of water possible in puddles (kgm\(^{-2}\)) REAL, INTENT(IN) :: ssnow_potev(mp) -!! latent heat flux associated potential evaporation (Wm\(^{-2}\)) -REAL, INTENT(IN) :: ssnow_evapfbl(mp) !! flux of water from soil surface (kg m\(^{-2}\)) + !! latent heat flux associated potential evaporation (Wm\(^{-2}\)) +REAL(KIND=r_2), INTENT(IN) :: ssnow_evapfbl(mp) + !! flux of water from soil surface (kg m\(^{-2}\)) REAL(KIND=r_2), INTENT(IN) :: ssnow_wb(mp) -!! water content in surface soil layer (m\(^{3}\) liquid water m\(^{-3}\) volume of soil) + !! water content in surface soil layer + !! (m\(^{3}\) liquid water m\(^{-3}\) volume of soil) REAL(KIND=r_2), INTENT(IN) :: ssnow_wbice(mp) -!! ice content in surface soil layer (m\(^{3}\) frozen water m\(^{-3}\) volume of soil) -REAL, INTENT(IN) :: ssnow_tss(mp) !! temperature of surface soil/snow layer (K) + !! ice content in surface soil layer + !! (m\(^{3}\) frozen water m\(^{-3}\) volume of soil) +REAL, INTENT(IN) :: ssnow_tss(mp) + !! temperature of surface soil/snow layer (K) REAL, DIMENSION(mp) :: & frescale, flower_limit, fupper_limit From fcce19cac697c9ddbb046b799671da501d5cc5c0 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 24 Jun 2025 16:01:41 +1000 Subject: [PATCH 04/28] (#560) - Function for the calculation of evapfbl --- src/science/canopy/cbl_dryLeaf.F90 | 109 +++++++++------------- src/science/soilsnow/cbl_remove_trans.F90 | 83 ++++++++++------ 2 files changed, 101 insertions(+), 91 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 75cb50a02..9a0e36891 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -15,33 +15,38 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & USE cable_def_types_mod USE cable_common_module -USE cbl_photosynthesis_module, ONLY : photosynthesis -USE cbl_fwsoil_module, ONLY : fwsoil_calc_std, fwsoil_calc_non_linear, & - fwsoil_calc_Lai_Ktaul, fwsoil_calc_sli -!data -USE cable_surface_types_mod, ONLY: evergreen_broadleaf, deciduous_broadleaf -USE cable_surface_types_mod, ONLY: evergreen_needleleaf, deciduous_needleleaf -USE cable_surface_types_mod, ONLY: c3_grassland, tundra, c3_cropland - -! maths & other constants -USE cable_other_constants_mod, ONLY : CLAI_THRESH => LAI_THRESH -! physical constants -USE cable_phys_constants_mod, ONLY : CTFRZ => TFRZ -USE cable_phys_constants_mod, ONLY : CDHEAT => DHEAT -USE cable_phys_constants_mod, ONLY : CRGAS => RGAS -USE cable_phys_constants_mod, ONLY : CCAPP => CAPP -USE cable_phys_constants_mod, ONLY : CRMAIR => RMAIR -USE cable_phys_constants_mod, ONLY : density_liq -! photosynthetic constants -USE cable_photo_constants_mod, ONLY : CMAXITER => MAXITER ! only integer here -USE cable_photo_constants_mod, ONLY : CTREFK => TREFK -USE cable_photo_constants_mod, ONLY : CGAM0 => GAM0 -USE cable_photo_constants_mod, ONLY : CGAM1 => GAM1 -USE cable_photo_constants_mod, ONLY : CGAM2 => GAM2 -USE cable_photo_constants_mod, ONLY : CRGSWC => RGSWC -USE cable_photo_constants_mod, ONLY : CRGBWC => RGBWC - -IMPLICIT NONE + USE cbl_photosynthesis_module, ONLY : photosynthesis + USE cbl_fwsoil_module, ONLY : fwsoil_calc_std, & + fwsoil_calc_non_linear, & + fwsoil_calc_Lai_Ktaul, & + fwsoil_calc_sli + + USE remove_trans_mod, ONLY : trans_soil_water + + !data + USE cable_surface_types_mod, ONLY: evergreen_broadleaf, deciduous_broadleaf + USE cable_surface_types_mod, ONLY: evergreen_needleleaf, deciduous_needleleaf + USE cable_surface_types_mod, ONLY: c3_grassland, tundra, c3_cropland + + ! maths & other constants + USE cable_other_constants_mod, ONLY : CLAI_THRESH => LAI_THRESH + ! physical constants + USE cable_phys_constants_mod, ONLY : CTFRZ => TFRZ + USE cable_phys_constants_mod, ONLY : CDHEAT => DHEAT + USE cable_phys_constants_mod, ONLY : CRGAS => RGAS + USE cable_phys_constants_mod, ONLY : CCAPP => CAPP + USE cable_phys_constants_mod, ONLY : CRMAIR => RMAIR + USE cable_phys_constants_mod, ONLY : density_liq + ! photosynthetic constants + USE cable_photo_constants_mod, ONLY : CMAXITER => MAXITER ! only integer here + USE cable_photo_constants_mod, ONLY : CTREFK => TREFK + USE cable_photo_constants_mod, ONLY : CGAM0 => GAM0 + USE cable_photo_constants_mod, ONLY : CGAM1 => GAM1 + USE cable_photo_constants_mod, ONLY : CGAM2 => GAM2 + USE cable_photo_constants_mod, ONLY : CRGSWC => RGSWC + USE cable_photo_constants_mod, ONLY : CRGBWC => RGBWC + + IMPLICIT NONE TYPE (radiation_type), INTENT(INOUT) :: rad TYPE (roughness_type), INTENT(INOUT) :: rough @@ -149,6 +154,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & ! For the calculation of the amount of transpired water REAL(r_2) :: xxd, xx REAL(r_2), DIMENSION(0:ms) :: diff + REAL(r_2), DIMENSION(mp) :: local_fevc REAL :: vpd, g1 ! Ticket #56 REAL, DIMENSION(mp,mf) :: & @@ -514,44 +520,19 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & ELSE - IF (ecx(i) > 0.0 .AND. canopy%fwet(i) < 1.0) THEN - evapfb(i) = ( 1.0 - canopy%fwet(i)) * REAL( ecx(i) ) *dels & - / air%rlam(i) - - xx = 0.; xxd = 0.; diff(:) = 0. - DO kk = 1,ms - - ! Root water extraction demand - xx = evapfb(i) * veg%froot(i,kk) + diff(kk-1) - ! Maximum water available at this soil layer - diff(kk) = MAX( 0.0_r_2, ssnow%wb(i,kk) - 1.1 * soil%swilt(i)) & ! m3/m3 - * soil%zse(kk)*density_liq - xxd = xx - diff(kk) - - ! ssnow%evapfbl is the water extracted from this layer - ! diff is the excess water demand that is transferred to the next layer - IF (xxd > 0.0) THEN - ssnow%evapfbl(i,kk) = diff(kk) - diff(kk) = xxd - ELSE - ssnow%evapfbl(i,kk) = xx - diff(kk) = 0.0 - END IF - ! ssnow%evapfbl(i,kk) = MIN( evapfb(i) * veg%froot(i,kk), & - ! MAX( 0.0, REAL( ssnow%wb(i,kk) ) - & - ! 1.1 * soil%swilt(i) ) * & - ! soil%zse(kk) * density_liq ) - - ENDDO - IF (cable_user%soil_struc=='default') THEN - canopy%fevc(i) = SUM(ssnow%evapfbl(i,:))*air%rlam(i)/dels - ecx(i) = canopy%fevc(i) / (1.0-canopy%fwet(i)) - ELSEIF (cable_user%soil_struc=='sli') THEN - canopy%fevc(i) = ecx(i)*(1.0-canopy%fwet(i)) - ENDIF - - ENDIF + local_fevc(i) = ( 1.0 - canopy%fwet(i)) * REAL( ecx(i) ) + IF (local_fevc(i) > 0.0_r_2) THEN + + ssnow%evapfbl(i,:) = trans_soil_water(dels, soil%swilt(i), & + veg%froot(i,:), soil%zse, local_fevc(i), ssnow%wb(i,:)) + IF (cable_user%soil_struc=='default') THEN + canopy%fevc(i) = SUM(ssnow%evapfbl(i,:))*air%rlam(i)/dels + ecx(i) = canopy%fevc(i) / (1.0-canopy%fwet(i)) + ELSE IF (cable_user%soil_struc=='sli') THEN + canopy%fevc(i) = ecx(i)*(1.0-canopy%fwet(i)) + END IF + END IF ENDIF ! Update canopy sensible heat flux: hcx(i) = (sum_rad_rniso(i)-ecx(i) & diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index f4a49dc98..3418a6da2 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -8,7 +8,7 @@ MODULE remove_trans_mod SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) - USE cable_common_module, ONLY : redistrb, cable_user + USE cable_common_module, ONLY : cable_user ! Removes transpiration water from soil. REAL, INTENT(IN) :: dels ! integration time step (s) @@ -16,35 +16,17 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow TYPE(soil_parameter_type), INTENT(INOUT) :: soil TYPE(veg_parameter_type), INTENT(INOUT) :: veg - REAL(r_2), DIMENSION(mp,0:ms) :: diff - REAL(r_2), DIMENSION(mp) :: xx,xxd,evap_cur - INTEGER k + INTEGER i, k IF (cable_user%FWSOIL_switch.NE.'Haverd2013') THEN - xx = 0.; xxd = 0.; diff(:,:) = 0. - DO k = 1,ms - - ! Removing transpiration from soil: - WHERE (canopy%fevc > 0.0 ) ! convert to mm/dels - - ! Calculate the amount (perhaps moisture/ice limited) - ! 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) - 1.1 * soil%swilt) & ! m3/m3 - * 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)*Cdensity_liq) - diff(:,k) = xxd - ELSEWHERE - ssnow%wb(:,k) = ssnow%wb(:,k) - xx / (soil%zse(k)*Cdensity_liq) - diff(:,k) = 0.0 - ENDWHERE - - END WHERE - - END DO + ssnow%evapfbl = 0.0_r_2 + + DO i = 1, mp + ssnow%evapfbl(i,:) = trans_soil_water(dels, soil%swilt(i), & + veg%froot(i,:), soil%zse, canopy%fevc(i), ssnow%wb(i,:)) + ssnow%wb(i,:) = ssnow%wb(i,:) - ssnow%evapfbl(i,:) / (soil%zse(:)*Cdensity_liq) + END DO ELSE WHERE (canopy%fevc .LT. 0.0_r_2) @@ -60,4 +42,51 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) END SUBROUTINE remove_trans + FUNCTION trans_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) + + !! Calculates the amount of water removed from the soil by transpiration. + ! + REAL, INTENT(IN) :: dels + !! integration time step (s) + REAL, INTENT(IN) :: swilt + !! wilting point (m3/m3) + REAL(r_2), INTENT(IN) :: fevc + !! transpiration (kg/m2/s) + REAL, DIMENSION(ms), INTENT(IN) :: froot + !! root fraction (-) + REAL, DIMENSION(ms), INTENT(IN) :: zse + !! soil depth (m) + REAL(r_2), DIMENSION(ms), INTENT(IN) :: wb + !! water balance (m3/m3) + + ! Local variables + REAL(r_2), DIMENSION(ms) :: evapfbl + REAL(r_2), DIMENSION(0:ms) :: diff + REAL(r_2) :: xx,xxd + INTEGER k + + xx = 0.; xxd = 0.; diff(:) = 0. + IF (fevc > 0.0) THEN + DO k = 1,ms + ! Removing transpiration from soil: + + ! Calculate the amount (perhaps moisture/ice limited) + ! that can be removed: + xx = fevc * dels / CHL * froot(k) + diff(k-1) ! kg/m2 + diff(k) = MAX( 0.0_r_2, wb(k) - 1.1 * swilt) * zse(k)*Cdensity_liq + xxd = xx - diff(k) + + IF ( xxd > 0.0 ) THEN + evapfbl(k) = diff(k) + diff(k) = xxd + ELSE + evapfbl(k) = xx + diff(k) = 0.0 + END IF + + END DO + END IF + + END FUNCTION trans_soil_water + END MODULE remove_trans_mod From 332b27284b73e4b048f26bfeedfdb2231005b8c1 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Wed, 25 Jun 2025 12:23:44 +1000 Subject: [PATCH 05/28] (#560) - Update fevc from soilsnow --- src/science/soilsnow/cbl_remove_trans.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 3418a6da2..28f75a5c0 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -23,9 +23,12 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) ssnow%evapfbl = 0.0_r_2 DO i = 1, mp - ssnow%evapfbl(i,:) = trans_soil_water(dels, soil%swilt(i), & - veg%froot(i,:), soil%zse, canopy%fevc(i), ssnow%wb(i,:)) - ssnow%wb(i,:) = ssnow%wb(i,:) - ssnow%evapfbl(i,:) / (soil%zse(:)*Cdensity_liq) + ssnow%evapfbl(i,:) = trans_soil_water(dels, soil%swilt(i), & + veg%froot(i,:), soil%zse, canopy%fevc(i), ssnow%wb(i,:)) + ssnow%wb(i,:) = ssnow%wb(i,:) - ssnow%evapfbl(i,:) / (soil%zse(:)*Cdensity_liq) + IF ( canopy%fevc(i) > 0.0_r_2 ) THEN + canopy%fevc(i) = SUM(ssnow%evapfbl(i,:)) * CHL / dels + END IF END DO ELSE From defcb20b6083b3a82fc56e72e47e60e6cd858ed1 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Fri, 27 Jun 2025 10:09:40 +1000 Subject: [PATCH 06/28] (#560) Clean up and a bit more doc --- src/science/canopy/cbl_dryLeaf.F90 | 17 ++++++++++-- src/science/soilsnow/cbl_remove_trans.F90 | 34 ++++++++--------------- 2 files changed, 26 insertions(+), 25 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 9a0e36891..54b0efaf5 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -21,7 +21,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & fwsoil_calc_Lai_Ktaul, & fwsoil_calc_sli - USE remove_trans_mod, ONLY : trans_soil_water + USE remove_trans_mod, ONLY : transp_soil_water !data USE cable_surface_types_mod, ONLY: evergreen_broadleaf, deciduous_broadleaf @@ -523,7 +523,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & local_fevc(i) = ( 1.0 - canopy%fwet(i)) * REAL( ecx(i) ) IF (local_fevc(i) > 0.0_r_2) THEN - ssnow%evapfbl(i,:) = trans_soil_water(dels, soil%swilt(i), & + ssnow%evapfbl(i,:) = transp_soil_water(dels, soil%swilt(i), & veg%froot(i,:), soil%zse, local_fevc(i), ssnow%wb(i,:)) IF (cable_user%soil_struc=='default') THEN @@ -667,8 +667,19 @@ END SUBROUTINE dryLeaf SUBROUTINE getrex_1d(theta, rex, fws, Fs, thetaS, thetaw, Etrans, gamma, dx, dt, zr) + !! Root extraction : Haverd et al. 2013 + !! **Warning**: This subroutine has diverged from the other `getrex_1d` subroutine + !! in cable_sli_roots.F90. Considering this subroutine predates the other one, + !! it is likely this is an older version and should be updated. Although no, + !! tests has been done to quantify the differences. + !! Changes identified as of 27/06/2025: + !! - `theta` and `thetas` arguments instead ot `thetaS` and `S=theta/thetaS` + !! - `rex` is INTENT(OUT) in sli_roots module. Correct for this version as well. + !! - Condition `WHERE (Fs(:) > zero .AND. layer_depth < zr ) ` changed to + !! `WHERE (Fs(:) > zero` (zr unused in sli_roots) + !! - `IF (ANY(((rex*dt) > MAX((theta(:)-thetaw(:)),zero)*dx(:)) .AND. (Etrans > zero))) THEN` + !! changed to `IF (ANY(((rex*dt) > (theta(:)-thetaw(:))*dx(:)) .AND. ((rex*dt) > zero))) THEN` - ! root extraction : Haverd et al. 2013 USE cable_def_types_mod, ONLY: r_2 IMPLICIT NONE diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 28f75a5c0..f3f603d95 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -7,10 +7,12 @@ MODULE remove_trans_mod CONTAINS SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) + !! Removes transpiration water from soil. + !! For Haverd2013, it also deals with negative canopy + !! transpiration. USE cable_common_module, ONLY : cable_user - ! Removes transpiration water from soil. REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(canopy_type), INTENT(INOUT) :: canopy TYPE(soil_snow_type), INTENT(INOUT) :: ssnow @@ -18,34 +20,22 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) TYPE(veg_parameter_type), INTENT(INOUT) :: veg INTEGER i, k - IF (cable_user%FWSOIL_switch.NE.'Haverd2013') THEN + IF (cable_user%FWSOIL_switch == 'Haverd2013') THEN - ssnow%evapfbl = 0.0_r_2 - - DO i = 1, mp - ssnow%evapfbl(i,:) = trans_soil_water(dels, soil%swilt(i), & - veg%froot(i,:), soil%zse, canopy%fevc(i), ssnow%wb(i,:)) - ssnow%wb(i,:) = ssnow%wb(i,:) - ssnow%evapfbl(i,:) / (soil%zse(:)*Cdensity_liq) - IF ( canopy%fevc(i) > 0.0_r_2 ) THEN - canopy%fevc(i) = SUM(ssnow%evapfbl(i,:)) * CHL / dels - END IF - END DO - - ELSE - WHERE (canopy%fevc .LT. 0.0_r_2) + WHERE (canopy%fevc < 0.0_r_2) canopy%fevw = canopy%fevw+canopy%fevc canopy%fevc = 0.0_r_2 END WHERE - DO k = 1,ms - ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq) - ENDDO - + END IF - ENDIF + DO k = 1,ms + ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq) + END DO END SUBROUTINE remove_trans - FUNCTION trans_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) + + FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) !! Calculates the amount of water removed from the soil by transpiration. ! @@ -90,6 +80,6 @@ FUNCTION trans_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) END DO END IF - END FUNCTION trans_soil_water + END FUNCTION transp_soil_water END MODULE remove_trans_mod From f60d767a5f127b888f18539947c39bfec5bb0541 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Fri, 27 Jun 2025 10:14:38 +1000 Subject: [PATCH 07/28] (#560) - Use liquid soil water only for transpiration --- src/science/canopy/cbl_dryLeaf.F90 | 2 +- src/science/soilsnow/cbl_remove_trans.F90 | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 54b0efaf5..90f9d6b5a 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -524,7 +524,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & IF (local_fevc(i) > 0.0_r_2) THEN ssnow%evapfbl(i,:) = transp_soil_water(dels, soil%swilt(i), & - veg%froot(i,:), soil%zse, local_fevc(i), ssnow%wb(i,:)) + veg%froot(i,:), soil%zse, local_fevc(i), ssnow%wbliq(i,:)) IF (cable_user%soil_struc=='default') THEN canopy%fevc(i) = SUM(ssnow%evapfbl(i,:))*air%rlam(i)/dels diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index f3f603d95..5568676e2 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -29,13 +29,14 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) END IF DO k = 1,ms - ssnow%wb(:,k) = ssnow%wb(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq) + ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq) + ssnow%wb(:,k) = ssnow%wbliq(:,k) + ssnow%wbice(:,k) END DO END SUBROUTINE remove_trans - FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) + FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wbliq) RESULT(evapfbl) !! Calculates the amount of water removed from the soil by transpiration. ! @@ -49,8 +50,8 @@ FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) !! root fraction (-) REAL, DIMENSION(ms), INTENT(IN) :: zse !! soil depth (m) - REAL(r_2), DIMENSION(ms), INTENT(IN) :: wb - !! water balance (m3/m3) + REAL(r_2), DIMENSION(ms), INTENT(IN) :: wbliq + !! liquid soil water (m3/m3) ! Local variables REAL(r_2), DIMENSION(ms) :: evapfbl @@ -66,7 +67,7 @@ FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wb) RESULT(evapfbl) ! Calculate the amount (perhaps moisture/ice limited) ! that can be removed: xx = fevc * dels / CHL * froot(k) + diff(k-1) ! kg/m2 - diff(k) = MAX( 0.0_r_2, wb(k) - 1.1 * swilt) * zse(k)*Cdensity_liq + diff(k) = MAX( 0.0_r_2, wbliq(k) - 1.1 * swilt) * zse(k)*Cdensity_liq xxd = xx - diff(k) IF ( xxd > 0.0 ) THEN From 7bd21f4f2225bac715097d9c57b159f8a1bc9f5d Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Fri, 27 Jun 2025 16:32:09 +1000 Subject: [PATCH 08/28] (#560) - Use swilt_vec and zse_vec for all applications --- src/science/canopy/cbl_dryLeaf.F90 | 4 ++-- src/science/soilsnow/cbl_remove_trans.F90 | 16 ++++++++++++---- src/science/soilsnow/cbl_soilsnow_data.F90 | 2 ++ 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 90f9d6b5a..8185b72d3 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -523,8 +523,8 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & local_fevc(i) = ( 1.0 - canopy%fwet(i)) * REAL( ecx(i) ) IF (local_fevc(i) > 0.0_r_2) THEN - ssnow%evapfbl(i,:) = transp_soil_water(dels, soil%swilt(i), & - veg%froot(i,:), soil%zse, local_fevc(i), ssnow%wbliq(i,:)) + ssnow%evapfbl(i,:) = transp_soil_water(dels, soil%swilt_vec(i,:), & + veg%froot(i,:), soil%zse_vec(i,:), local_fevc(i), ssnow%wbliq(i,:)) IF (cable_user%soil_struc=='default') THEN canopy%fevc(i) = SUM(ssnow%evapfbl(i,:))*air%rlam(i)/dels diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 5568676e2..c521c8317 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -29,10 +29,18 @@ SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) END IF DO k = 1,ms - ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/(soil%zse(k)*Cdensity_liq) + ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/ & + (soil%zse_vec(:,k)*Cdensity_liq) ssnow%wb(:,k) = ssnow%wbliq(:,k) + ssnow%wbice(:,k) END DO + IF (cable_user%gw_model) THEN + ssnow%wb = ssnow%wbliq + den_rat * ssnow%wbice + ssnow%wmliq = ssnow%wbliq * zse_vec * Cdensity_liq !mass + ssnow%wmtot = ssnow%wmliq + ssnow%wmice !mass + + END IF + END SUBROUTINE remove_trans @@ -42,13 +50,13 @@ FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wbliq) RESULT(evapfbl) ! REAL, INTENT(IN) :: dels !! integration time step (s) - REAL, INTENT(IN) :: swilt + REAL(r_2), DIMENSION(ms), INTENT(IN) :: swilt !! wilting point (m3/m3) REAL(r_2), INTENT(IN) :: fevc !! transpiration (kg/m2/s) REAL, DIMENSION(ms), INTENT(IN) :: froot !! root fraction (-) - REAL, DIMENSION(ms), INTENT(IN) :: zse + REAL(r_2), DIMENSION(ms), INTENT(IN) :: zse !! soil depth (m) REAL(r_2), DIMENSION(ms), INTENT(IN) :: wbliq !! liquid soil water (m3/m3) @@ -67,7 +75,7 @@ FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wbliq) RESULT(evapfbl) ! Calculate the amount (perhaps moisture/ice limited) ! that can be removed: xx = fevc * dels / CHL * froot(k) + diff(k-1) ! kg/m2 - diff(k) = MAX( 0.0_r_2, wbliq(k) - 1.1 * swilt) * zse(k)*Cdensity_liq + diff(k) = MAX( 0.0_r_2, wbliq(k) - 1.1 * swilt(k)) * zse(k)*Cdensity_liq xxd = xx - diff(k) IF ( xxd > 0.0 ) THEN diff --git a/src/science/soilsnow/cbl_soilsnow_data.F90 b/src/science/soilsnow/cbl_soilsnow_data.F90 index 88ac8adc5..738432c4b 100644 --- a/src/science/soilsnow/cbl_soilsnow_data.F90 +++ b/src/science/soilsnow/cbl_soilsnow_data.F90 @@ -21,4 +21,6 @@ MODULE cbl_ssnow_data_mod IMPLICIT NONE +REAL, PARAMETER :: den_rat = Cdensity_ice / Cdensity_liq + END MODULE cbl_ssnow_data_mod From 43c8ffe294a37f33a4de5401555fda658d0e5134 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Fri, 27 Jun 2025 16:32:25 +1000 Subject: [PATCH 09/28] (#560) - Correct initialisation of zse_vec --- src/offline/cable_parameters.F90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 0c553c89b..44ebc6b9d 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -1405,9 +1405,6 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & soil%watr(landpt(e)%cstart:landpt(e)%cend,klev) = & REAL(inWatr(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%zse_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(soil%zse(landpt(e)%cstart:landpt(e)%cend), r_2) - soil%css_vec(landpt(e)%cstart:landpt(e)%cend, klev) = & REAL(incss(landpt(e)%ilon, landpt(e)%ilat), r_2) @@ -1684,12 +1681,13 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & soil%GWdz = MAX(1.0,MIN(20.0,soil%GWdz - SUM(soil%zse,dim=1))) !set vectorized versions as same as defaut for now - soil%swilt_vec(:,:) = REAL(SPREAD(soil%swilt(:),2,ms),r_2) - soil%sfc_vec(:,:) = REAL(SPREAD(soil%sfc(:),2,ms),r_2) + soil%swilt_vec(:,:) = REAL(SPREAD(soil%swilt(:),2,ms),r_2) + soil%sfc_vec(:,:) = REAL(SPREAD(soil%sfc(:),2,ms),r_2) soil%sucs_vec(:,:) = REAL(SPREAD(soil%sucs(:),2,ms),r_2) - soil%bch_vec(:,:) = REAL(SPREAD(soil%bch(:),2,ms),r_2) + soil%bch_vec(:,:) = REAL(SPREAD(soil%bch(:),2,ms),r_2) soil%ssat_vec(:,:) = REAL(SPREAD(soil%ssat(:),2,ms),r_2) soil%hyds_vec(:,:) = REAL(SPREAD(soil%hyds(:),2,ms),r_2) + soil%zse_vec(:,:) = REAL(SPREAD(soil%zse(:),1,mp),r_2) END SUBROUTINE write_default_params !============================================================================= From 4d0093fc694847c6d199e19e88632f16066ca129 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Mon, 30 Jun 2025 15:37:03 +1000 Subject: [PATCH 10/28] (#560) - Doc formatting --- src/science/canopy/cbl_dryLeaf.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 8185b72d3..8a07acd8f 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -668,16 +668,19 @@ END SUBROUTINE dryLeaf SUBROUTINE getrex_1d(theta, rex, fws, Fs, thetaS, thetaw, Etrans, gamma, dx, dt, zr) !! Root extraction : Haverd et al. 2013 + !! !! **Warning**: This subroutine has diverged from the other `getrex_1d` subroutine !! in cable_sli_roots.F90. Considering this subroutine predates the other one, - !! it is likely this is an older version and should be updated. Although no, - !! tests has been done to quantify the differences. + !! it is likely this is an older version and should be updated. No + !! tests have been done to quantify the differences. + !! !! Changes identified as of 27/06/2025: - !! - `theta` and `thetas` arguments instead ot `thetaS` and `S=theta/thetaS` - !! - `rex` is INTENT(OUT) in sli_roots module. Correct for this version as well. - !! - Condition `WHERE (Fs(:) > zero .AND. layer_depth < zr ) ` changed to + !! + !! - `theta` and `thetas` arguments instead ot `thetaS` and `S=theta/thetaS` + !! - `rex` is INTENT(OUT) in sli_roots module. Correct for this version as well. + !! - Condition `WHERE (Fs(:) > zero .AND. layer_depth < zr ) ` changed to !! `WHERE (Fs(:) > zero` (zr unused in sli_roots) - !! - `IF (ANY(((rex*dt) > MAX((theta(:)-thetaw(:)),zero)*dx(:)) .AND. (Etrans > zero))) THEN` + !! - `IF (ANY(((rex*dt) > MAX((theta(:)-thetaw(:)),zero)*dx(:)) .AND. (Etrans > zero))) THEN` !! changed to `IF (ANY(((rex*dt) > (theta(:)-thetaw(:))*dx(:)) .AND. ((rex*dt) > zero))) THEN` USE cable_def_types_mod, ONLY: r_2 From 4c1c9244b31cd0b220878694a81f66725471634b Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Mon, 30 Jun 2025 16:01:19 +1000 Subject: [PATCH 11/28] (#560) - Removed unused variables and arguments --- src/science/canopy/cbl_dryLeaf.F90 | 2 -- src/science/soilsnow/cbl_remove_trans.F90 | 4 +--- src/science/soilsnow/cbl_soilsnow_main.F90 | 2 +- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index 8a07acd8f..c3a389459 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -152,8 +152,6 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & INTEGER :: i, j, k, kk ! iteration count ! For the calculation of the amount of transpired water - REAL(r_2) :: xxd, xx - REAL(r_2), DIMENSION(0:ms) :: diff REAL(r_2), DIMENSION(mp) :: local_fevc REAL :: vpd, g1 ! Ticket #56 diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index c521c8317..d0fe90aab 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -6,18 +6,16 @@ MODULE remove_trans_mod CONTAINS -SUBROUTINE remove_trans(dels, soil, ssnow, canopy, veg) +SUBROUTINE remove_trans(soil, ssnow, canopy) !! Removes transpiration water from soil. !! For Haverd2013, it also deals with negative canopy !! transpiration. USE cable_common_module, ONLY : cable_user - REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(canopy_type), INTENT(INOUT) :: canopy TYPE(soil_snow_type), INTENT(INOUT) :: ssnow TYPE(soil_parameter_type), INTENT(INOUT) :: soil - TYPE(veg_parameter_type), INTENT(INOUT) :: veg INTEGER i, k IF (cable_user%FWSOIL_switch == 'Haverd2013') THEN diff --git a/src/science/soilsnow/cbl_soilsnow_main.F90 b/src/science/soilsnow/cbl_soilsnow_main.F90 index c7fd4d7db..1751a5982 100644 --- a/src/science/soilsnow/cbl_soilsnow_main.F90 +++ b/src/science/soilsnow/cbl_soilsnow_main.F90 @@ -131,7 +131,7 @@ SUBROUTINE soil_snow(dels, soil, ssnow, canopy, met, bal, veg) ! Add new snow melt to global snow melt variable: ssnow%smelt = ssnow%smelt + snowmlt - CALL remove_trans(dels, soil, ssnow, canopy, veg) + CALL remove_trans(soil, ssnow, canopy) CALL soilfreeze(dels, soil, ssnow, soil%heat_cap_lower_limit ) From 55200244af845538c9481c336dffb87e83da323b Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 1 Jul 2025 11:00:58 +1000 Subject: [PATCH 12/28] (#560) - Use remove_trans for GW soil --- src/science/gw_hydro/cable_gw_hydro.F90 | 2 +- src/science/soilsnow/cbl_remove_trans.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index 8fb9f90ac..429db06b1 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -1090,7 +1090,7 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) 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) + CALL remove_trans(soil, ssnow, canopy) !> 12. Snow freezes and melts. CALL GWsoilfreeze(dels, soil, ssnow) diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index d0fe90aab..d7ed0cf3c 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -16,7 +16,7 @@ SUBROUTINE remove_trans(soil, ssnow, canopy) TYPE(canopy_type), INTENT(INOUT) :: canopy TYPE(soil_snow_type), INTENT(INOUT) :: ssnow TYPE(soil_parameter_type), INTENT(INOUT) :: soil - INTEGER i, k + INTEGER k IF (cable_user%FWSOIL_switch == 'Haverd2013') THEN @@ -34,7 +34,7 @@ SUBROUTINE remove_trans(soil, ssnow, canopy) IF (cable_user%gw_model) THEN ssnow%wb = ssnow%wbliq + den_rat * ssnow%wbice - ssnow%wmliq = ssnow%wbliq * zse_vec * Cdensity_liq !mass + ssnow%wmliq = ssnow%wbliq * soil%zse_vec * Cdensity_liq !mass ssnow%wmtot = ssnow%wmliq + ssnow%wmice !mass END IF From f097c1e66b58f153911ad8915889761fc878b4b0 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 1 Jul 2025 13:03:52 +1000 Subject: [PATCH 13/28] (#560) - Fix MPI for evapfbl --- src/offline/cable_mpimaster.F90 | 4 ++-- src/offline/cable_mpiworker.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 74aefa282..15b9ff4bc 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -1930,7 +1930,7 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& bidx = bidx + 1 CALL MPI_Get_address (ssnow%evapfbl(off,1), displs(bidx), ierr) - CALL MPI_Type_create_hvector (ms, r1len, r1stride, MPI_BYTE, & + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & & types(bidx), ierr) blen(bidx) = 1 @@ -4864,7 +4864,7 @@ SUBROUTINE master_outtypes (comm,met,canopy,ssnow,rad,bal,air,soil,veg) midx = midx + 1 ! REAL(r_1) CALL MPI_Get_address (ssnow%evapfbl(off,1), maddr(midx), ierr) ! 12 - CALL MPI_Type_create_hvector (ms, r1len, r1stride, MPI_BYTE, & + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & & mat_t(midx, rank), ierr) CALL MPI_Type_commit (mat_t(midx, rank), ierr) midx = midx + 1 diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index f2066beb1..8358b21f8 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -1198,7 +1198,7 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& bidx = bidx + 1 CALL MPI_Get_address (ssnow%evapfbl, displs(bidx), ierr) - blen(bidx) = ms * r1len + blen(bidx) = ms * r2len bidx = bidx + 1 CALL MPI_Get_address (ssnow%qstss, displs(bidx), ierr) @@ -3780,7 +3780,7 @@ SUBROUTINE worker_outtype (comm,met,canopy,ssnow,rad,bal,air,soil,veg) bidx = bidx + 1 CALL MPI_Get_address (ssnow%evapfbl(off,1), displs(bidx), ierr) - blocks(bidx) = r1len * ms + blocks(bidx) = r2len * ms !midx = midx + 1 ! REAL(r_1) From e752a90a6c339f084cdd7c2e670987532c40e582 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 1 Jul 2025 16:30:29 +1000 Subject: [PATCH 14/28] (#560) - Change fwsoil calculations to use wbliq and *_vec soil properties --- src/science/canopy/cbl_fwsoil.F90 | 36 ++++++++++++------------------- 1 file changed, 14 insertions(+), 22 deletions(-) diff --git a/src/science/canopy/cbl_fwsoil.F90 b/src/science/canopy/cbl_fwsoil.F90 index 4e25e0854..f83608bf3 100644 --- a/src/science/canopy/cbl_fwsoil.F90 +++ b/src/science/canopy/cbl_fwsoil.F90 @@ -19,22 +19,13 @@ SUBROUTINE fwsoil_calc_std(fwsoil, soil, ssnow, veg) REAL, INTENT(OUT), DIMENSION(:):: fwsoil ! soil water modifier of stom. cond REAL, DIMENSION(mp) :: rwater ! soil water availability - !note even though swilt_vec is defined in default model it is r_2 - !and even using real(_vec) gives results different from trunk (rounding - !errors) + ! Moving to use *_vec variables even outside of the groundwater + ! option for simplicity. It introduces small rounding differences. - IF (.NOT.cable_user%gw_model) THEN + rwater = MAX(1.0e-9, & + SUM(veg%froot * MAX(1.0e-9,MIN(1.0, REAL((ssnow%wbliq - & + soil%swilt_vec)/(soil%sfc_vec-soil%swilt_vec)) )),2) ) - rwater = MAX(1.0e-9, & - SUM(veg%froot * MAX(1.0e-9,MIN(1.0, REAL(ssnow%wb) - & - SPREAD(soil%swilt, 2, ms))),2) /(soil%sfc-soil%swilt)) - - ELSE - rwater = MAX(1.0e-9, & - SUM(veg%froot * MAX(1.0e-9,MIN(1.0, REAL((ssnow%wbliq - & - soil%swilt_vec)/(soil%sfc_vec-soil%swilt_vec)) )),2) ) - - ENDIF ! Remove vbeta #56 IF(cable_user%GS_SWITCH == 'medlyn') THEN @@ -59,8 +50,8 @@ SUBROUTINE fwsoil_calc_non_linear(fwsoil, soil, ssnow, veg) INTEGER :: j rwater = MAX(1.0e-9, & - SUM(veg%froot * MAX(0.0,MIN(1.0, REAL(ssnow%wb) - & - SPREAD(soil%swilt, 2, ms))),2) /(soil%sfc-soil%swilt)) + SUM(veg%froot * MAX(0.0,MIN(1.0, REAL((ssnow%wbliq) - & + SPREAD(soil%swilt, 2, ms)))),2) /(soil%sfc-soil%swilt)) fwsoil = 1. @@ -112,14 +103,15 @@ SUBROUTINE fwsoil_calc_Lai_Ktaul(fwsoil, soil, ssnow, veg) DO ns=1,ms - dummy(:) = rootgamma/MAX(1.0e-3_r_2,ssnow%wb(:,ns)-soil%swilt(:)) + dummy(:) = rootgamma/MAX(1.0e-3_r_2,ssnow%wbliq(:,ns)-soil%swilt_vec(:,ns)) - frwater(:,ns) = MAX(1.0e-4_r_2,((ssnow%wb(:,ns)-soil%swilt(:))/soil%ssat(:)) & - ** dummy) + frwater(:,ns) = MAX(1.0e-4_r_2,((ssnow%wbliq(:,ns)-soil%swilt_vec(:,ns))/ & + soil%ssat_vec(:,ns))** dummy) fwsoil(:) = MIN(1.0,MAX(fwsoil(:),frwater(:,ns))) - normFac(:) = normFac(:) + frwater(:,ns) * veg%froot(:,ns) + ! normFac unused. + !normFac(:) = normFac(:) + frwater(:,ns) * veg%froot(:,ns) ENDDO @@ -134,8 +126,8 @@ SUBROUTINE fwsoil_calc_sli(fwsoil, soil, ssnow, veg) REAL, INTENT(OUT), DIMENSION(:):: fwsoil ! soil water modifier of stom. cond REAL, DIMENSION(mp,ms):: tmp2d1, tmp2d2, delta_root, alpha2a_root, alpha2_root ! Lai and Katul formulation for root efficiency function vh 17/07/09 - alpha2a_root = MAX(ssnow%wb-soil%swilt_vec, 0.001_r_2)/(soil%ssat_vec) - tmp2d1 = ssnow%wb -soil%swilt_vec + alpha2a_root = MAX(ssnow%wbliq-soil%swilt_vec, 0.001_r_2)/(soil%ssat_vec) + tmp2d1 = ssnow%wbliq -soil%swilt_vec tmp2d2 = SPREAD(veg%gamma,2,ms)/tmp2d1*LOG(alpha2a_root) WHERE ((tmp2d1>0.001) .AND. (tmp2d2 > -10.0)) alpha2_root = EXP(tmp2d2) From 71447e5d340f9bada81901401f8331ebcad1f2fc Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Wed, 2 Jul 2025 12:14:22 +1000 Subject: [PATCH 15/28] (#560) - Attribute negative fevc to fevw in all cases --- src/science/soilsnow/cbl_remove_trans.F90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index d7ed0cf3c..8096a6446 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -8,8 +8,8 @@ MODULE remove_trans_mod SUBROUTINE remove_trans(soil, ssnow, canopy) !! Removes transpiration water from soil. - !! For Haverd2013, it also deals with negative canopy - !! transpiration. + !! We also attribute the negative canopy transpiration (dew) + !! to the wet canopy flux. USE cable_common_module, ONLY : cable_user @@ -18,13 +18,11 @@ SUBROUTINE remove_trans(soil, ssnow, canopy) TYPE(soil_parameter_type), INTENT(INOUT) :: soil INTEGER k - IF (cable_user%FWSOIL_switch == 'Haverd2013') THEN - WHERE (canopy%fevc < 0.0_r_2) - canopy%fevw = canopy%fevw+canopy%fevc - canopy%fevc = 0.0_r_2 - END WHERE - END IF + WHERE (canopy%fevc < 0.0_r_2) + canopy%fevw = canopy%fevw+canopy%fevc + canopy%fevc = 0.0_r_2 + END WHERE DO k = 1,ms ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/ & From f6b58b33db9ccc1b5ead9dd1435b37767ea19893 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Wed, 2 Jul 2025 12:25:09 +1000 Subject: [PATCH 16/28] (#560) - improve comments --- src/science/soilsnow/cbl_remove_trans.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/science/soilsnow/cbl_remove_trans.F90 b/src/science/soilsnow/cbl_remove_trans.F90 index 8096a6446..d7ce9f6a7 100644 --- a/src/science/soilsnow/cbl_remove_trans.F90 +++ b/src/science/soilsnow/cbl_remove_trans.F90 @@ -66,10 +66,15 @@ FUNCTION transp_soil_water(dels, swilt, froot, zse, fevc, wbliq) RESULT(evapfbl) xx = 0.; xxd = 0.; diff(:) = 0. IF (fevc > 0.0) THEN DO k = 1,ms - ! Removing transpiration from soil: ! Calculate the amount (perhaps moisture/ice limited) ! that can be removed: + ! xx: water demand from the transpiration and above soil layers + ! diff(k-1): excess demand from higher soil layers + ! diff(k): maximum water amount available for this layer (supply) + ! xxd: demand minus supply. If the demand is larger (xxd>0), + ! evapfbl is limited by the supply and the excess demand is shifted + ! to the next layer. xx = fevc * dels / CHL * froot(k) + diff(k-1) ! kg/m2 diff(k) = MAX( 0.0_r_2, wbliq(k) - 1.1 * swilt(k)) * zse(k)*Cdensity_liq xxd = xx - diff(k) From 24ef0425ad27ded453ebe7463be8ed1708559d2a Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Wed, 2 Jul 2025 14:20:37 +1000 Subject: [PATCH 17/28] (#560) - fix broken link in doc. New domain for ACCESS-Hive docs --- .../docs/developer_guide/documentation_guidelines/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/docs/developer_guide/documentation_guidelines/index.md b/documentation/docs/developer_guide/documentation_guidelines/index.md index f960508e3..296a95c0a 100644 --- a/documentation/docs/developer_guide/documentation_guidelines/index.md +++ b/documentation/docs/developer_guide/documentation_guidelines/index.md @@ -29,5 +29,5 @@ To help you find the file corresponding to a page, on the rendered [documentatio [cheat-sheets]: ../other_resources/cheat_sheets.md [api-guidelines]: science_doc.md [doc-pages]: https://cable.readthedocs.io/en/latest -[Hive-contribute]: https://access-hive.org.au/about/contribute/ +[Hive-contribute]: https://docs.access-hive.org.au/about/contribute/ [cable-lsm-join]: https://github.com/CABLE-LSM/CABLE/issues/110 From 7e83ba75035e0f64d819a3040a15260dd5c39488 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Mon, 7 Jul 2025 09:27:48 +1000 Subject: [PATCH 18/28] (#560) - Remove remove_transGW, unused now. --- src/science/gw_hydro/cable_gw_hydro.F90 | 82 ------------------------- 1 file changed, 82 deletions(-) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index 429db06b1..7bc56a0d8 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -244,88 +244,6 @@ 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 - - ! Removes transpiration water from soil. - REAL, INTENT(IN) :: dels ! integration time step (s) - TYPE(canopy_type), INTENT(INOUT) :: canopy - TYPE(soil_snow_type), INTENT(INOUT) :: ssnow - 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 - REAL(r_2), DIMENSION(mp,ms) :: zse_mp_mm - INTEGER :: k,i - - DO k=1,ms - DO i=1,mp - zse_mp_mm(i,k) = REAL(soil%zse_vec(i,k)*Cdensity_liq,r_2) - END DO - END DO - - IF (cable_user%FWSOIL_switch.NE.'Haverd2013') THEN - - xx(:) = 0._r_2 - xxd(:) = 0._r_2 - diff(:,:) = 0._r_2 - - DO k = 1,ms - - DO i=1,mp - - IF (canopy%fevc(i) .GT. 0._r_2) THEN - - 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) - - IF (xxd(i) .GT. 0._r_2) THEN - ssnow%wbliq(i,k) = ssnow%wbliq(i,k) - diff(i,k)/zse_mp_mm(i,k) - diff(i,k) = xxd(i) - ELSE - ssnow%wbliq(i,k) = ssnow%wbliq(i,k) - xx(i)/zse_mp_mm(i,k) - diff(i,k) = 0._r_2 - END IF - - - END IF !fvec > 0 - - END DO !mp - END DO !ms - - ELSE - - WHERE (canopy%fevc .LT. 0.0_r_2) - canopy%fevw = canopy%fevw+canopy%fevc - canopy%fevc = 0.0_r_2 - END WHERE - DO k = 1,ms - ssnow%wbliq(:,k) = ssnow%wbliq(:,k) - ssnow%evapfbl(:,k)/(soil%zse_vec(:,k)*m2mm) - ENDDO - - ENDIF - - DO k=1,ms - 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) + den_rat * ssnow%wbice(i,k) !volume ! MMY - END DO - END DO - - - END SUBROUTINE remove_transGW - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!MD GW code from here on!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !---------------------------------------------------------------------- From 84ed8197220632c87febbdf6b44e6ce844898f38 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Wed, 9 Jul 2025 11:36:00 +1000 Subject: [PATCH 19/28] (#560) - Add use statement for remove_trans in gw_hydro. --- src/science/gw_hydro/cable_gw_hydro.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/science/gw_hydro/cable_gw_hydro.F90 b/src/science/gw_hydro/cable_gw_hydro.F90 index 7bc56a0d8..dd14fadd7 100644 --- a/src/science/gw_hydro/cable_gw_hydro.F90 +++ b/src/science/gw_hydro/cable_gw_hydro.F90 @@ -840,6 +840,7 @@ SUBROUTINE soil_snow_gw(dels, soil, ssnow, canopy, met, bal, veg) USE cable_common_module USE snow_processes_soil_thermal_mod, ONLY : snow_processes_soil_thermal ! inserted by rk4417 - phase2 + USE remove_trans_mod, ONLY : remove_trans REAL , INTENT(IN) :: dels ! integration time step (s) TYPE(soil_parameter_type), INTENT(INOUT) :: soil From aa9c084692adeb9169f7960cf3b6549ce97451ca Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Wed, 9 Jul 2025 15:41:11 +1000 Subject: [PATCH 20/28] Update comment for evapfbl to reflect change to r_2. Co-authored-by: Sean Bryan <39685865+SeanBryan51@users.noreply.github.com> --- src/offline/cable_mpimaster.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 15b9ff4bc..a4529f9b7 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -4862,7 +4862,7 @@ SUBROUTINE master_outtypes (comm,met,canopy,ssnow,rad,bal,air,soil,veg) & mat_t(midx, rank), ierr) CALL MPI_Type_commit (mat_t(midx, rank), ierr) midx = midx + 1 - ! REAL(r_1) + ! REAL(r_2) CALL MPI_Get_address (ssnow%evapfbl(off,1), maddr(midx), ierr) ! 12 CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & & mat_t(midx, rank), ierr) From d6acbe3ac510ff95898b42a118304d8ece85403b Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Thu, 10 Jul 2025 16:15:46 +1000 Subject: [PATCH 21/28] (#560) - Add variable definition changes in coupled --- src/coupled/AM3/control/cable/CM3/cbl_canopy_type.F90 | 7 ------- src/coupled/AM3/control/cable/CM3/ssnow_type.F90 | 4 ++-- src/coupled/shared/cable_canopy_type_mod.F90 | 2 -- src/coupled/shared/cable_soilsnow_type_mod.F90 | 2 +- 4 files changed, 3 insertions(+), 12 deletions(-) diff --git a/src/coupled/AM3/control/cable/CM3/cbl_canopy_type.F90 b/src/coupled/AM3/control/cable/CM3/cbl_canopy_type.F90 index 615642fa7..9a310d437 100644 --- a/src/coupled/AM3/control/cable/CM3/cbl_canopy_type.F90 +++ b/src/coupled/AM3/control/cable/CM3/cbl_canopy_type.F90 @@ -63,7 +63,6 @@ MODULE cable_canopy_type_mod REAL, ALLOCATABLE, PUBLIC :: fwet(:) ! fraction of canopy wet REAL, ALLOCATABLE, PUBLIC :: fns_cor(:) ! correction to net rad avail to soil (W/m2) REAL, ALLOCATABLE, PUBLIC :: ga_cor(:) ! correction to ground heat flux (W/m2) - REAL, ALLOCATABLE, PUBLIC :: evapfbl(:,:) ! REAL, ALLOCATABLE, PUBLIC :: gswx(:,:) ! stom cond for water REAL, ALLOCATABLE, PUBLIC :: zetar(:,:) ! stability parameter (ref height) REAL, ALLOCATABLE, PUBLIC :: zetash(:,:) ! stability parameter (shear height) @@ -144,7 +143,6 @@ MODULE cable_canopy_type_mod REAL, POINTER, PUBLIC :: fwet(:) ! fraction of canopy wet REAL, POINTER, PUBLIC :: fns_cor(:) ! correction to net rad avail to soil (W/m2) REAL, POINTER, PUBLIC :: ga_cor(:) ! correction to ground heat flux (W/m2) - REAL, POINTER, PUBLIC :: evapfbl(:,:) ! REAL, POINTER, PUBLIC :: gswx(:,:) ! stom cond for water REAL, POINTER, PUBLIC :: zetar(:,:) ! stability parameter (ref height) REAL, POINTER, PUBLIC :: zetash(:,:) ! stability parameter (shear height) @@ -230,7 +228,6 @@ SUBROUTINE alloc_canopy_type(var, mp) ALLOCATE( var% fwet(mp) ) ALLOCATE( var% fns_cor(mp) ) !REV_CORR variable ALLOCATE( var% ga_cor(mp) ) !REV_CORR variable -ALLOCATE( var % evapfbl(mp,nsl) ) ALLOCATE( var% epot(mp) ) ALLOCATE( var% fnpp(mp) ) ALLOCATE( var% fevw_pot(mp) ) @@ -308,7 +305,6 @@ SUBROUTINE alloc_canopy_type(var, mp) var % fwet(:) = 0.0 var % fns_cor(:) = 0.0 var % ga_cor(:) = 0.0 -var % evapfbl(:,:) = 0.0 var % gswx(:,:) = 0.0 var % zetar(:,:) = 0.0 var % zetash(:,:) = 0.0 @@ -382,7 +378,6 @@ SUBROUTINE dealloc_canopy_type(var) DEALLOCATE( var% fwet ) DEALLOCATE( var% fns_cor ) !REV_CORR variable DEALLOCATE( var% ga_cor ) !REV_CORR variable - DEALLOCATE ( var % evapfbl ) DEALLOCATE( var% epot ) DEALLOCATE( var% fnpp ) DEALLOCATE( var% fevw_pot ) @@ -473,7 +468,6 @@ SUBROUTINE assoc_canopy_type(canopy, canopy_data ) canopy% fwet => canopy_data% fwet canopy% fns_cor => canopy_data% fns_cor canopy% ga_cor => canopy_data% ga_cor -canopy% evapfbl => canopy_data% evapfbl canopy% gswx => canopy_data% gswx canopy% zetar => canopy_data% zetar canopy% zetash => canopy_data% zetash @@ -557,7 +551,6 @@ SUBROUTINE nullify_canopy_cbl( var ) NULLIFY( var% fwet ) NULLIFY( var% fns_cor ) !REV_CORR variable NULLIFY( var% ga_cor ) !REV_CORR variable - NULLIFY( var % evapfbl ) NULLIFY( var% epot ) NULLIFY( var% fnpp ) NULLIFY( var% fevw_pot ) diff --git a/src/coupled/AM3/control/cable/CM3/ssnow_type.F90 b/src/coupled/AM3/control/cable/CM3/ssnow_type.F90 index c107987a6..de57ec639 100644 --- a/src/coupled/AM3/control/cable/CM3/ssnow_type.F90 +++ b/src/coupled/AM3/control/cable/CM3/ssnow_type.F90 @@ -74,11 +74,11 @@ MODULE cable_soil_snow_type_mod REAL, ALLOCATABLE :: tggsn (:,:) ! snow temperature in K REAL, ALLOCATABLE :: dtmlt (:,:) ! water flux to the soil REAL, ALLOCATABLE :: albsoilsn (:,:) ! soil + snow reflectance - REAL, ALLOCATABLE :: evapfbl (:,:) ! REAL, ALLOCATABLE :: tilefrac (:,:) ! factor for latent heat REAL(r_2), ALLOCATABLE :: wbtot (:) ! total soil water (mm) + REAL(r_2), ALLOCATABLE :: evapfbl (:,:) ! REAL(r_2), ALLOCATABLE :: gammzz (:,:) ! heat capacity for each soil layer REAL(r_2), ALLOCATABLE :: wb (:,:) ! volumetric soil moisture (solid+liq) REAL(r_2), ALLOCATABLE :: wbice (:,:) ! soil ice @@ -209,11 +209,11 @@ MODULE cable_soil_snow_type_mod REAL, POINTER :: tggsn (:,:) ! snow temperature in K REAL, POINTER :: dtmlt (:,:) ! water flux to the soil REAL, POINTER :: albsoilsn (:,:) ! soil + snow reflectance - REAL, POINTER :: evapfbl (:,:) ! REAL, POINTER :: tilefrac (:,:) ! factor for latent heat REAL(r_2), POINTER :: wbtot (:) ! total soil water (mm) + REAL(r_2), POINTER :: evapfbl (:,:) ! REAL(r_2), POINTER :: gammzz (:,:) ! heat capacity for each soil layer REAL(r_2), POINTER :: wb (:,:) ! volumetric soil moisture (solid+liq) REAL(r_2), POINTER :: wbice (:,:) ! soil ice diff --git a/src/coupled/shared/cable_canopy_type_mod.F90 b/src/coupled/shared/cable_canopy_type_mod.F90 index 2cf050de5..6e519f5b4 100644 --- a/src/coupled/shared/cable_canopy_type_mod.F90 +++ b/src/coupled/shared/cable_canopy_type_mod.F90 @@ -62,7 +62,6 @@ MODULE cable_canopy_type_mod ga_cor ! correction to ground heat flux (W/m2) REAL, DIMENSION(:,:), POINTER :: & - evapfbl, & gswx, & ! stom cond for water zetar, & ! stability parameter (ref height) ! vh_js ! @@ -158,7 +157,6 @@ SUBROUTINE alloc_canopy_type(var, mp) ALLOCATE( var% fwet(mp) ) ALLOCATE( var% fns_cor(mp) ) !REV_CORR variable ALLOCATE( var% ga_cor(mp) ) !REV_CORR variable -ALLOCATE ( var % evapfbl(mp,ms) ) ALLOCATE( var% epot(mp) ) ALLOCATE( var% fnpp(mp) ) ALLOCATE( var% fevw_pot(mp) ) diff --git a/src/coupled/shared/cable_soilsnow_type_mod.F90 b/src/coupled/shared/cable_soilsnow_type_mod.F90 index 789c6a085..85dc668ef 100644 --- a/src/coupled/shared/cable_soilsnow_type_mod.F90 +++ b/src/coupled/shared/cable_soilsnow_type_mod.F90 @@ -72,7 +72,6 @@ MODULE cable_soil_snow_type_mod tggsn, & ! snow temperature in K dtmlt, & ! water flux to the soil albsoilsn, & ! soil + snow reflectance - evapfbl, & ! tilefrac ! factor for latent heat @@ -80,6 +79,7 @@ MODULE cable_soil_snow_type_mod wbtot ! total soil water (mm) REAL(r_2), DIMENSION(:,:), POINTER :: & + evapfbl, & ! gammzz, & ! heat capacity for each soil layer wb, & ! volumetric soil moisture (solid+liq) wbice, & ! soil ice From 5f352945bb8f013302cd7c9ad070b2c1b662f071 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:49:14 +1000 Subject: [PATCH 22/28] Done science/soilsnow --- src/science/soilsnow/cbl_GW.F90 | 314 +++++++++++++----- src/science/soilsnow/cbl_smoisturev.F90 | 8 +- src/science/soilsnow/cbl_snowAccum.F90 | 4 + src/science/soilsnow/cbl_snowCheck.F90 | 10 +- src/science/soilsnow/cbl_snowDensity.F90 | 8 +- src/science/soilsnow/cbl_snowMelt.F90 | 4 +- src/science/soilsnow/cbl_snowl_adjust.F90 | 6 +- src/science/soilsnow/cbl_soilfreeze.F90 | 6 +- .../soilsnow/cbl_soilsnow_init_special.F90 | 2 +- src/science/soilsnow/cbl_stempv.F90 | 4 +- src/science/soilsnow/cbl_surfbv.F90 | 11 +- src/science/soilsnow/cbl_thermal.F90 | 3 +- 12 files changed, 284 insertions(+), 96 deletions(-) diff --git a/src/science/soilsnow/cbl_GW.F90 b/src/science/soilsnow/cbl_GW.F90 index 6511b3017..8910d383b 100644 --- a/src/science/soilsnow/cbl_GW.F90 +++ b/src/science/soilsnow/cbl_GW.F90 @@ -7,6 +7,10 @@ MODULE GWstempv_mod CONTAINS SUBROUTINE GWstempv(dels, canopy, ssnow, soil) + +!*## Purpose +! updates soil temp and ground heat flux + USE cable_common_module, ONLY: cable_user USE total_soil_conductivity_mod, ONLY: total_soil_conductivity USE old_soil_conductivity_mod, ONLY: old_soil_conductivity @@ -19,18 +23,22 @@ SUBROUTINE GWstempv(dels, canopy, ssnow, soil) TYPE(soil_parameter_type), INTENT(INOUT) :: soil - REAL, DIMENSION(mp) :: & +! REAL, DIMENSION(mp) :: & ! rk4417 - phase2 + REAL(r_2), DIMENSION(mp) :: & coefa, coefb, & ! sgamm ! - REAL, DIMENSION(mp) :: & +! REAL, DIMENSION(mp) :: & ! rk4417 - phase2 + REAL(r_2), DIMENSION(mp) :: & dtg, & ! - ew, & ! - xx, & ! - wblfsp ! +! ew, & ! not used ! rk4417 - phase2 + xx !, & ! +! wblfsp ! not used ! rk4417 - phase2 REAL(r_2), DIMENSION(mp,ms) :: & - ccnsw ! soil thermal conductivity (incl water/ice) + ccnsw,& ! soil thermal conductivity (incl water/ice) + gammzz_snow + REAL(r_2), DIMENSION(mp, -2:ms) :: & at, bt, ct, rhs ! @@ -39,15 +47,28 @@ SUBROUTINE GWstempv(dels, canopy, ssnow, soil) REAL(r_2), DIMENSION(mp,ms+3) :: tmp_mat ! temp. matrix for tggsn & tgg - INTEGER :: j,k - REAL :: exp_arg + INTEGER :: j,k,i + REAL(r_2) :: exp_arg,dels_r2 LOGICAL :: direct2min = .FALSE. - at = 0.0 - bt = 1.0 - ct = 0.0 - coeff = 0.0 + dels_r2 = REAL(dels,r_2) ! rk4417 - phase2 + + at = 0._r_2 ! MMY@23Apr2023 are these taken from CABLE-GW? + bt = 1._r_2 ! accept them gammzz_snow is used the eq later + ct = 0._r_2 + coeff = 0._r_2 + + ssnow%otgg(:,:) = ssnow%tgg(:,:) ! MMY??? ssnow%otgg has gotten value in SUBROUTINE soil_snow_gw before call snow_processes_soil_thermal + + gammzz_snow(:,:) = 0._r_2 + + k=1 + DO i=1,mp + IF (ssnow%isflag(i) /= 0) THEN + gammzz_snow(i,k) = REAL(Ccgsnow * ssnow%snowd(i),r_2) + END IF + END DO IF (cable_user%soil_thermal_fix) THEN @@ -60,76 +81,162 @@ SUBROUTINE GWstempv(dels, canopy, ssnow, soil) xx(:) = 0. + +! WHERE(ssnow%isflag == 0) ! rk4417 - phase2 +! xx(:) = MAX( 0., ssnow%snowd / ssnow%ssdnn ) +! ccnsw(:,1) = ( ccnsw(:,1) - 0.2 ) * ( soil%zse(1) / ( soil%zse(1) + xx(:) ) & +! ) + 0.2 +! END WHERE + WHERE(ssnow%isflag == 0) - xx(:) = MAX( 0., ssnow%snowd / ssnow%ssdnn ) - ccnsw(:,1) = ( ccnsw(:,1) - 0.2 ) * ( soil%zse(1) / ( soil%zse(1) + xx(:) ) & - ) + 0.2 + xx(:) = MAX( 0._r_2, real(ssnow%snowd / ssnow%ssdnn,r_2) ) + ccnsw(:,1) = ( ccnsw(:,1) - 0.2_r_2 ) * ( soil%zse_vec(:,1) / ( soil%zse_vec(:,1) + xx(:) ) & + ) + 0.2_r_2 END WHERE - DO k = 3, ms + +! DO k = 3, ms ! rk4417 - phase2 +! WHERE (ssnow%isflag == 0) +! coeff(:,k) = 2.0 / ( soil%zse(k-1) / ccnsw(:,k-1) + soil%zse(k) / & +! ccnsw(:,k) ) +! END WHERE +! END DO + DO k = 3, ms WHERE (ssnow%isflag == 0) - coeff(:,k) = 2.0 / ( soil%zse(k-1) / ccnsw(:,k-1) + soil%zse(k) / & + coeff(:,k) = 2.0 / ( soil%zse_vec(:,k-1) / ccnsw(:,k-1) + soil%zse_vec(:,k) / & ccnsw(:,k) ) END WHERE END DO + + +! k = 1 ! rk4417 - phase2 +! WHERE( ssnow%isflag == 0 ) +! coeff(:,2) = 2.0 / ( ( soil%zse(1) + xx(:) ) / ccnsw(:,1) + soil%zse(2) / & +! ccnsw(:,2) ) +! coefa = 0.0 +! coefb = REAL( coeff(:,2) ) +! +! wblfsp = ssnow%wblf(:,k) +! +! ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)), & +! ( 1.0 - soil%ssat_vec(:,k) ) * & +! soil%css_vec(:,k) * soil%rhosoil_vec(:,k) & +! + soil%ssat_vec(:,k) * ( wblfsp * Ccs_rho_wat + & +! ssnow%wbfice(:,k) * Ccs_rho_ice ) ) & +! * soil%zse_vec(:,k) +! +! ssnow%gammzz(:,k) = ssnow%gammzz(:,k) + Ccgsnow * ssnow%snowd +! +! dtg = dels / ssnow%gammzz(:,k) +! +! at(:,k) = - dtg * coeff(:,k) +! ct(:,k) = - dtg * coeff(:,k+1) ! c3(ms)=0 & not really used +! bt(:,k) = 1.0 - at(:,k) - ct(:,k) +! END WHERE + k = 1 WHERE( ssnow%isflag == 0 ) - coeff(:,2) = 2.0 / ( ( soil%zse(1) + xx(:) ) / ccnsw(:,1) + soil%zse(2) / & - ccnsw(:,2) ) - coefa = 0.0 - coefb = REAL( coeff(:,2) ) - - wblfsp = ssnow%wblf(:,k) + coeff(:,2) = 2._r_2 / ( ( soil%zse_vec(:,1) + xx(:) ) / ccnsw(:,1) + soil%zse_vec(:,2) / & + ccnsw(:,2) ) + coefa = 0._r_2 + coefb = coeff(:,2) + ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)), & ( 1.0 - soil%ssat_vec(:,k) ) * & soil%css_vec(:,k) * soil%rhosoil_vec(:,k) & - + soil%ssat_vec(:,k) * ( wblfsp * Ccs_rho_wat + & - ssnow%wbfice(:,k) * Ccs_rho_ice ) ) & - * soil%zse_vec(:,k) - - ssnow%gammzz(:,k) = ssnow%gammzz(:,k) + Ccgsnow * ssnow%snowd - - dtg = dels / ssnow%gammzz(:,k) - + + ssnow%wbliq(:,k)*real(Ccswat*Cdensity_liq,r_2) & + !+ ssnow%wbice(:,k)*real(C%csice*C%density_liq*0.9,r_2) ) & ! MMY + + ssnow%wbice(:,k)*real(Ccsice*Cdensity_ice,r_2) ) & ! MMY + * soil%zse_vec(:,k) + gammzz_snow(:,k) + + dtg = dels_r2 / ssnow%gammzz(:,k) + at(:,k) = - dtg * coeff(:,k) ct(:,k) = - dtg * coeff(:,k+1) ! c3(ms)=0 & not really used bt(:,k) = 1.0 - at(:,k) - ct(:,k) - END WHERE - DO k = 2, ms - WHERE( ssnow%isflag == 0 ) - wblfsp = ssnow%wblf(:,k) +! DO k = 2, ms ! rk4417 - phase2 +! +! WHERE( ssnow%isflag == 0 ) +! +! wblfsp = ssnow%wblf(:,k) +! +! ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)), & +! ( 1.0 - soil%ssat_vec(:,k) ) * & +! soil%css_vec(:,k) * soil%rhosoil_vec(:,k) & +! + soil%ssat_vec(:,k) * ( wblfsp * Ccs_rho_wat + & +! ssnow%wbfice(:,k) * Ccs_rho_ice ) ) & +! * soil%zse_vec(:,k) +! +! dtg = dels / ssnow%gammzz(:,k) +! at(:,k) = - dtg * coeff(:,k) +! ct(:,k) = - dtg * coeff(:,k+1) ! c3(ms)=0 & not really used +! bt(:,k) = 1.0 - at(:,k) - ct(:,k) +! +! END WHERE +! +! END DO + + DO k = 2, ms + WHERE( ssnow%isflag == 0 ) + ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)), & ( 1.0 - soil%ssat_vec(:,k) ) * & soil%css_vec(:,k) * soil%rhosoil_vec(:,k) & - + soil%ssat_vec(:,k) * ( wblfsp * Ccs_rho_wat + & - ssnow%wbfice(:,k) * Ccs_rho_ice ) ) & - * soil%zse_vec(:,k) - - dtg = dels / ssnow%gammzz(:,k) + + ssnow%wbliq(:,k)*real(Ccswat*Cdensity_liq,r_2) & + !+ ssnow%wbice(:,k)*real(C%csice*C%density_liq*0.9,r_2) ) & ! MMY + + ssnow%wbice(:,k)*real(Ccsice*Cdensity_ice,r_2) ) & ! MMY + * soil%zse_vec(:,k) + gammzz_snow(:,k) + + dtg = dels_r2 / ssnow%gammzz(:,k) + at(:,k) = - dtg * coeff(:,k) ct(:,k) = - dtg * coeff(:,k+1) ! c3(ms)=0 & not really used bt(:,k) = 1.0 - at(:,k) - ct(:,k) - + END WHERE END DO + +! WHERE( ssnow%isflag == 0 ) ! rk4417 - phase2 +! bt(:,1) = bt(:,1) - canopy%dgdtg * dels / ssnow%gammzz(:,1) +! ssnow%tgg(:,1) = ssnow%tgg(:,1) + ( canopy%ga - ssnow%tgg(:,1) & +! * REAL( canopy%dgdtg ) ) * dels / REAL( ssnow%gammzz(:,1) ) +! END WHERE + WHERE( ssnow%isflag == 0 ) - bt(:,1) = bt(:,1) - canopy%dgdtg * dels / ssnow%gammzz(:,1) - ssnow%tgg(:,1) = ssnow%tgg(:,1) + ( canopy%ga - ssnow%tgg(:,1) & - * REAL( canopy%dgdtg ) ) * dels / REAL( ssnow%gammzz(:,1) ) + bt(:,1) = bt(:,1) - canopy%dgdtg * dels_r2 / ssnow%gammzz(:,1) + ssnow%tgg(:,1) = ssnow%tgg(:,1) + real(( real(canopy%ga,r_2) - real(ssnow%tgg(:,1),r_2) & + * REAL( canopy%dgdtg ) ) * dels_r2 / ssnow%gammzz(:,1) ) END WHERE coeff(:,1-3) = 0.0 ! coeff(:,-2) +! ! 3-layer snow points done here ! rk4417 - phase2 +! WHERE( ssnow%isflag /= 0 ) +! +! ssnow%sconds(:,1) = MAX( 0.2, MIN( 2.876e-6 * ssnow%ssdn(:,1)**2 & +! + 0.074, max_sconds ) ) +! ssnow%sconds(:,2) = MAX(0.2, MIN(2.876e-6 * ssnow%ssdn(:,2)**2 & +! & + 0.074, max_sconds) ) +! ssnow%sconds(:,3) = MAX(0.2, MIN(2.876e-6 * ssnow%ssdn(:,3)**2 & +! & + 0.074, max_sconds) ) +! coeff(:,-1) = 2.0 / (ssnow%sdepth(:,1) / ssnow%sconds(:,1) & +! & + ssnow%sdepth(:,2) / ssnow%sconds(:,2) ) +! coeff(:,0) = 2.0 / (ssnow%sdepth(:,2) / ssnow%sconds(:,2) & +! & + ssnow%sdepth(:,3) / ssnow%sconds(:,3) ) +! coeff(:,1) = 2.0 / (ssnow%sdepth(:,3) / ssnow%sconds(:,3) & +! & + soil%zse(1) / ccnsw (:,1) ) +! END WHERE + ! 3-layer snow points done here WHERE( ssnow%isflag /= 0 ) @@ -139,32 +246,55 @@ SUBROUTINE GWstempv(dels, canopy, ssnow, soil) & + 0.074, max_sconds) ) ssnow%sconds(:,3) = MAX(0.2, MIN(2.876e-6 * ssnow%ssdn(:,3)**2 & & + 0.074, max_sconds) ) - coeff(:,-1) = 2.0 / (ssnow%sdepth(:,1) / ssnow%sconds(:,1) & - & + ssnow%sdepth(:,2) / ssnow%sconds(:,2) ) - coeff(:,0) = 2.0 / (ssnow%sdepth(:,2) / ssnow%sconds(:,2) & - & + ssnow%sdepth(:,3) / ssnow%sconds(:,3) ) - coeff(:,1) = 2.0 / (ssnow%sdepth(:,3) / ssnow%sconds(:,3) & - & + soil%zse(1) / ccnsw (:,1) ) + coeff(:,-1) = 2._r_2 / (real(ssnow%sdepth(:,1) / ssnow%sconds(:,1),r_2) & + & + real(ssnow%sdepth(:,2) / ssnow%sconds(:,2),r_2) ) + coeff(:,0) = 2._r_2 / (real(ssnow%sdepth(:,2) / ssnow%sconds(:,2),r_2) & + & + real(ssnow%sdepth(:,3) / ssnow%sconds(:,3),r_2) ) + coeff(:,1) = 2._r_2 / (real(ssnow%sdepth(:,3) / ssnow%sconds(:,3),r_2) & + & + soil%zse_vec(:,1) / ccnsw (:,1) ) END WHERE + +! DO k = 2, ms ! rk4417 - phase2 +! WHERE( ssnow%isflag /= 0 ) & +! coeff(:,k) = 2.0 / ( soil%zse(k-1) / ccnsw(:,k-1) + soil%zse(k) / & +! ccnsw(:,k) ) +! END DO + DO k = 2, ms WHERE( ssnow%isflag /= 0 ) & - coeff(:,k) = 2.0 / ( soil%zse(k-1) / ccnsw(:,k-1) + soil%zse(k) / & + coeff(:,k) = 2._r_2 / ( soil%zse_vec(:,k-1) / ccnsw(:,k-1) + soil%zse_vec(:,k) / & ccnsw(:,k) ) END DO +! WHERE( ssnow%isflag /= 0 ) ! rk4417 - phase2 +! coefa = REAL( coeff (:,-1) ) +! coefb = REAL( coeff (:,1) ) +! END WHERE + WHERE( ssnow%isflag /= 0 ) - coefa = REAL( coeff (:,-1) ) - coefb = REAL( coeff (:,1) ) + coefa = coeff (:,-1) + coefb = coeff (:,1) END WHERE + +! DO k = 1, 3 ! rk4417 - phase2 +! WHERE( ssnow%isflag /= 0 ) +! sgamm = ssnow%ssdn(:,k) * Ccgsnow * ssnow%sdepth(:,k) +! dtg = dels / sgamm +! at(:,k-3) = - dtg * coeff(:,k-3) +! ct(:,k-3) = - dtg * coeff(:,k-2) +! bt(:,k-3) = 1.0 - at(:,k-3) - ct(:,k-3) +! END WHERE +! END DO + DO k = 1, 3 WHERE( ssnow%isflag /= 0 ) - sgamm = ssnow%ssdn(:,k) * Ccgsnow * ssnow%sdepth(:,k) - dtg = dels / sgamm + sgamm = real(ssnow%ssdn(:,k) * Ccgsnow * ssnow%sdepth(:,k),r_2) + dtg = dels_r2 / sgamm at(:,k-3) = - dtg * coeff(:,k-3) ct(:,k-3) = - dtg * coeff(:,k-2) bt(:,k-3) = 1.0 - at(:,k-3) - ct(:,k-3) @@ -172,36 +302,64 @@ SUBROUTINE GWstempv(dels, canopy, ssnow, soil) END DO - DO k = 1, ms + +! DO k = 1, ms ! rk4417 - phase2 +! WHERE( ssnow%isflag /= 0 ) +! wblfsp = ssnow%wblf(:,k) +! +! ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)),& +! ( 1.0 - soil%ssat_vec(:,k) ) * soil%css_vec(:,k) * & +! soil%rhosoil_vec(:,k) + soil%ssat_vec(:,k) * ( wblfsp * Ccs_rho_wat +& +! ssnow%wbfice(:,k) * Ccs_rho_ice)) * & +! soil%zse_vec(:,k) +! +! dtg = dels / ssnow%gammzz(:,k) +! at(:,k) = - dtg * coeff(:,k) +! ct(:,k) = - dtg * coeff(:,k + 1) ! c3(ms)=0 & not really used +! bt(:,k) = 1.0 - at(:,k) - ct(:,k) +! END WHERE +! END DO + DO k = 1, ms WHERE( ssnow%isflag /= 0 ) - wblfsp = ssnow%wblf(:,k) - - ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)),& - ( 1.0 - soil%ssat_vec(:,k) ) * soil%css_vec(:,k) * & - soil%rhosoil_vec(:,k) + soil%ssat_vec(:,k) * ( wblfsp * Ccs_rho_wat +& - ssnow%wbfice(:,k) * Ccs_rho_ice)) * & - soil%zse_vec(:,k) - dtg = dels / ssnow%gammzz(:,k) + ssnow%gammzz(:,k) = MAX((soil%heat_cap_lower_limit(:,k)), & + ( 1.0 - soil%ssat_vec(:,k) ) * & + soil%css_vec(:,k) * soil%rhosoil_vec(:,k) & + + ssnow%wbliq(:,k)*real(Ccswat*Cdensity_liq,r_2) & + !+ ssnow%wbice(:,k)*real(C%csice*C%density_liq*0.9,r_2) ) & ! MMY + + ssnow%wbice(:,k)*real(Ccsice*Cdensity_ice,r_2) ) & ! MMY + * soil%zse_vec(:,k) + gammzz_snow(:,k) + + dtg = dels_r2 / ssnow%gammzz(:,k) at(:,k) = - dtg * coeff(:,k) ct(:,k) = - dtg * coeff(:,k + 1) ! c3(ms)=0 & not really used - bt(:,k) = 1.0 - at(:,k) - ct(:,k) - + bt(:,k) = 1._r_2 - at(:,k) - ct(:,k) END WHERE - END DO - WHERE( ssnow%isflag /= 0 ) - sgamm = ssnow%ssdn(:,1) * Ccgsnow * ssnow%sdepth(:,1) - - bt(:,-2) = bt(:,-2) - canopy%dgdtg * dels / sgamm - - ssnow%tggsn(:,1) = ssnow%tggsn(:,1) + ( canopy%ga - ssnow%tggsn(:,1 ) & - * REAL( canopy%dgdtg ) ) * dels / sgamm + +! WHERE( ssnow%isflag /= 0 ) ! rk4417 - phase2 +! sgamm = ssnow%ssdn(:,1) * Ccgsnow * ssnow%sdepth(:,1) +! +! bt(:,-2) = bt(:,-2) - canopy%dgdtg * dels / sgamm +! +! ssnow%tggsn(:,1) = ssnow%tggsn(:,1) + ( canopy%ga - ssnow%tggsn(:,1 ) & +! * REAL( canopy%dgdtg ) ) * dels / sgamm +! +! rhs(:,1-3) = ssnow%tggsn(:,1) +! END WHERE + WHERE( ssnow%isflag /= 0 ) + sgamm = real(ssnow%ssdn(:,1) * Ccgsnow * ssnow%sdepth(:,1),r_2) + + bt(:,-2) = bt(:,-2) - canopy%dgdtg * dels_r2 / sgamm + + ssnow%tggsn(:,1) = ssnow%tggsn(:,1) +real( ( real(canopy%ga,r_2) - real(ssnow%tggsn(:,1),r_2) & + * (canopy%dgdtg) * dels_r2) / sgamm ) + rhs(:,1-3) = ssnow%tggsn(:,1) - END WHERE + END WHERE ! note in the following that tgg and tggsn are processed together @@ -212,8 +370,10 @@ SUBROUTINE GWstempv(dels, canopy, ssnow, soil) ssnow%tggsn = REAL( tmp_mat(:,1:3) ) ssnow%tgg = REAL( tmp_mat(:,4:(ms+3)) ) - canopy%sghflux = coefa * ( ssnow%tggsn(:,1) - ssnow%tggsn(:,2) ) - canopy%ghflux = coefb * ( ssnow%tgg(:,1) - ssnow%tgg(:,2) ) ! +ve downwards +! canopy%sghflux = coefa * ( ssnow%tggsn(:,1) - ssnow%tggsn(:,2) ) ! rk4417 - phase2 +! canopy%ghflux = coefb * ( ssnow%tgg(:,1) - ssnow%tgg(:,2) ) ! +ve downwards + canopy%sghflux = real(coefa) * ( ssnow%tggsn(:,1) - ssnow%tggsn(:,2) ) + canopy%ghflux = real(coefb) * ( ssnow%tgg(:,1) - ssnow%tgg(:,2) ) ! +ve downwards END SUBROUTINE GWstempv diff --git a/src/science/soilsnow/cbl_smoisturev.F90 b/src/science/soilsnow/cbl_smoisturev.F90 index 635e4f06b..93e4c0fe9 100644 --- a/src/science/soilsnow/cbl_smoisturev.F90 +++ b/src/science/soilsnow/cbl_smoisturev.F90 @@ -309,9 +309,13 @@ SUBROUTINE smoisturev (dels,ssnow,soil,veg) ct(:,k) = dtt(:,k) * ( - z2(:,k+1) * 0.5 * soil%zse(k) & / soil%zshh (k+1) - z3(:,k+1) ) +! bt(:,k) = 1.0 + dtt(:,k) * ( - z2(:,k+1) * 0.5 * soil%zse(k+1) & ! replaced by rk4417 - phase2 +! / soil%zshh (k+1) + z2(:,k) * 0.5 * soil%zse(k) & +! / soil%zshh (k) + z3(:,k+1) + z3(:,k) ) + bt(:,k) = 1.0 + dtt(:,k) * ( - z2(:,k+1) * 0.5 * soil%zse(k+1) & - / soil%zshh (k+1) + z2(:,k) * 0.5 * soil%zse(k) & - / soil%zshh (k) + z3(:,k+1) + z3(:,k) ) + / soil%zshh (k+1) + z2(:,k) * 0.5 * soil%zse( MAX( k-1, & + 1 ) ) / soil%zshh (k) + z3(:,k+1) + z3(:,k) ) END DO diff --git a/src/science/soilsnow/cbl_snowAccum.F90 b/src/science/soilsnow/cbl_snowAccum.F90 index 27ed28d03..adcf8285a 100644 --- a/src/science/soilsnow/cbl_snowAccum.F90 +++ b/src/science/soilsnow/cbl_snowAccum.F90 @@ -8,6 +8,10 @@ MODULE snow_accum_mod SUBROUTINE snow_accum ( dels, canopy, met, ssnow, soil ) + !*## Purpose + ! calculate snowfall and snow evaporation and update snow depth, snow temp, + ! snow mass, snow density + USE cable_common_module IMPLICIT NONE diff --git a/src/science/soilsnow/cbl_snowCheck.F90 b/src/science/soilsnow/cbl_snowCheck.F90 index 2b4d64800..6059453d8 100644 --- a/src/science/soilsnow/cbl_snowCheck.F90 +++ b/src/science/soilsnow/cbl_snowCheck.F90 @@ -7,10 +7,12 @@ MODULE snowcheck_mod CONTAINS SUBROUTINE snowcheck(dels, ssnow, soil, met ) - + !*## Purpose + ! Set up snow depth, snow mass, snow temp and snow layer used + USE cable_common_module -IMPLICIT NONE + IMPLICIT NONE REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow @@ -23,6 +25,7 @@ SUBROUTINE snowcheck(dels, ssnow, soil, met ) DO j=1,mp IF( ssnow%snowd(j) <= 0.0 ) THEN + ! using a single snow layer but there is no snow yet ssnow%isflag(j) = 0 ssnow%ssdn(j,:) = 120.0 @@ -39,6 +42,7 @@ SUBROUTINE snowcheck(dels, ssnow, soil, met ) ! in loop: IF( ssnow%snowd(j) <= 0.0 ) THEN ELSEIF( ssnow%snowd(j) < snmin * ssnow%ssdnn(j) ) THEN + ! snow depth is between 0 and 1*snow density IF( ssnow%isflag(j) == 1 ) THEN ssnow%ssdn(j,1) = ssnow%ssdnn(j) @@ -61,7 +65,6 @@ SUBROUTINE snowcheck(dels, ssnow, soil, met ) ssnow%ssdn(j,:) = ssnow%ssdnn(j) - ELSE ! in loop: IF( ssnow%snowd(j) <= 0.0 ) THEN ! sufficient snow now for 3 layer snowpack @@ -72,7 +75,6 @@ SUBROUTINE snowcheck(dels, ssnow, soil, met ) ssnow%ssdn(j,2) = ssnow%ssdn(j,1) ssnow%ssdn(j,3) = ssnow%ssdn(j,1) - ssnow%sdepth(j,1) = ssnow%t_snwlr(j) ssnow%smass(j,1) = ssnow%t_snwlr(j) * ssnow%ssdn(j,1) diff --git a/src/science/soilsnow/cbl_snowDensity.F90 b/src/science/soilsnow/cbl_snowDensity.F90 index e7674c68c..854bc6045 100644 --- a/src/science/soilsnow/cbl_snowDensity.F90 +++ b/src/science/soilsnow/cbl_snowDensity.F90 @@ -7,8 +7,10 @@ MODULE snowdensity_mod CONTAINS SUBROUTINE snowdensity (dels, ssnow, soil) - -IMPLICIT NONE + !*## Purpose + ! Calculate snow density for either single snow layer and three snow layers + + IMPLICIT NONE REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow @@ -16,7 +18,7 @@ SUBROUTINE snowdensity (dels, ssnow, soil) TYPE(soil_parameter_type), INTENT(INOUT) :: soil REAL, DIMENSION(mp) :: ssnow_tgg_min1 - REAL, DIMENSION(mp,3) :: ssnow_tgg_min +! REAL, DIMENSION(mp,3) :: ssnow_tgg_min ! not used, commented out by rk4417 - phase2 ssnow_tgg_min1 = MIN( CTFRZ, ssnow%tgg(:,1) ) diff --git a/src/science/soilsnow/cbl_snowMelt.F90 b/src/science/soilsnow/cbl_snowMelt.F90 index 12ac79caa..3911c9ae8 100644 --- a/src/science/soilsnow/cbl_snowMelt.F90 +++ b/src/science/soilsnow/cbl_snowMelt.F90 @@ -7,9 +7,11 @@ MODULE snow_melting_mod CONTAINS SUBROUTINE snow_melting (dels, snowmlt, ssnow, soil ) + !*## Purpose + ! Snow melting USE cable_common_module -IMPLICIT NONE + IMPLICIT NONE REAL, INTENT(IN) :: dels ! integration time step (s) diff --git a/src/science/soilsnow/cbl_snowl_adjust.F90 b/src/science/soilsnow/cbl_snowl_adjust.F90 index 023f2f50c..9a07b0c11 100644 --- a/src/science/soilsnow/cbl_snowl_adjust.F90 +++ b/src/science/soilsnow/cbl_snowl_adjust.F90 @@ -7,8 +7,12 @@ MODULE snowl_adjust_mod CONTAINS SUBROUTINE snowl_adjust(dels, ssnow, canopy ) + !*## Purpose + ! Adjust levels in the snowpack due to snow accumulation/melting, + ! snow aging etc... -IMPLICIT NONE + IMPLICIT NONE + REAL, INTENT(IN) :: dels ! integration time step (s) TYPE(soil_snow_type), INTENT(INOUT) :: ssnow diff --git a/src/science/soilsnow/cbl_soilfreeze.F90 b/src/science/soilsnow/cbl_soilfreeze.F90 index a1b1a3cde..567f4dfcf 100644 --- a/src/science/soilsnow/cbl_soilfreeze.F90 +++ b/src/science/soilsnow/cbl_soilfreeze.F90 @@ -15,9 +15,9 @@ SUBROUTINE soilfreeze(dels, soil, ssnow,heat_cap_lower_limit) REAL(r_2), DIMENSION(mp) :: sicefreeze REAL(r_2), DIMENSION(mp) :: sicemelt REAL, DIMENSION(mp) :: xx -INTEGER :: i,k -REAL :: heat_cap_lower_limit(mp,ms) -REAL :: max_arg1, max_arg2 + INTEGER :: i,k + REAL :: heat_cap_lower_limit(mp,ms) ! best to declare INTENT - rk4417 - phase2 + REAL :: max_arg1, max_arg2 xx = 0. DO k = 1, ms !loop over soil levels diff --git a/src/science/soilsnow/cbl_soilsnow_init_special.F90 b/src/science/soilsnow/cbl_soilsnow_init_special.F90 index 69e25503e..f5a4a93db 100644 --- a/src/science/soilsnow/cbl_soilsnow_init_special.F90 +++ b/src/science/soilsnow/cbl_soilsnow_init_special.F90 @@ -27,7 +27,7 @@ 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 -REAL :: heat_cap_lower_limit(mp,ms) +REAL :: heat_cap_lower_limit(mp,ms) ! best to declare INTENT - rk4417 - phase2 ! since CMIP5 only ever (potentially) TRUE offline if special initialization ! requested and then only on first timestep diff --git a/src/science/soilsnow/cbl_stempv.F90 b/src/science/soilsnow/cbl_stempv.F90 index cd5a280b2..ad081012b 100644 --- a/src/science/soilsnow/cbl_stempv.F90 +++ b/src/science/soilsnow/cbl_stempv.F90 @@ -47,13 +47,15 @@ SUBROUTINE stempv(dels, canopy, ssnow, soil,heat_cap_lower_limit) INTEGER :: j,k REAL :: exp_arg LOGICAL :: direct2min = .FALSE. -REAL :: heat_cap_lower_limit(mp,ms) + REAL :: heat_cap_lower_limit(mp,ms) ! best to declare INTENT here - rk4417 - phase2 at = 0.0 bt = 1.0 ct = 0.0 coeff = 0.0 + ssnow%otgg(:,:) = ssnow%tgg ! FEEDBACK (OK to insert this line as per MMY code?) --rk4417 + IF (cable_user%soil_thermal_fix) THEN ccnsw = total_soil_conductivity(ssnow,soil) ELSE diff --git a/src/science/soilsnow/cbl_surfbv.F90 b/src/science/soilsnow/cbl_surfbv.F90 index b6117bb70..c5c0f2bd7 100644 --- a/src/science/soilsnow/cbl_surfbv.F90 +++ b/src/science/soilsnow/cbl_surfbv.F90 @@ -27,7 +27,8 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) TYPE(soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters !jhan:cable.nml - INTEGER :: nglacier = 0 ! 0 original, 1 off, 2 new Eva +! INTEGER :: nglacier = 0 ! 0 original, 1 off, 2 new Eva + INTEGER :: nglacier ! 0 original, 1 off, 2 new Eva ! rk4417 - phase2 REAL, DIMENSION(mp) :: & rnof5, & ! @@ -42,7 +43,13 @@ SUBROUTINE surfbv (dels, met, ssnow, soil, veg, canopy ) REAL :: wb_lake_T, rnof2_T, ratio INTEGER :: k,j - IF( .NOT. cable_runtime%UM .OR. cable_runtime%esm15 ) THEN +! IF( .NOT. cable_runtime%UM .OR. cable_runtime%esm15 ) THEN ! rk4417 - phase2 +! nglacier = 2 +! ENDIF + + IF( cable_runtime%UM ) THEN + nglacier = 0 + ELSE nglacier = 2 ENDIF diff --git a/src/science/soilsnow/cbl_thermal.F90 b/src/science/soilsnow/cbl_thermal.F90 index cb682df1f..74fd4ec16 100644 --- a/src/science/soilsnow/cbl_thermal.F90 +++ b/src/science/soilsnow/cbl_thermal.F90 @@ -48,7 +48,8 @@ SUBROUTINE snow_processes_soil_thermal(dels,ssnow,soil,veg,canopy,met,bal) ! snow aging etc... CALL snowl_adjust(dels, ssnow, canopy ) - IF (cable_user%gw_model) CALL GWstempv(dels, canopy, ssnow, soil) +! IF (cable_user%gw_model) CALL GWstempv(dels, canopy, ssnow, soil) ! replaced by rk4417 - phase2 + CALL GWstempv(dels, canopy, ssnow, soil) !do the soil and snow melting, freezing prior to water movement DO i=1,mp From d7568331737c0a5f5ea2d7dae2a41f9cb5ba09cb Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:49:17 +1000 Subject: [PATCH 23/28] Done with science/canopy --- src/science/canopy/cable_canopy.F90 | 259 ++++++++++++---------- src/science/canopy/cbl_SurfaceWetness.F90 | 38 ++-- src/science/canopy/cbl_dryLeaf.F90 | 9 + src/science/canopy/cbl_fwsoil.F90 | 4 +- src/science/canopy/cbl_latent_heat.F90 | 42 +++- src/science/canopy/cbl_pot_evap_snow.F90 | 86 ++++--- src/science/canopy/cbl_wetleaf.F90 | 2 +- 7 files changed, 272 insertions(+), 168 deletions(-) diff --git a/src/science/canopy/cable_canopy.F90 b/src/science/canopy/cable_canopy.F90 index 96073e6dc..c7f1f4cb5 100644 --- a/src/science/canopy/cable_canopy.F90 +++ b/src/science/canopy/cable_canopy.F90 @@ -19,6 +19,10 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima USE cbl_dryLeaf_module, ONLY : dryLeaf USE cable_within_canopy_module, ONLY : within_canopy USE cbl_SurfaceWetness_module, ONLY : Surf_wetness_fact +USE cable_psm, ONLY: or_soil_evap_resistance,rtevap_max,& ! inserted by rk4417 - phase2 + rt_Dff,update_or_soil_resis +USE cable_gw_hydro_module, ONLY : pore_space_relative_humidity, den_rat ! Use public variable den_rat to avoid calling subroutine set_den_rat repeatedly ! line inserted by rk4417 - phase2 + ! data USE cable_common_module, ONLY: cable_runtime, cable_user @@ -65,11 +69,11 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima USE grid_constants_mod_cbl, ONLY: ICE_SoilType - USE cable_def_types_mod - USE cable_air_module - USE cable_roughness_module +USE cable_def_types_mod +USE cable_air_module +USE cable_roughness_module - IMPLICIT NONE +IMPLICIT NONE TYPE (balances_type), INTENT(INOUT) :: bal TYPE (radiation_type), INTENT(INOUT) :: rad @@ -83,8 +87,8 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima TYPE (soil_parameter_type), INTENT(INOUT) :: soil TYPE (veg_parameter_type), INTENT(INOUT) :: veg -REAL :: reducedLAIdue2snow(mp) -LOGICAL :: sunlit_veg_mask(mp) + REAL :: reducedLAIdue2snow(mp) + LOGICAL :: sunlit_veg_mask(mp) REAL, INTENT(IN) :: dels ! integration time setp (s) INTEGER :: & iter, & ! iteration # @@ -220,6 +224,13 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ssnow%tss = REAL((1-ssnow%isflag))*ssnow%tgg(:,1) + & REAL(ssnow%isflag)*ssnow%tggsn(:,1) ENDIF + + IF (cable_user%gw_model) THEN ! IF block inserted by rk4417 - phase2 + ssnow%wbliq(:,:) = ssnow%wb(:,:) - den_rat*ssnow%wbice(:,:) + ELSE + ssnow%wbliq(:,:) = ssnow%wb(:,:) - ssnow%wbice(:,:) + ENDIF + tss4 = ssnow%tss**4 canopy%fes = 0. canopy%fess = 0. @@ -364,8 +375,10 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima IF (cable_user%or_evap) THEN - write(6,*) "GW or ORevepis not an option right now" - !H! call or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) +! write(6,*) "GW or ORevepis not an option right now" - commented out by rk4417 - phase2 +!H! call or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) +! line above replaced by below - rk4417 - phase2 + CALL or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) END IF ! Vegetation boundary-layer conductance (mol/m2/s) @@ -430,7 +443,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima rad%lwabv(j) = CCAPP * Crmair * ( tlfy(j) - met%tk(j) ) * & sum_rad_gradis(j) ! vh_js ! - + IF ( (rad%lwabv(j) / (2.0*(1.0-rad%transd(j)) & * CSBOLTZ*CEMLEAF)+met%tvrad(j)**4) .GT. 0.0) THEN @@ -446,7 +459,6 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima canopy%tv(j) = met%tvrad(j) ENDIF - ENDDO @@ -460,9 +472,11 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ! Saturation specific humidity at soil/snow surface temperature: CALL qsatfjh(mp, ssnow%qstss, CRMH2o, Crmair, CTETENA, CTETENB, CTETENC,ssnow%tss-CTfrz,met%pmb) - if (cable_user%gw_model .OR. cable_user%or_evap) & - write(6,*) "GW or ORevepis not an option right now" + IF (cable_user%gw_model .OR. cable_user%or_evap) & +! write(6,*) "GW or ORevepis not an option right now" - commented out by rk4417 - phase2 !H! call pore_space_relative_humidity(ssnow,soil,veg) +! line above replaced by below - rk4417 - phase2 + CALL pore_space_relative_humidity(ssnow,soil,veg) IF (cable_user%soil_struc=='default') THEN @@ -500,7 +514,8 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima REAL(canopy%DvLitt), & ssnow%isflag, REAL(ssnow%satfrac),ssnow%rtsoil, & REAL(ssnow%rtevap_sat), REAL(ssnow%rtevap_unsat), & - ssnow%snowd, ssnow%tgg(:,1) ) + ssnow%snowd, ssnow%tgg(:,1), & + veg%iveg, rtevap_max, canopy%sublayer_dz, rt_Dff ) ! inserted by rk4417 - phase2 ENDIF @@ -512,12 +527,19 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ssnow%snowd, ssnow%wb(:,1), ssnow%wbice(:,1), & ssnow%pudsto, ssnow%pudsmx, ssnow%potev, & ssnow%wetfac, ssnow%evapfbl(:,1), ssnow%cls, & - ssnow%tss, canopy%fes, canopy%fess, canopy%fesp ) + ssnow%tss, canopy%fes, canopy%fess, canopy%fesp, & + cable_user%gw_model, den_rat, soil%watr(:,1) ) ! line inserted by rk4417 - phase2 ! Calculate soil sensible heat: ! INH: I think this should be - met%tvair !canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tk) /ssnow%rtsoil - IF (cable_user%gw_model .OR. cable_user%or_evap) THEN + +! IF (cable_user%gw_model .OR. cable_user%or_evap) THEN +! MMY comment out since rt_qh_sublayer is only given value when or-on +! rather than depends on cable_user%gw_model +! line above replaced by below - rk4417 - phase2 + IF (cable_user%or_evap) THEN ! MMY + canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tk) / & (ssnow%rtsoil + ssnow%rt_qh_sublayer) !note if or_evap and litter are true then litter resistance is @@ -534,11 +556,10 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ELSE -write(6,*) "SLI is not an option right now" + write(6,*) "SLI is not an option right now" - commented out by rk4417 - phase2 ! SLI SEB to get canopy%fhs, canopy%fess, canopy%ga ! (Based on old Tsoil, new canopy%tv, new canopy%fns) !H!CALL sli_main(1,dels,veg,soil,ssnow,met,canopy,air,rad,1) - ENDIF @@ -573,7 +594,8 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima REAL(canopy%DvLitt), & ssnow%isflag, REAL(ssnow%satfrac),ssnow%rtsoil, & REAL(ssnow%rtevap_sat), REAL(ssnow%rtevap_unsat), & - ssnow%snowd, ssnow%tgg(:,1) ) + ssnow%snowd, ssnow%tgg(:,1), & + veg%iveg, rtevap_max, canopy%sublayer_dz, rt_Dff ) ! inserted by rk4417 - phase2 ENDIF @@ -584,12 +606,18 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ssnow%snowd, ssnow%wb(:,1), ssnow%wbice(:,1), & ssnow%pudsto, ssnow%pudsmx, ssnow%potev, & ssnow%wetfac, ssnow%evapfbl(:,1), ssnow%cls, & - ssnow%tss, canopy%fes, canopy%fess, canopy%fesp ) + ssnow%tss, canopy%fes, canopy%fess, canopy%fesp,& + cable_user%gw_model, den_rat, soil%watr(:,1) ) ! line inserted by rk4417 - phase2 ! Soil sensible heat: !canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) /ssnow%rtsoil - IF (cable_user%gw_model .OR. cable_user%or_evap) THEN + +! IF (cable_user%gw_model .OR. cable_user%or_evap) THEN ! MMY comment out since rt_qh_sublayer is only given value when or-on + ! rather than depends on cable_user%gw_model +! line above replaced by below - rk4417 - phase2 + IF (cable_user%or_evap) THEN ! MMY + canopy%fhs = air%rho*CCAPP*(ssnow%tss - met%tvair) / & (ssnow%rtsoil + REAL(ssnow%rt_qh_sublayer)) @@ -610,7 +638,8 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima canopy%ga = canopy%fns-canopy%fhs-canopy%fes ! *ssnow%cls ELSE -write(6,*) "SLI is not an option right now" + write(6,*) "SLI is not an option right now" - commented out by rk4417 - phase2 + ! SLI SEB to get canopy%fhs, canopy%fess, canopy%ga ! (Based on old Tsoil, new canopy%tv, new canopy%fns) !H!CALL sli_main(1,dels,veg,soil,ssnow,met,canopy,air,rad,1) @@ -810,7 +839,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima MIN(1., ( (r_sc(j)+rhlitt(j)*canopy%us(j)) / MAX( 1., & rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j) & + rt1usc(j) + rhlitt(j)*canopy%us(j) )) ) - Ctfrz - ELSEIF (cable_user%or_evap .OR. cable_user%gw_model) THEN + ELSE IF (cable_user%or_evap) THEN canopy%tscrn(j) = ssnow%tss(j) + (met%tk(j) - ssnow%tss(j)) * & MIN(1., ( (ssnow%rt_qh_sublayer(j)*canopy%us(j) + r_sc(j) ) / & MAX( 1., rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j) & @@ -820,7 +849,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima MIN(1., (r_sc(j) / MAX( 1., & rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j) & + rt1usc(j))) ) - Ctfrz - ENDIF + END IF ENDIF @@ -851,7 +880,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima rough%rt0us(j) + rough%rt1usa(j) + rough%rt1usb(j) & + rt1usc(j) + relitt(j)*canopy%us(j) )) ) - ELSEIF (cable_user%or_evap .OR. cable_user%gw_model) THEN + ELSEIF (cable_user%or_evap) THEN !using alpm1 as a dumy variable alpm1(j) = REAL(& ssnow%satfrac(j)/(REAL(ssnow%rtsoil(j),r_2)+& @@ -908,110 +937,110 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ! d(canopy%fhs)/d(ssnow%tgg) ! d(canopy%fes)/d(dq) !IF (cable_user%soil_struc=='default') THEN - ssnow%dfn_dtg = (-1.)*4.*CEMSOIL*CSBOLTZ*tss4/ssnow%tss + ssnow%dfn_dtg = (-1.)*4.*CEMSOIL*CSBOLTZ*tss4/ssnow%tss - !INH: REV_CORR revised sensitivity terms working variable - rttsoil = ssnow%rtsoil - IF (cable_user%L_REV_CORR) THEN - WHERE (canopy%vlaiw > CLAI_THRESH) + !INH: REV_CORR revised sensitivity terms working variable + rttsoil = ssnow%rtsoil + IF (cable_user%L_REV_CORR) THEN + WHERE (canopy%vlaiw > CLAI_THRESH) !if %vlaiw<=%LAI_THRESH then %rt1 already added to %rtsoil rttsoil = rttsoil + rough%rt1 - ENDWHERE - ENDIF + ENDWHERE + ENDIF -IF (cable_user%gw_model .or. cable_user%or_evap) THEN + IF (cable_user%or_evap) THEN - ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ real(ssnow%rt_qh_sublayer)) +! ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ real(ssnow%rt_qh_sublayer)) - ! INH simplifying code for legibility - !ssnow%dfe_ddq = real(ssnow%satfrac)*air%rho*air%rlam*ssnow%cls/ & - ! (ssnow%rtsoil+ real(ssnow%rtevap_sat)) + - ! (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf)*& - ! air%rho*air%rlam*ssnow%cls/ (ssnow%rtsoil+ - ! real(ssnow%rtevap_unsat) ) - ssnow%dfe_ddq = real(ssnow%satfrac)/(ssnow%rtsoil+ real(ssnow%rtevap_sat)) & - + (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf) & - / (ssnow%rtsoil+ real(ssnow%rtevap_unsat) ) - - !mrd561 fixes. Do same thing as INH but has been tested. - IF (cable_user%L_REV_CORR) THEN - alpm1 = real(ssnow%satfrac/(real(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) + & - (1.0-ssnow%satfrac) / (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) ) - beta2 = real(ssnow%satfrac/(real(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) + & +! ! INH simplifying code for legibility +! !ssnow%dfe_ddq = real(ssnow%satfrac)*air%rho*air%rlam*ssnow%cls/ & +! ! (ssnow%rtsoil+ real(ssnow%rtevap_sat)) + +! ! (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf)*& +! ! air%rho*air%rlam*ssnow%cls/ (ssnow%rtsoil+ +! ! real(ssnow%rtevap_unsat) ) +! ssnow%dfe_ddq = real(ssnow%satfrac)/(ssnow%rtsoil+ real(ssnow%rtevap_sat)) & +! + (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf) & +! / (ssnow%rtsoil+ real(ssnow%rtevap_unsat) ) + + !mrd561 fixes. Do same thing as INH but has been tested. + IF (cable_user%L_REV_CORR) THEN + alpm1 = REAL(ssnow%satfrac/(REAL(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) + & + (1.0-ssnow%satfrac) / (REAL(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) ) + beta2 = REAL(ssnow%satfrac/(REAL(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) + & (1.0-ssnow%satfrac) * ssnow%rh_srf & - / (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) ) - WHERE (canopy%vlaiw > CLAI_THRESH) - alpm1 = alpm1 + 1._r_2/real(rough%rt1,r_2) - beta_div_alpm = beta2 / alpm1 !might need limit here - rttsoil = ssnow%rtsoil + rough%rt1 - ELSEWHERE!if there is no canopy then qa should not change - beta_div_alpm=0.0 !do not divide by aplm1 prevent issues - rttsoil = ssnow%rtsoil - ENDWHERE - ssnow%dfh_dtg = air%rho*CCAPP/(rttsoil + & - real(ssnow%rt_qh_sublayer)) - ssnow%dfe_ddq = real(ssnow%satfrac*(1.0-real(beta_div_alpm,r_2)) / & - (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) + & - (1.0-ssnow%satfrac)* (ssnow%rh_srf - real(beta_div_alpm,r_2)) / & - (real(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) ) - - ELSE ! IF (cable_user%L_REV_CORR) THEN - ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ real(ssnow%rt_qh_sublayer)) - ssnow%dfe_ddq = real(ssnow%satfrac)/(ssnow%rtsoil+ real(ssnow%rtevap_sat)) & - + (1.0-real(ssnow%satfrac))*real(ssnow%rh_srf) & - / (ssnow%rtsoil+ real(ssnow%rtevap_unsat) ) - ENDIF ! IF (cable_user%L_REV_CORR) THEN + / (REAL(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) ) + WHERE (canopy%vlaiw > CLAI_THRESH) + alpm1 = alpm1 + 1._r_2/REAL(rough%rt1,r_2) + beta_div_alpm = beta2 / alpm1 !might need limit here + rttsoil = ssnow%rtsoil + rough%rt1 + ELSEWHERE!if there is no canopy then qa should not change + beta_div_alpm=0.0 !do not divide by aplm1 prevent issues + rttsoil = ssnow%rtsoil + ENDWHERE + ssnow%dfh_dtg = air%rho*CCAPP/(rttsoil + & + REAL(ssnow%rt_qh_sublayer)) + ssnow%dfe_ddq = REAL(ssnow%satfrac*(1.0-REAL(beta_div_alpm,r_2)) / & + (REAL(ssnow%rtsoil,r_2)+ ssnow%rtevap_sat) + & + (1.0-ssnow%satfrac)* (ssnow%rh_srf - REAL(beta_div_alpm,r_2)) / & + (REAL(ssnow%rtsoil,r_2)+ ssnow%rtevap_unsat ) ) + + ELSE ! IF (cable_user%L_REV_CORR) THEN + ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ REAL(ssnow%rt_qh_sublayer)) + ssnow%dfe_ddq = REAL(ssnow%satfrac)/(ssnow%rtsoil+ REAL(ssnow%rtevap_sat)) & + + (1.0-real(ssnow%satfrac))*REAL(ssnow%rh_srf) & + / (ssnow%rtsoil+ REAL(ssnow%rtevap_unsat) ) + ENDIF ! IF (cable_user%L_REV_CORR) THEN - !cls applies for both REV_CORR false and true - ssnow%dfe_ddq = ssnow%dfe_ddq*air%rho*air%rlam*ssnow%cls + !cls applies for both REV_CORR false and true + ssnow%dfe_ddq = ssnow%dfe_ddq*air%rho*air%rlam*ssnow%cls - !REV_CORR: factor %wetfac needed for potev>0. and gw_model &/or snow cover - !NB %wetfac=1. if or_evap - IF (cable_user%L_REV_CORR) THEN - WHERE (ssnow%potev >= 0.) - ssnow%dfe_ddq = ssnow%dfe_ddq*ssnow%wetfac - ENDWHERE - ENDIF - - -ELSEIF (cable_user%litter) THEN ! IF (cable_user%gw_model .or. cable_user%or_evap) THEN - !vh_js! INH simplifying code for legibility and REV_CORR - !ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ & - ! real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP)) - !ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ & - ! (ssnow%rtsoil+ real((1-ssnow%isflag))*veg%clitt*0.003/canopy%DvLitt) + !REV_CORR: factor %wetfac needed for potev>0. and gw_model &/or snow cover + !NB %wetfac=1. if or_evap + IF (cable_user%L_REV_CORR) THEN + WHERE (ssnow%potev >= 0.) + ssnow%dfe_ddq = ssnow%dfe_ddq*ssnow%wetfac + ENDWHERE + ENDIF + + + ELSEIF (cable_user%litter) THEN ! IF (cable_user%gw_model .or. cable_user%or_evap) THEN + !vh_js! INH simplifying code for legibility and REV_CORR + !ssnow%dfh_dtg = air%rho*CCAPP/(ssnow%rtsoil+ & + ! real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP)) + !ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ & + ! (ssnow%rtsoil+ real((1-ssnow%isflag))*veg%clitt*0.003/canopy%DvLitt) - !recalculated - probably not needed - rhlitt = real((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP) - relitt = real((1-ssnow%isflag))*veg%clitt*0.003/canopy%DvLitt + !recalculated - probably not needed + rhlitt = REAL((1-ssnow%isflag))*veg%clitt*0.003/canopy%kthLitt/(air%rho*CCAPP) + relitt = REAL((1-ssnow%isflag))*veg%clitt*0.003/canopy%DvLitt - !incorporates REV_CORR changes - ssnow%dfh_dtg = air%rho*CCAPP/(rttsoil+rhlitt) - ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/(rttsoil+relitt) - - !REV_CORR: factor ssnow%wetfac is not applied if dew/frost i.e. potev<0 - IF (cable_user%L_REV_CORR) THEN - WHERE (ssnow%potev < 0.) - ssnow%dfe_ddq = air%rho*air%rlam*ssnow%cls/(rttsoil+relitt) - ENDWHERE - ENDIF - -ELSE ! i.e. NOT (%gw_model .or. %or_evap or SLI) - !ssnow%dfh_dtg = air%rho*CCAPP/ssnow%rtsoil - !ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ssnow%rtsoil + !incorporates REV_CORR changes + ssnow%dfh_dtg = air%rho*CCAPP/(rttsoil+rhlitt) + ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/(rttsoil+relitt) + + !REV_CORR: factor ssnow%wetfac is not applied if dew/frost i.e. potev<0 + IF (cable_user%L_REV_CORR) THEN + WHERE (ssnow%potev < 0.) + ssnow%dfe_ddq = air%rho*air%rlam*ssnow%cls/(rttsoil+relitt) + ENDWHERE + ENDIF + + ELSE ! i.e. NOT (%gw_model .or. %or_evap or SLI) + !ssnow%dfh_dtg = air%rho*CCAPP/ssnow%rtsoil + !ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/ssnow%rtsoil - !incorporates REV_CORR changes - ssnow%dfh_dtg = air%rho*CCAPP/rttsoil - ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/rttsoil + !incorporates REV_CORR changes + ssnow%dfh_dtg = air%rho*CCAPP/rttsoil + ssnow%dfe_ddq = ssnow%wetfac*air%rho*air%rlam*ssnow%cls/rttsoil - !REV_CORR: factor ssnow%wetfac is not applied if dew/frost i.e. potev<0 - IF (cable_user%L_REV_CORR) THEN - WHERE (ssnow%potev < 0.) - ssnow%dfe_ddq = air%rho*air%rlam*ssnow%cls/(rttsoil+relitt) - ENDWHERE - ENDIF - -ENDIF ! IF (cable_user%gw_model .or. cable_user%or_evap) THEN + !REV_CORR: factor ssnow%wetfac is not applied if dew/frost i.e. potev<0 + IF (cable_user%L_REV_CORR) THEN + WHERE (ssnow%potev < 0.) + ssnow%dfe_ddq = air%rho*air%rlam*ssnow%cls/(rttsoil+relitt) + ENDWHERE + ENDIF + + ENDIF ! IF (cable_user%gw_model .or. cable_user%or_evap) THEN ssnow%ddq_dtg = (Crmh2o/Crmair) /met%pmb * CTETENA*CTETENB * CTETENC & / ( ( CTETENC + ssnow%tss-Ctfrz )**2 )*EXP( CTETENB * & diff --git a/src/science/canopy/cbl_SurfaceWetness.F90 b/src/science/canopy/cbl_SurfaceWetness.F90 index b35530cf6..04581e689 100644 --- a/src/science/canopy/cbl_SurfaceWetness.F90 +++ b/src/science/canopy/cbl_SurfaceWetness.F90 @@ -15,6 +15,8 @@ SUBROUTINE Surf_wetness_fact( cansat, canopy, ssnow,veg, met, soil, dels ) USE cable_phys_constants_mod, ONLY : CTFRZ => TFRZ !H!USE cable_gw_hydro_module, ONLY : calc_srf_wet_fraction +! line above uncommented by rk4417 - phase2 +USE cable_gw_hydro_module, ONLY : calc_srf_wet_fraction USE cable_init_wetfac_mod, ONLY: initialize_wetfac TYPE (veg_parameter_type), INTENT(INOUT) :: veg @@ -56,24 +58,30 @@ SUBROUTINE Surf_wetness_fact( cansat, canopy, ssnow,veg, met, soil, dels ) !calc the surface wetness for soil evap in this routine !include the default wetfac when or_evap and gw_model are not used -!H!gw n/a here and so copied default below -!H! CALL calc_srf_wet_fraction(ssnow,soil,met,veg) -!H! ELSE !Default formulation - !call saturated_fraction(ssnow,soil,veg) - ssnow%satfrac(:) = 1.0e-8 - ssnow%rh_srf(:) = 1.0 + + CALL calc_srf_wet_fraction(ssnow,soil,met,veg) + +! line above replaces remaining code below - rk4417 - phase2 + +! !H!gw n/a here and so copied default below +! !H! CALL calc_srf_wet_fraction(ssnow,soil,met,veg) +! !H! ELSE !Default formulation +! + ! !call saturated_fraction(ssnow,soil,veg) + ! ssnow%satfrac(:) = 1.0e-8 + ! ssnow%rh_srf(:) = 1.0 - !This is updating wetfac iusing same calc as initialization - !originally code in canopy used 1e-6 as MIN - CALL initialize_wetfac( mp, ssnow%wetfac, soil%swilt, soil%sfc, & - ssnow%wb(:,1), ssnow%wbice(:,1), ssnow%snowd, & - veg%iveg, met%tk, Ctfrz ) + ! !This is updating wetfac iusing same calc as initialization + ! !originally code in canopy used 1e-6 as MIN + ! CALL initialize_wetfac( mp, ssnow%wetfac, soil%swilt, soil%sfc, & + ! ssnow%wb(:,1), ssnow%wbice(:,1), ssnow%snowd, & + ! veg%iveg, met%tk, Ctfrz ) - ! owetfac introduced to reduce sharp changes in dry regions, - ! especially in offline runs in which there may be discrepancies b/n - ! timing of precip and temperature change (EAK apr2009) - ssnow%wetfac = 0.5*(ssnow%wetfac + ssnow%owetfac) + ! ! owetfac introduced to reduce sharp changes in dry regions, + ! ! especially in offline runs in which there may be discrepancies b/n + ! ! timing of precip and temperature change (EAK apr2009) + ! ssnow%wetfac = 0.5*(ssnow%wetfac + ssnow%owetfac) RETURN END SUBROUTINE Surf_wetness_fact diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index c3a389459..9ca9d786c 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -149,6 +149,10 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & REAL, DIMENSION(mp,2) :: gsw_term, lower_limit2 ! local temp var +! Two lines inserted below - rk4417 - phase2 + REAL, DIMENSION(0:ms+1) :: diff ! MMY ! Martin's fix on water extraction from soil + REAL :: xx,xxd ! MMY + INTEGER :: i, j, k, kk ! iteration count ! For the calculation of the amount of transpired water @@ -246,6 +250,11 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & gwwet(i) = 1.075 * sum_gbh(i) ghrwet(i) = sum_rad_gradis(i) + ghwet(i) +! I checked with Claire...appears fine commented out - rk4417 - phase2 +! ! Calculate fraction of canopy which is wet: ! inserted by rk4417 - phase2 +! canopy%fwet(i) = MAX( 0.0, MIN( 1.0, 0.8 * canopy%cansto(i)/ MAX( & +! cansat(i),0.01 ) ) ) + ! Calculate lat heat from wet canopy, may be neg. ! if dew on wet canopy to avoid excessive evaporation: ccfevw(i) = MIN(canopy%cansto(i) * air%rlam(i) / dels, & diff --git a/src/science/canopy/cbl_fwsoil.F90 b/src/science/canopy/cbl_fwsoil.F90 index f83608bf3..d30fbf57b 100644 --- a/src/science/canopy/cbl_fwsoil.F90 +++ b/src/science/canopy/cbl_fwsoil.F90 @@ -49,8 +49,8 @@ SUBROUTINE fwsoil_calc_non_linear(fwsoil, soil, ssnow, veg) REAL, DIMENSION(mp,3) :: xi, ti, si INTEGER :: j - rwater = MAX(1.0e-9, & - SUM(veg%froot * MAX(0.0,MIN(1.0, REAL((ssnow%wbliq) - & + rwater = MAX(1.0e-9, & + SUM(veg%froot * MAX(1.0-9,MIN(1.0, REAL((ssnow%wbliq) - & SPREAD(soil%swilt, 2, ms)))),2) /(soil%sfc-soil%swilt)) fwsoil = 1. diff --git a/src/science/canopy/cbl_latent_heat.F90 b/src/science/canopy/cbl_latent_heat.F90 index 7329c2edc..0fdf85ff6 100644 --- a/src/science/canopy/cbl_latent_heat.F90 +++ b/src/science/canopy/cbl_latent_heat.F90 @@ -17,7 +17,8 @@ SUBROUTINE Latent_heat_flux( mp, CTFRZ, dels, soil_zse, soil_swilt, & ssnow_snowd, ssnow_wb, ssnow_wbice, & ssnow_pudsto, ssnow_pudsmx, ssnow_potev, & ssnow_wetfac, ssnow_evapfbl, ssnow_cls, & - ssnow_tss, canopy_fes, canopy_fess, canopy_fesp ) + ssnow_tss, canopy_fes, canopy_fess, canopy_fesp, & + cable_user_gw_model, den_rat, soil_watr ) ! line inserted by rk4417 - phase2 ) !*## Purpose ! @@ -106,6 +107,15 @@ SUBROUTINE Latent_heat_flux( mp, CTFRZ, dels, soil_zse, soil_swilt, & REAL, INTENT(IN) :: ssnow_tss(mp) !! temperature of surface soil/snow layer (K) +REAL(KIND=r_2), INTENT(IN) :: soil_watr(mp) ! line inserted by rk4417 - phase2 +!! residual water content of the soil (mm\(^{3}\)/mm\(^{3}\)) + +LOGICAL , INTENT(IN) :: cable_user_gw_model ! line inserted by rk4417 - phase2 +!! NAMELIST switch for gw model + +REAL(r_2) :: den_rat ! line inserted by rk4417 - phase2 + + REAL, DIMENSION(mp) :: & frescale, flower_limit, fupper_limit @@ -202,12 +212,29 @@ SUBROUTINE Latent_heat_flux( mp, CTFRZ, dels, soil_zse, soil_swilt, & ! switch. ! The options differ in the amount of water that remains at the end of the time step. ! - IF (.NOT.cable_user_l_new_reduce_soilevp) THEN - flower_limit(j) = REAL(ssnow_wb(j))-soil_swilt(j)/2.0 - ELSE + +! IF (.NOT.cable_user_l_new_reduce_soilevp) THEN +! flower_limit(j) = REAL(ssnow_wb(j))-soil_swilt(j)/2.0 +! ELSE +! ! E.Kowalczyk 2014 - reduces the soil evaporation +! flower_limit(j) = REAL(ssnow_wb(j))-soil_swilt(j) +! ENDIF + +! replaced IF block above by below - rk4417 - phase2 + + IF (.NOT.cable_user_l_new_reduce_soilevp) THEN + IF (cable_user_gw_model) THEN ! MMY + flower_limit(j) = REAL(ssnow_wb(j))-REAL(soil_watr(j)) ! MMY + ! MMY watr is better than swilt/2., as it has a clear physical meaning + ELSE ! MMY + flower_limit(j) = REAL(ssnow_wb(j))-soil_swilt(j)/2.0 + END IF ! MMY + ELSE ! E.Kowalczyk 2014 - reduces the soil evaporation flower_limit(j) = REAL(ssnow_wb(j))-soil_swilt(j) - ENDIF + END IF + + fupper_limit(j) = MAX( 0., & flower_limit(j) * frescale(j) & - ssnow_evapfbl(j)*air_rlam(j)/dels) @@ -221,8 +248,9 @@ SUBROUTINE Latent_heat_flux( mp, CTFRZ, dels, soil_zse, soil_swilt, & ! soil latent flux. **WARNING** frozen_limit=0.85 is hard coded - if it is changed ! then the corresponding limit in [[cbl_soilsnow]] must also be changed. - fupper_limit(j) = REAL(ssnow_wb(j)-ssnow_wbice(j)/frozen_limit)*frescale(j) - +! fupper_limit(j) = REAL(ssnow_wb(j)-ssnow_wbice(j)/0.85)*frescale(j) +! replaced line above by below - rk4417 - phase2 + fupper_limit(j) = REAL(ssnow_wb(j)-ssnow_wbice(j)*den_rat/0.85)*frescale(j) ! MMY keep fupper_limit consistent fupper_limit(j) = MAX(REAL(fupper_limit(j),r_2),0.) canopy_fess(j) = MIN(canopy_fess(j), REAL(fupper_limit(j),r_2)) diff --git a/src/science/canopy/cbl_pot_evap_snow.F90 b/src/science/canopy/cbl_pot_evap_snow.F90 index 0cc8cccd8..1ece5dea3 100644 --- a/src/science/canopy/cbl_pot_evap_snow.F90 +++ b/src/science/canopy/cbl_pot_evap_snow.F90 @@ -83,8 +83,13 @@ FUNCTION Humidity_deficit_method( mp, Ctfrz, veg_clitt,cable_user_or_evap, & canopy_DvLitt, & ssnow_isflag, ssnow_satfrac, ssnow_rtsoil, & ssnow_rtevap_sat, ssnow_rtevap_unsat, & - ssnow_snowd, ssnow_tgg & + ssnow_snowd, ssnow_tgg, & + veg_iveg, rtevap_max, canopy_sublayer_dz, & ! inserted by rk4417 - phase2 + rt_Dff & ! inserted by rk4417 - phase2 ) RESULT(ssnowpotev) + +USE cable_def_types_mod, ONLY : r_2 ! inserted by rk4417 - phase2 + IMPLICIT NONE INTEGER :: mp @@ -107,6 +112,11 @@ FUNCTION Humidity_deficit_method( mp, Ctfrz, veg_clitt,cable_user_or_evap, & REAL :: ssnow_rtevap_sat(mp) ! REAL :: ssnow_rtevap_unsat(mp) ! +INTEGER :: veg_iveg(mp) ! inserted by rk4417 - phase2 +REAL(r_2) :: rtevap_max ! inserted by rk4417 - phase2 +REAL(r_2) :: canopy_sublayer_dz(mp) ! inserted by rk4417 - phase2 +REAL(r_2) :: rt_Dff ! inserted by rk4417 - phase2 + !local vars INTEGER :: j REAL, DIMENSION(mp) :: q_air @@ -129,39 +139,59 @@ FUNCTION Humidity_deficit_method( mp, Ctfrz, veg_clitt,cable_user_or_evap, & ENDIF ENDDO -IF (cable_user_or_evap .or. cable_user_gw_model) then - write(6,*) "GW or ORevepis not an option right now" - !H! IF (cable_user_or_evap) THEN - !H! do j=1,mp - !H! - !H! if (veg_iveg(j) .lt. 16 .and. ssnow_snowd(j) .lt. 1e-7) THEN - !H! - !H! if (dq(j) .le. 0.0) THEN - !H! ssnow_rtevap_sat(j) = min(rtevap_max,canopy_sublayer_dz(j)/rt_Dff) - !H! end if - !H! - !H! if (dqu(j) .le. 0.0) THEN - !H! ssnow_rtevap_unsat(j) = min(rtevap_max,canopy_sublayer_dz(j)/rt_Dff) - !H! end if - !H! - !H! end if - !H! - !H! end do - !H! - !H! END IF - - ssnowpotev = air_rho * air_rlam * ( & +!IF (cable_user_or_evap .or. cable_user_gw_model) then +! write(6,*) "GW or ORevepis not an option right now" +! !H! IF (cable_user_or_evap) THEN +! !H! do j=1,mp +! !H! +! !H! if (veg_iveg(j) .lt. 16 .and. ssnow_snowd(j) .lt. 1e-7) THEN +! !H! +! !H! if (dq(j) .le. 0.0) THEN +! !H! ssnow_rtevap_sat(j) = min(rtevap_max,canopy_sublayer_dz(j)/rt_Dff) +! !H! end if +! !H! +! !H! if (dqu(j) .le. 0.0) THEN +! !H! ssnow_rtevap_unsat(j) = min(rtevap_max,canopy_sublayer_dz(j)/rt_Dff) +! !H! end if +! !H! +! !H! end if +! !H! +! !H! end do +! !H! +! !H! END IF + +! block above replaced by below - rk4417 - phase2 + +IF (cable_user_or_evap) THEN ! .or. cable_user_gw_model) then ! MMY + + DO j=1,mp + + IF (veg_iveg(j) .LT. 16 .AND. ssnow_snowd(j) .LT. 1e-7) THEN + + IF (dq(j) .LE. 0.0) THEN + ssnow_rtevap_sat(j) = MIN(rtevap_max,canopy_sublayer_dz(j)/rt_Dff) + END IF + + IF (dqu(j) .LE. 0.0) THEN + ssnow_rtevap_unsat(j) = MIN(rtevap_max,canopy_sublayer_dz(j)/rt_Dff) + END IF + + END IF + + END DO + + ssnowpotev = air_rho * air_rlam * ( & REAL(ssnow_satfrac) * dq /(ssnow_rtsoil + REAL(ssnow_rtevap_sat)) + & (1.0 - REAL(ssnow_satfrac))* dqu/( & ssnow_rtsoil + REAL(ssnow_rtevap_unsat)) ) - ELSEIF (cable_user_litter) THEN + ELSE IF (cable_user_litter) THEN ! vh_js ! - ssnowpotev = air_rho * air_rlam * dq /( ssnow_rtsoil + & + ssnowpotev = air_rho * air_rlam * dq /( ssnow_rtsoil + & REAL((1-ssnow_isflag))* veg_clitt*0.003/canopy_DvLitt) - ELSE - ssnowpotev = air_rho * air_rlam * dq / ssnow_rtsoil - ENDIF + ELSE + ssnowpotev = air_rho * air_rlam * dq / ssnow_rtsoil +END IF RETURN END FUNCTION Humidity_deficit_method diff --git a/src/science/canopy/cbl_wetleaf.F90 b/src/science/canopy/cbl_wetleaf.F90 index a2f223b9e..08f7d1262 100644 --- a/src/science/canopy/cbl_wetleaf.F90 +++ b/src/science/canopy/cbl_wetleaf.F90 @@ -56,7 +56,7 @@ SUBROUTINE wetLeaf( dels, cansat, tlfy, & !i sums, terms of convenience/readability REAL, DIMENSION(mp) :: & - sum_gbh, xx1 + sum_gbh !, xx1 ! xx1 not used - rk4417 - phase2 INTEGER :: j From 8069df269aca7c5e55cf9d087003d430615166eb Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:49:17 +1000 Subject: [PATCH 24/28] Done science/gw_hydro --- src/science/gw_hydro/cable_psm.F90 | 59 +++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 13 deletions(-) diff --git a/src/science/gw_hydro/cable_psm.F90 b/src/science/gw_hydro/cable_psm.F90 index ac4d201e1..ff0dda8b2 100644 --- a/src/science/gw_hydro/cable_psm.F90 +++ b/src/science/gw_hydro/cable_psm.F90 @@ -5,6 +5,8 @@ MODULE cable_psm roughness_type USE cable_common_module, ONLY : cable_user +! USE cable_surface_types_mod, ONLY : lakes_cable ! should include this line to replace 16 everywhere below with lakes_cable - rk4417 + IMPLICIT NONE @@ -14,6 +16,8 @@ MODULE cable_psm litter_thermal_diff=2.7e-5 !param based on vh thermal diffusivity REAL(r_2), PARAMETER :: rtevap_max = 10000.0 + ! these precomputed values are taken by the sample code in Wikipedia, + ! and the sample itself takes them from the GNU Scientific Library ! comment inserted by rk4417 - phase2 REAL(r_2), DIMENSION(0:8), PARAMETER :: gamma_pre = & (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, & 771.32342877765313, -176.61502916214059, 12.507343278686905, & @@ -69,7 +73,8 @@ SUBROUTINE or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) - REAL(r_2), DIMENSION(mp) :: sublayer_dz, eddy_shape,eddy_mod,soil_moisture_mod, & +! REAL(r_2), DIMENSION(mp) :: sublayer_dz, eddy_shape,eddy_mod,soil_moisture_mod, & + REAL(r_2), DIMENSION(mp) :: eddy_shape,eddy_mod,soil_moisture_mod, & ! sublayer_dz is not used - rk4417 - phase2 soil_moisture_mod_sat, wb_liq, & pore_size,pore_radius, rel_s,hk_zero,hk_zero_sat,time_scale !note pore_size in m @@ -81,22 +86,46 @@ SUBROUTINE or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) INTEGER :: i,j,k - IF (cable_user%litter) THEN - litter_dz(:) = veg%clitt*0.003 - ELSE - litter_dz(:) = 0.0 - ENDIF +! canopy%sublayer_dz(:) = 0.005 ! this line appears in MMY code but will leave commented for now -- rk4417 ! MMY@23Apr2023 need to check whether it should be commented + +! IF (cable_user%litter) THEN +! litter_dz(:) = veg%clitt*0.003 +! ELSE +! litter_dz(:) = 0.0 +! ENDIF + +! replaced above block by below - rk4417 - phase2 + + litter_dz(:) = 0.0 + IF (cable_user%litter) THEN + WHERE (ssnow%isflag == 0 .OR. ssnow%snowd <= 0.1) + litter_dz(:) = veg%clitt*0.003 + END WHERE + END IF pore_radius(:) = 0.148 / (1000.0*9.81*ABS(soil%sucs_vec(:,1))/1000.0) !should replace 0.148 with surface tension, unit coversion, and angle pore_size(:) = pore_radius(:)*SQRT((pi_r_2)) !scale ustar according to the exponential wind profile, assuming we are a mm from the surface - eddy_shape = 0.3*met%ua/ MAX(1.0e-4,MAX(1.0e-3,canopy%us)*& - EXP(-rough%coexp*(1.0-canopy%sublayer_dz/MAX(1e-2,rough%hruff)))) - int_eddy_shape = FLOOR(eddy_shape) - eddy_mod(:) = 0.0 + +! eddy_shape = 0.3*met%ua/ MAX(1.0e-4,MAX(1.0e-3,canopy%us)*& +! EXP(-rough%coexp*(1.0-canopy%sublayer_dz/MAX(1e-2,rough%hruff)))) +! int_eddy_shape = FLOOR(eddy_shape) +! eddy_mod(:) = 0.0 + +! replaced above block by below - rk4417 - phase2 + + eddy_mod(:) = 0.0 + eddy_shape(:) = 1.0 + int_eddy_shape(:) = 0 + DO i=1,mp - IF (veg%iveg(i) .LT. 16) THEN + IF (veg%iveg(i) .LT. 16) THEN ! should probably replace 16 with lakes_cable - rk4417 +! 2 statements below inserted by rk4417 - phase2 + eddy_shape(i) = 0.3*met%ua(i)/ max(1.0e-4,max(1.0e-3,canopy%us(i))*& + exp(-rough%coexp(i)*(1.0-canopy%sublayer_dz(i)/max(1e-2,rough%hruff(i))))) + int_eddy_shape(i) = floor(eddy_shape(i)) + eddy_mod(i) = 2.2*SQRT(112.0*(pi_r_2)) / (2.0**(eddy_shape(i)+1.0) * SQRT(eddy_shape(i)+1.0)) IF (int_eddy_shape(i) .GT. 0) THEN @@ -117,8 +146,9 @@ SUBROUTINE or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) END DO DO i=1,mp - IF (veg%iveg(i) .LT. 16) THEN + IF (veg%iveg(i) .LT. 16) THEN ! should probably replace 16 with lakes_cable - rk4417 + ! ccc Should include den_rat before ssnow%wbice for gw_model. wb_liq(i) = REAL(MAX(0.0001,MIN((pi_r_2)/4.0, & (ssnow%wb(i,1)-ssnow%wbice(i,1) - ssnow%satfrac(i)*soil%ssat_vec(i,1)) / & MAX((1._r_2 - ssnow%satfrac(i)),1e-5) ) ) ) @@ -145,6 +175,8 @@ SUBROUTINE or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) END IF IF (canopy%sublayer_dz(i) .GE. 1.0e-7) THEN +! ccc Possible other formulation? Found in MMY groundwater code branch. +! if (canopy%sublayer_dz(i) .ge. 1.0e-7 .and. hk_zero(i) .lt. 1.0e14) then ssnow%rtevap_unsat(i) = MIN(rtevap_max, & rough%z0soil(i)/canopy%sublayer_dz(i) * (lm/ (4.0*hk_zero(i)) +& (canopy%sublayer_dz(i) + pore_size(i) * soil_moisture_mod(i)) / rt_Dff)) @@ -161,6 +193,7 @@ SUBROUTINE or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) lm/ (4.0*hk_zero_sat(i)) + (canopy%sublayer_dz(i) + pore_size(i) * soil_moisture_mod_sat(i)) / rt_Dff) ssnow%rt_qh_sublayer(i) = 0.0 +! ssnow%rt_qh_sublayer(i) = canopy%sublayer_dz(i) / litter_thermal_diff ! line replaces above one in MMY code -- rk4417 ! MMY@23Apr2023 need to test, I have no idea why it is changed END IF ELSE @@ -168,7 +201,7 @@ SUBROUTINE or_soil_evap_resistance(soil,air,met,canopy,ssnow,veg,rough) ssnow%rtevap_unsat(i) = 0.0 ssnow%rt_qh_sublayer(i) = 0.0 ssnow%satfrac(i) = 0.5 - IF (veg%iveg(i) .EQ. 16 .AND. met%tk(i) .LT. 268.15 ) & + IF (veg%iveg(i) .EQ. 16 .AND. met%tk(i) .LT. 268.15 ) & ! should probably replace 16 with lakes_cable - rk4417 ssnow%rtevap_sat(i) = 0.41*ssnow%rtsoil(i) END IF From 7ffe3aef3d2d4b961bd1ddad1c24e21c7e240284 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:49:17 +1000 Subject: [PATCH 25/28] Done src/util --- src/util/cable_common.F90 | 2 ++ src/util/cable_runtime_opts_mod.F90 | 4 +++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/util/cable_common.F90 b/src/util/cable_common.F90 index 20106ae42..4f3d9299a 100644 --- a/src/util/cable_common.F90 +++ b/src/util/cable_common.F90 @@ -91,6 +91,8 @@ MODULE cable_common_module !! name of file for IGBP soil map CHARACTER(LEN=500) :: gw_elev !! name of file for gw/elevation data + CHARACTER(LEN=500) :: gw_soils + !! tiled/layered soil parameters CHARACTER(LEN=500) :: fxpft !! filename for PFT fraction and transition,wood harvest, secondary harvest CHARACTER(LEN=500) :: fxluh2cable diff --git a/src/util/cable_runtime_opts_mod.F90 b/src/util/cable_runtime_opts_mod.F90 index 8719cbc06..4e746d8c2 100644 --- a/src/util/cable_runtime_opts_mod.F90 +++ b/src/util/cable_runtime_opts_mod.F90 @@ -110,7 +110,7 @@ MODULE cable_runtime_opts_mod LOGICAL :: or_evap = .FALSE. LOGICAL :: test_new_gw = .FALSE. LOGICAL :: sync_nc_file = .FALSE. - INTEGER :: max_spins = -1 + INTEGER :: max_spins = -1 ! MMY-phase2 change from -1 to 30 - rk4417 LOGICAL :: fix_access_roots = .FALSE. !use pft dependent roots in ACCESS !ACCESS roots LOGICAL :: access13roots = .FALSE. !switch to use ACCESS1.3 %froot @@ -118,6 +118,8 @@ MODULE cable_runtime_opts_mod LOGICAL :: l_limit_labile = .FALSE. ! #237: limit Labile in spinup LOGICAL :: NtilesThruMetFile = .FALSE. ! #199: Specify Ntiles thru met file ! #338 https://github.com/CABLE-LSM/CABLE/issues/338 + ! rk4417 - phase2 + INTEGER :: force_npatches_as=-1 LOGICAL :: l_ice_consistency = .FALSE. END TYPE kbl_user_switches From 4410fc38a8d8b088912fc97381a26d578f27ef26 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:49:18 +1000 Subject: [PATCH 26/28] Done src/offline --- src/offline/cable_checks.F90 | 16 +- src/offline/cable_define_types.F90 | 3 + src/offline/cable_initialise.F90 | 88 ++++++- src/offline/cable_input.F90 | 5 + src/offline/cable_iovars.F90 | 15 +- src/offline/cable_output.F90 | 370 +++++++++++++++++++++++------ 6 files changed, 411 insertions(+), 86 deletions(-) diff --git a/src/offline/cable_checks.F90 b/src/offline/cable_checks.F90 index 5ad32517c..ac1930d47 100644 --- a/src/offline/cable_checks.F90 +++ b/src/offline/cable_checks.F90 @@ -84,7 +84,8 @@ MODULE cable_checks_module ESoil = [-0.0015, 0.0015], & TVeg = [-0.0003, 0.0003], & ECanop = [-0.0003, 0.0003], & - PotEvap = [-0.005, 0.005], & !note should encompass Evap + PotEvap = [-0.005, 0.005], & !note should encompass Evap ! I have PotEvap = [-0.0006, 0.0006] if it makes a difference - rk4417 + ACond = [0.0, 1.0], & SoilWet = [-0.4, 1.2], & Albedo = [0.0, 1.0], & @@ -133,12 +134,13 @@ MODULE cable_checks_module css = [700.0, 2200.0], & rhosoil = [300.0, 3000.0], & hyds = [5.0E-7, 8.5E-3], & ! vh_js ! sep14 +! MMY 8.5E-3->8.5 since hyds uses m/s, but hyds_vec uses mm/s rs20 = [0.0, 10.0], & sand = [0.0, 1.0], & sfc = [0.1, 0.5], & silt = [0.0, 1.0], & ssat = [0.35, 0.5], & - sucs = [-0.8, -0.03], & + sucs = [-0.8, -0.03], & ! MMY@23Apr2023 works for CABLE non-GW swilt = [0.05, 0.4], & froot = [0.0, 1.0], & zse = [0.0, 5.0], & @@ -498,7 +500,10 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met, ! Local variables REAL(r_2), DIMENSION(:),POINTER, SAVE :: owb ! volumetric soil moisture ! at previous time step + REAL(r_2), DIMENSION(mp) :: delwb ! change in soilmoisture + ! To remove qrecharge or not from the mass balance + REAL(r_2), DIMENSION(mp) :: qrecharge_opt ! b/w tsteps REAL, DIMENSION(mp) :: canopy_wbal !canopy water balance INTEGER :: j, k ! do loop counter @@ -521,9 +526,14 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met, ! it's included in change in canopy storage calculation)) ! rml 28/2/11 ! BP changed rnof1+rnof2 to ssnow%runoff which also included rnof5 ! which is used when nglacier=2 in soilsnow routines (BP feb2011) + qrecharge_opt = ssnow%qrecharge + IF ( cable_user%GW_MODEL ) THEN + qrecharge_opt = 0.0_r_2 + ENDIF + bal%wbal = REAL(met%precip - canopy%delwc - ssnow%snowd+ssnow%osnowd & - ssnow%runoff-(canopy%fevw+canopy%fevc & - + canopy%fes/ssnow%cls)*dels/air%rlam - delwb - ssnow%qrecharge) + + canopy%fes/ssnow%cls)*dels/air%rlam - delwb - qrecharge_opt) ! Canopy water balance: precip-change.can.storage-throughfall-evap+dew canopy_wbal = REAL(met%precip-canopy%delwc-canopy%through & diff --git a/src/offline/cable_define_types.F90 b/src/offline/cable_define_types.F90 index 5358b6fbe..43602ab8e 100644 --- a/src/offline/cable_define_types.F90 +++ b/src/offline/cable_define_types.F90 @@ -214,6 +214,9 @@ MODULE cable_def_types_mod TYPE soil_snow_type INTEGER, DIMENSION(:), POINTER :: isflag ! 0 => no snow 1 => snow + !* Flag (0/1) to indicate if snow is present + ! For isflag=1, if there is enough snow, it will use a 3 layer snowpack + ! If there isn't enough snow, it will use a single layer snowpack REAL, DIMENSION(:), POINTER :: & iantrct, & ! pointer to Antarctic land points diff --git a/src/offline/cable_initialise.F90 b/src/offline/cable_initialise.F90 index 766e6bcce..e4fc8e65e 100644 --- a/src/offline/cable_initialise.F90 +++ b/src/offline/cable_initialise.F90 @@ -141,6 +141,21 @@ SUBROUTINE get_default_inits(met,soil,ssnow,canopy,logn, EMSOIL) canopy%fhs = 0.0 ! sensible heat flux from soil (W/m2) canopy%us = 0.1 ! friction velocity (needed in roughness before first call to canopy: should in be in restart?) +! block below added by rk4417 - phase2 +! ccc - These need to be set elsewhere since this subroutine is deprecated +! (it is never called) and should be removed. + canopy%sublayer_dz = 0.001 + ssnow%rtevap_sat = 0.0 + ssnow%rtevap_unsat = 0.0 + ssnow%satfrac = 1.0e-12 + ssnow%qhz = 0.0 + ssnow%wtd = 1000.0 + ssnow%wb_hys = 0.99*soil%ssat_vec + ssnow%smp_hys = -0.99*soil%sucs_vec + ssnow%ssat_hys = gw_params%ssat_wet_factor*soil%ssat_vec + ssnow%watr_hys = soil%watr + ssnow%hys_fac = 1.0 + END SUBROUTINE get_default_inits !============================================================================== @@ -408,15 +423,74 @@ SUBROUTINE get_restart_data(logn,ssnow,canopy,rough,bgc, & CALL readpar(ncid_rin,'runoff',dummy,ssnow%runoff,filename%restart_in, & max_vegpatches,'def',from_restart,mp) - !MD - ok = NF90_INQ_VARID(ncid_rin,'GWwb',parID) - IF(ok == NF90_NOERR) THEN - CALL readpar(ncid_rin,'GWwb',dummy,ssnow%GWwb,filename%restart_in, & - max_vegpatches,'def',from_restart,mp) - ELSE - ssnow%GWwb = 0.95*soil%ssat + IF (cable_user%gw_model) THEN ! inserted by rk4417 - phase2 + ok = NF90_INQ_VARID(ncid_rin,'GWwb',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'GWwb',dummy,ssnow%GWwb,filename%restart_in, & + max_vegpatches,'def',from_restart,mp) + ELSE + ssnow%GWwb = 0.95*soil%ssat + END IF + +! below part until END IF inserted by rk4417 - phase2 + + ok = NF90_INQ_VARID(ncid_rin,'wb_hys',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'wb_hys',dummy,ssnow%wb_hys,filename%restart_in,& + max_vegpatches,'msd',from_restart,mp) + ELSE + ssnow%wb_hys = 0.99*soil%ssat_vec + END IF + + ok = NF90_INQ_VARID(ncid_rin,'smp_hys',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'smp_hys',dummy,ssnow%smp_hys,filename%restart_in,& + max_vegpatches,'msd',from_restart,mp) + ELSE + ssnow%smp_hys = -1.0*abs(soil%sucs_vec)*0.99 + END IF + + ok = NF90_INQ_VARID(ncid_rin,'ssat_hys',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'ssat_hys',dummy,ssnow%ssat_hys,filename%restart_in,& + max_vegpatches,'msd',from_restart,mp) + ELSE + ssnow%ssat_hys = soil%ssat_vec + END IF + + ok = NF90_INQ_VARID(ncid_rin,'watr_hys',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'watr_hys',dummy,ssnow%watr_hys,filename%restart_in,& + max_vegpatches,'msd',from_restart,mp) + ELSE + ssnow%watr_hys = soil%watr + END IF + + + ok = NF90_INQ_VARID(ncid_rin,'hys_fac',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'hys_fac',dummy,ssnow%hys_fac,filename%restart_in,& + max_vegpatches,'msd',from_restart,mp) + ELSE + ssnow%hys_fac = 1.0 + END IF END IF + IF (cable_user%or_evap) THEN + ok = NF90_INQ_VARID(ncid_rin,'sublayer_dz',parID) + IF(ok == NF90_NOERR) THEN + CALL readpar(ncid_rin,'sublayer_dz',dummy,canopy%sublayer_dz,filename%restart_in,& + max_vegpatches,'def',from_restart,mp) + ELSE + canopy%sublayer_dz(:) = 0.01 + END IF + + IF (ANY(canopy%sublayer_dz < 0.0) .OR. ANY(canopy%sublayer_dz > 0.5)) THEN + WRITE(*,*) 'problem with sublayer_dz and restart. check restart values!' + END IF + END IF + + !$ IF(cable_user%SOIL_STRUC=='sli'.or.cable_user%FWSOIL_SWITCH=='Haverd2013') THEN !$ CALL readpar(ncid_rin,'gamma',dummy,veg%gamma,filename%restart_in, & !$ max_vegpatches,'def',from_restart,mp) diff --git a/src/offline/cable_input.F90 b/src/offline/cable_input.F90 index 5f0aa20f7..5325eba37 100644 --- a/src/offline/cable_input.F90 +++ b/src/offline/cable_input.F90 @@ -2832,6 +2832,11 @@ SUBROUTINE load_parameters(met,air,ssnow,veg,climate,bgc,soil,canopy,rough,rad, sum_flux,veg,mp) WRITE(logn,*) ' CABLE variables allocated with ', mp, ' patch(es).' +! call below inserted by rk4417 - phase2 + !CALL for gw_model false and true sets constants when false + CALL GWspatialParameters(logn,soil,ssnow) ! MMY gw_model = True read var from gridinfo + ! MMY gw_model = False use default values + IF (icycle > 0 .OR. CABLE_USER%CASA_DUMP_WRITE ) & CALL alloc_casavariable(casabiome,casapool,casaflux, & casamet,casabal,mp) diff --git a/src/offline/cable_iovars.F90 b/src/offline/cable_iovars.F90 index a2ca5247a..f2ce83915 100644 --- a/src/offline/cable_iovars.F90 +++ b/src/offline/cable_iovars.F90 @@ -160,7 +160,9 @@ MODULE cable_IO_vars_module veg_class,soil_class,mvtype,mstype,patchfrac, & WatSat,GWWatSat,SoilMatPotSat,GWSoilMatPotSat, & HkSat,GWHkSat,FrcSand,FrcClay,Clappb,Watr,GWWatr,sfc_vec,forg,swilt_vec, & - slope,slope_std,GWdz,SatFracmax,Qhmax,QhmaxEfold,HKefold,HKdepth + slope,slope_std,GWdz,SatFracmax,Qhmax,QhmaxEfold,HKefold,HKdepth, & + sand_vec,clay_vec,bch_vec,org_vec,elev,elev_std ! inserted by rk4417 - phase2 + INTEGER :: ishorizon,nhorizons,clitt, & zeta,fsatmax, & gamma,ZR,F10 @@ -382,6 +384,8 @@ MODULE cable_IO_vars_module isoil = .FALSE., & ! soil type from global index meth = .FALSE., & ! method for solving turbulence in canopy scheme za = .FALSE., & ! something to do with roughness ???? + elev = .false.,& !mean subgrid elev ! inserted 2 lines - rk4417 - phase2 + elev_std=.false.,& !stddev of subgrid elev slope = .FALSE.,& !mean subgrid slope slope_std=.FALSE.,& !stddev of subgrid slope GWdz=.FALSE.,& !aquifer thickness @@ -389,7 +393,14 @@ MODULE cable_IO_vars_module Qhmax=.FALSE.,& QhmaxEfold=.FALSE.,& HKefold=.FALSE.,& - HKdepth +! HKdepth ! commented out by rk4417 - phase2 + HKdepth=.false.,& ! added this block - rk4417 - phase2 + SMP=.false.,& + SMP_hys=.false.,& + WB_hys=.false.,& + SSAT_hys=.false.,& + WATR_hys=.false.,& + hys_fac=.false. END TYPE output_inclusion_type diff --git a/src/offline/cable_output.F90 b/src/offline/cable_output.F90 index b7ada0c57..f8b5600a6 100644 --- a/src/offline/cable_output.F90 +++ b/src/offline/cable_output.F90 @@ -69,7 +69,9 @@ MODULE cable_output_module PlantTurnover, PlantTurnoverLeaf, PlantTurnoverFineRoot, & PlantTurnoverWood, PlantTurnoverWoodDist, PlantTurnoverWoodCrowding, & PlantTurnoverWoodResourceLim, dCdt, Area, LandUseFlux, patchfrac, & - vcmax,hc,WatTable,GWMoist,SatFrac,Qrecharge + vcmax,hc,WatTable,GWMoist,SatFrac,Qrecharge, & + SMP,SMP_hys,WB_hys,SSAT_hys, WATR_hys,hys_fac ! added by rk4417 - phase2 + END TYPE out_varID_type TYPE(out_varID_type) :: ovid ! netcdf variable IDs for output variables TYPE(parID_type) :: opid ! netcdf variable IDs for output variables @@ -229,6 +231,17 @@ MODULE cable_output_module REAL(KIND=4), POINTER, DIMENSION(:) :: RootResp ! autotrophic root respiration [umol/m2/s] REAL(KIND=4), POINTER, DIMENSION(:) :: StemResp ! autotrophic stem respiration [umol/m2/s] + +! remainder of TYPE block below added by rk4417 - phase2 + + REAL(KIND=4), POINTER, DIMENSION(:,:) :: SMP ! soil pressure [m] + REAL(KIND=4), POINTER, DIMENSION(:,:) :: SMP_hys ! soil pressure [m] + REAL(KIND=4), POINTER, DIMENSION(:,:) :: WB_hys ! soil pressure [m] + + REAL(KIND=4), POINTER, DIMENSION(:,:) :: SSAT_hys ! soil pressure [m] + REAL(KIND=4), POINTER, DIMENSION(:,:) :: WATR_hys ! soil pressure [m] + REAL(KIND=4), POINTER, DIMENSION(:,:) :: hys_fac ! soil pressure [m] + END TYPE output_temporary_type TYPE output_var_settings_type @@ -627,6 +640,44 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) ALLOCATE(out%SoilTemp(mp,ms)) out%SoilTemp = 0.0 ! initialise END IF + +! block below inserted by rk4417 - phase2 + IF(output%soil .OR. output%SMP) THEN + CALL define_ovar(ncid_out, ovid%SMP, 'SMP', 'm', & + 'Average layer soil pressure', patchout%SMP, & + 'soil', xID, yID, zID, landID, patchID, soilID, tID) + ALLOCATE(out%SMP(mp,ms)) + out%SMP = 0.0 ! initialise + ENDIF + if (cable_user%gw_model .and. gw_params%BC_hysteresis) then + CALL define_ovar(ncid_out, ovid%SMP_hys, 'SMP_hys', 'm', & + 'Average layer soil pressure at hys trans', patchout%SMP_hys, & + 'soil', xID, yID, zID, landID, patchID, soilID, tID) + ALLOCATE(out%SMP_hys(mp,ms)) + out%SMP_hys = 0.0 ! initialise + CALL define_ovar(ncid_out, ovid%WB_hys, 'WB_hys', 'm', & + 'wb at wet/dry or dry/wet transition', patchout%WB_hys, & + 'soil', xID, yID, zID, landID, patchID, soilID, tID) + ALLOCATE(out%WB_hys(mp,ms)) + out%WB_hys = 0.0 ! initialise + CALL define_ovar(ncid_out, ovid%SSAT_hys, 'SSAT_hys', 'm', & + 'hysteresis adj ssat', patchout%SSAT_hys, & + 'soil', xID, yID, zID, landID, patchID, soilID, tID) + ALLOCATE(out%SSAT_hys(mp,ms)) + out%SSAT_hys = 0.0 ! initialise + CALL define_ovar(ncid_out, ovid%WATR_hys, 'WATR_hys', 'm', & + 'hysteresis adj watr', patchout%WATR_hys, & + 'soil', xID, yID, zID, landID, patchID, soilID, tID) + ALLOCATE(out%WATR_hys(mp,ms)) + out%WATR_hys = 0.0 ! initialise + CALL define_ovar(ncid_out, ovid%hys_fac, 'hys_fac', 'm', & + '1.0 wet 0.5 dry', patchout%hys_fac, & + 'soil', xID, yID, zID, landID, patchID, soilID, tID) + ALLOCATE(out%hys_fac(mp,ms)) + out%hys_fac = 0.0 ! initialise + END IF +! end of block - rk4417 - phase2 + IF(output%BaresoilT) THEN CALL define_ovar(ncid_out, ovid%BaresoilT, 'BaresoilT', & 'K', 'Bare soil temperature', patchout%BaresoilT, & @@ -1096,39 +1147,121 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) IF(output%isoil) CALL define_ovar(ncid_out, opid%isoil, & 'isoil', '-', 'Soil type', patchout%isoil, 'integer', & xID, yID, zID, landID, patchID) - IF(output%bch) CALL define_ovar(ncid_out, opid%bch, & - 'bch', '-', 'Parameter b, Campbell eqn 1985', patchout%bch, 'real', & + IF(output%bch) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%bch, & + 'bch', '-', 'Parameter b, Campbell eqn 1985', patchout%bch, soilID,'soil', & + xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%bch, & + 'bch', '-', 'Parameter b, Campbell eqn 1985', patchout%bch, 'real', & + xID, yID, zID, landID, patchID) + END IF + END IF + IF(output%clay) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%clay, & + 'clay', '-', 'Fraction of soil which is clay', patchout%clay, soilID,'soil', & xID, yID, zID, landID, patchID) - IF(output%clay) CALL define_ovar(ncid_out, opid%clay, & + ELSE + CALL define_ovar(ncid_out, opid%clay, & 'clay', '-', 'Fraction of soil which is clay', patchout%clay, 'real', & xID, yID, zID, landID, patchID) - IF(output%sand) CALL define_ovar(ncid_out, opid%sand, & + END IF + END IF + IF(output%sand) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%sand, & + 'sand', '-', 'Fraction of soil which is sand', patchout%sand, soilID,'soil', & + xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%sand, & 'sand', '-', 'Fraction of soil which is sand', patchout%sand, 'real', & xID, yID, zID, landID, patchID) - IF(output%silt) CALL define_ovar(ncid_out, opid%silt, & - 'silt', '-', 'Fraction of soil which is silt', patchout%silt, 'real', & + END IF + END IF + IF(output%silt) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%silt, & + 'silt', '-', 'Fraction of soil which is silt', patchout%silt, soilID,'soil', & xID, yID, zID, landID, patchID) - IF(output%ssat) CALL define_ovar(ncid_out, opid%ssat, & - 'ssat', '-', 'Fraction of soil volume which is water @ saturation', & - patchout%ssat, 'real', xID, yID, zID, landID, patchID) - IF(output%sfc) CALL define_ovar(ncid_out, opid%sfc, & - 'sfc', '-', 'Fraction of soil volume which is water @ field capacity', & - patchout%sfc, 'real', xID, yID, zID, landID, patchID) - IF(output%swilt) CALL define_ovar(ncid_out, opid%swilt, & - 'swilt', '-', 'Fraction of soil volume which is water @ wilting point', & - patchout%swilt, 'real', xID, yID, zID, landID, patchID) - IF(output%hyds) CALL define_ovar(ncid_out, opid%hyds, & - 'hyds', 'm/s', 'Hydraulic conductivity @ saturation', & - patchout%hyds, 'real', xID, yID, zID, landID, patchID) - IF(output%sucs) CALL define_ovar(ncid_out, opid%sucs, & - 'sucs', 'm', 'Suction @ saturation', & + ELSE + CALL define_ovar(ncid_out, opid%silt, & + 'silt', '-', 'Fraction of soil which is silt', patchout%silt, 'real', & + xID, yID, zID, landID, patchID) + END IF + END IF + IF(output%ssat) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%ssat, & + 'ssat', '-', 'Fraction of soil volume which is water @ saturation', & + patchout%ssat, soilID,'soil', xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%ssat, & + 'ssat', '-', 'Fraction of soil volume which is water @ saturation', & + patchout%ssat, 'real', xID, yID, zID, landID, patchID) + END IF + END IF + IF(output%sfc) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%sfc, & + 'sfc', '-', 'Fraction of soil volume which is water @ field capacity', & + patchout%sfc, soilID,'soil', xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%sfc, & + 'sfc', '-', 'Fraction of soil volume which is water @ field capacity', & + patchout%sfc, 'real', xID, yID, zID, landID, patchID) + END IF + END IF + IF(output%swilt) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%swilt, & + 'swilt', '-', 'Fraction of soil volume which is water @ wilting point', & + patchout%swilt, soilID,'soil', xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%swilt, & + 'swilt', '-', 'Fraction of soil volume which is water @ wilting point', & + patchout%swilt, 'real', xID, yID, zID, landID, patchID) + END IF + END IF + IF(output%hyds) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, opid%hyds, & + 'hyds', 'm/s', 'Hydraulic conductivity @ saturation', & + patchout%hyds, soilID,'soil', xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%hyds, & + 'hyds', 'm/s', 'Hydraulic conductivity @ saturation', & + patchout%hyds, 'real', xID, yID, zID, landID, patchID) + END IF + END IF + IF(output%sucs) THEN + IF (cable_user%gw_model) THEN + ! gw_model uses sucs_vec that is a 2D variable + CALL define_ovar(ncid_out, opid%sucs, & + 'sucs', 'm', 'Suction @ saturation', & + patchout%sucs, soilID,'soil', xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, opid%sucs, & + 'sucs', 'm', 'Suction @ saturation', & patchout%sucs, 'real', xID, yID, zID, landID, patchID) + END IF + END IF IF(output%css) CALL define_ovar(ncid_out, opid%css, & 'css', 'J/kg/C', 'Heat capacity of soil minerals', & - patchout%css, 'real', xID, yID, zID, landID, patchID) - IF(output%rhosoil) CALL define_ovar(ncid_out, & - opid%rhosoil, 'rhosoil', 'kg/m^3', 'Density of soil minerals', & - patchout%rhosoil, 'real', xID, yID, zID, landID, patchID) +! patchout%css, 'real', xID, yID, zID, landID, patchID) + patchout%css, soilID,'soil', xID, yID, zID, landID, patchID) + IF(output%rhosoil) THEN + IF (cable_user%gw_model) THEN + CALL define_ovar(ncid_out, & + opid%rhosoil, 'rhosoil', 'kg/m^3', 'Density of soil minerals', & + patchout%rhosoil, soilID,'soil', xID, yID, zID, landID, patchID) + ELSE + CALL define_ovar(ncid_out, & + opid%rhosoil, 'rhosoil', 'kg/m^3', 'Density of soil minerals', & + patchout%rhosoil, 'real', xID, yID, zID, landID, patchID) + END IF + END IF IF(output%rs20) CALL define_ovar(ncid_out, opid%rs20, & 'rs20', '-', 'Soil respiration coefficient at 20C', & patchout%rs20, 'real', xID, yID, zID, landID, patchID) @@ -1245,6 +1378,24 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) ! patchout%GWdz, 'real', xID, yID, zID, landID, patchID) ! IF(output%params .AND. cable_user%gw_model) THEN + +! block below inserted by rk4417 - phase2 + call define_ovar(ncid_out, opid%slope, & + 'slope', '-', 'mean subgrid topographic slope', & + patchout%slope, 'real', xid, yid, zid, landid, patchid) + call define_ovar(ncid_out, opid%elev, & + 'elev', '-', 'mean subgrid topographic elev', & + patchout%elev, 'real', xid, yid, zid, landid, patchid) + + CALL define_ovar(ncid_out, opid%slope_std, & + 'slope_std', '-', 'Mean subgrid topographic slope_std', & + patchout%slope_std, 'real', xID, yID, zID, landID, patchID) + + CALL define_ovar(ncid_out, opid%GWdz, & + 'GWdz', '-', 'Mean aquifer layer thickness ', & + patchout%GWdz, 'real', xID, yID, zID, landID, patchID) +! end of block - rk4417 - phase2 + CALL define_ovar(ncid_out, opid%Qhmax, & 'Qhmax', 'mm/s', 'Maximum subsurface drainage ', & patchout%Qhmax, 'real', xID, yID, zID, landID, patchID) @@ -1376,56 +1527,120 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) CALL check_and_write(opid%isoil, & 'isoil', REAL(soil%isoilm, 4), ranges%isoil, patchout%isoil, out_settings) END IF - out_settings%dimswitch = "real" - IF (output%bch) THEN - CALL check_and_write(opid%bch, & - 'bch', REAL(soil%bch, 4), ranges%bch, patchout%bch, out_settings) - END IF - IF (output%clay) THEN - CALL check_and_write(opid%clay, & - 'clay', REAL(soil%clay, 4), ranges%clay, patchout%clay, out_settings) - END IF - IF (output%sand) THEN - CALL check_and_write(opid%sand, & - 'sand', REAL(soil%sand, 4), ranges%sand, patchout%sand, out_settings) - END IF - IF (output%silt) THEN - CALL check_and_write(opid%silt, & - 'silt', REAL(soil%silt, 4), ranges%silt, patchout%silt, out_settings) - END IF - IF (output%css) THEN - CALL check_and_write(opid%css, & - 'css', REAL(soil%css, 4), ranges%css, patchout%css, out_settings) - END IF - IF (output%rhosoil) THEN - CALL check_and_write(opid%rhosoil, 'rhosoil',REAL(soil%rhosoil,4), & - ranges%rhosoil, patchout%rhosoil, out_settings) - END IF - IF (output%hyds) THEN - CALL check_and_write(opid%hyds, & - 'hyds', REAL(soil%hyds, 4), ranges%hyds, patchout%hyds, out_settings) - END IF - IF (output%sucs) THEN - CALL check_and_write(opid%sucs, & - 'sucs', REAL(soil%sucs, 4), ranges%sucs, patchout%sucs, out_settings) - END IF + + + IF (cable_user%gw_model) THEN + out_settings%dimswitch = "soil" + + IF (output%bch) THEN + CALL check_and_write(opid%bch, & + 'bch', REAL(soil%bch_vec, 4), ranges%bch, patchout%bch, & + out_settings) + END IF + IF (output%clay) THEN + CALL check_and_write(opid%clay, & + 'clay', REAL(soil%clay_vec, 4), ranges%clay, patchout%clay, & + out_settings) + END IF + IF (output%sand) THEN + CALL check_and_write(opid%sand, & + 'sand', REAL(soil%sand_vec, 4), ranges%sand, patchout%sand, & + out_settings) + END IF + IF (output%silt) THEN + CALL check_and_write(opid%silt, & + 'silt', REAL(soil%silt_vec, 4), ranges%silt, patchout%silt, & + out_settings) + END IF + IF (output%css) THEN + CALL check_and_write(opid%css, & + 'css', REAL(soil%css_vec, 4), ranges%css, patchout%css, & + out_settings) + END IF + IF (output%rhosoil) THEN + CALL check_and_write(opid%rhosoil, & + 'rhosoil', REAL(soil%rhosoil_vec, 4), ranges%rhosoil, & + patchout%rhosoil, out_settings) + END IF + IF (output%hyds) THEN + CALL check_and_write(opid%hyds, & + 'hyds', REAL(soil%hyds_vec, 4), ranges%hyds, patchout%hyds, & + out_settings) + END IF + IF (output%sucs) THEN + CALL check_and_write(opid%sucs, & + 'sucs', REAL(soil%sucs_vec, 4), ranges%sucs, patchout%sucs, & + out_settings) + END IF + IF (output%ssat) THEN + CALL check_and_write(opid%ssat, & + 'ssat', REAL(soil%ssat_vec, 4), ranges%ssat, patchout%ssat, & + out_settings) + END IF + IF (output%sfc) THEN + CALL check_and_write(opid%sfc, & + 'sfc', REAL(soil%sfc_vec, 4), ranges%sfc, patchout%sfc, & + out_settings) + END IF + IF (output%swilt) THEN + CALL check_and_write(opid%swilt, & + 'swilt', REAL(soil%swilt_vec, 4), ranges%swilt, patchout%swilt, & + out_settings) + END IF + ELSE + out_settings%dimswitch = "real" + IF (output%bch) THEN + CALL check_and_write(opid%bch, & + 'bch', REAL(soil%bch, 4), ranges%bch, patchout%bch, out_settings) + END IF + IF (output%clay) THEN + CALL check_and_write(opid%clay, & + 'clay', REAL(soil%clay, 4), ranges%clay, patchout%clay, out_settings) + END IF + IF (output%sand) THEN + CALL check_and_write(opid%sand, & + 'sand', REAL(soil%sand, 4), ranges%sand, patchout%sand, out_settings) + END IF + IF (output%silt) THEN + CALL check_and_write(opid%silt, & + 'silt', REAL(soil%silt, 4), ranges%silt, patchout%silt, out_settings) + END IF + IF (output%css) THEN + CALL check_and_write(opid%css, & + 'css', REAL(soil%css, 4), ranges%css, patchout%css, out_settings) + END IF + IF (output%rhosoil) THEN + CALL check_and_write(opid%rhosoil, 'rhosoil',REAL(soil%rhosoil,4), & + ranges%rhosoil, patchout%rhosoil, out_settings) + END IF + IF (output%hyds) THEN + CALL check_and_write(opid%hyds, & + 'hyds', REAL(soil%hyds, 4), ranges%hyds, patchout%hyds, out_settings) + END IF + IF (output%sucs) THEN + CALL check_and_write(opid%sucs, & + 'sucs', REAL(soil%sucs, 4), ranges%sucs, patchout%sucs, out_settings) + END IF + IF (output%ssat) THEN + CALL check_and_write(opid%ssat, & + 'ssat', REAL(soil%ssat, 4), ranges%ssat, patchout%ssat, out_settings) + END IF + IF (output%sfc) THEN + CALL check_and_write(opid%sfc, & + 'sfc', REAL(soil%sfc, 4), ranges%sfc, patchout%sfc, out_settings) + END IF + IF (output%swilt) THEN + CALL check_and_write(opid%swilt, & + 'swilt', REAL(soil%swilt, 4), ranges%swilt, patchout%swilt, out_settings) + END IF + END IF + + ! ccc - Currently the gw code does not include a vectorised rs20. IF (output%rs20) THEN CALL check_and_write(opid%rs20, & 'rs20', REAL(veg%rs20, 4), ranges%rs20, patchout%rs20, out_settings) ! 'rs20',REAL(soil%rs20,4),ranges%rs20,patchout%rs20,out_settings) END IF - IF (output%ssat) THEN - CALL check_and_write(opid%ssat, & - 'ssat', REAL(soil%ssat, 4), ranges%ssat, patchout%ssat, out_settings) - END IF - IF (output%sfc) THEN - CALL check_and_write(opid%sfc, & - 'sfc', REAL(soil%sfc, 4), ranges%sfc, patchout%sfc, out_settings) - END IF - IF (output%swilt) THEN - CALL check_and_write(opid%swilt, & - 'swilt', REAL(soil%swilt, 4), ranges%swilt, patchout%swilt, out_settings) - END IF ! IF (output%slope) THEN ! CALL check_and_write(ncid_out, opid%slope, & @@ -1447,10 +1662,17 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) END IF IF (output%zse) THEN - out_settings%dimswitch = "soil" - CALL check_and_write(opid%zse, & - 'zse', SPREAD(REAL(soil%zse, 4), 1, mp),ranges%zse, & - patchout%zse, out_settings)! no spatial dim at present + IF (cable_user%GW_MODEL) THEN + out_settings%dimswitch = "soil" + CALL check_and_write(opid%zse, & + 'zse', REAL(soil%zse_vec, 4),ranges%zse, & + patchout%zse, out_settings)! no spatial dim at present + ELSE + out_settings%dimswitch = "soil" + CALL check_and_write(opid%zse, & + 'zse', SPREAD(REAL(soil%zse, 4), 1, mp),ranges%zse, & + patchout%zse, out_settings)! no spatial dim at present + END IF END IF !~ Veg From 7d583236a6da50f70cde8c9ce1396b78e8874bda Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:50:14 +1000 Subject: [PATCH 27/28] Done src/offline --- src/offline/cable_initialise.F90 | 2 +- src/offline/cable_mpicommon.F90 | 9 +- src/offline/cable_mpimaster.F90 | 154 ++++++- src/offline/cable_mpiworker.F90 | 132 +++++- src/offline/cable_output.F90 | 106 ++++- src/offline/cable_parameters.F90 | 528 ++++++++++++----------- src/offline/cable_pft_params.F90 | 8 +- src/offline/cable_serial.F90 | 10 + src/offline/cbl_model_driver_offline.F90 | 17 +- src/science/canopy/cable_canopy.F90 | 4 +- 10 files changed, 657 insertions(+), 313 deletions(-) diff --git a/src/offline/cable_initialise.F90 b/src/offline/cable_initialise.F90 index e4fc8e65e..f12a8b9be 100644 --- a/src/offline/cable_initialise.F90 +++ b/src/offline/cable_initialise.F90 @@ -152,7 +152,7 @@ SUBROUTINE get_default_inits(met,soil,ssnow,canopy,logn, EMSOIL) ssnow%wtd = 1000.0 ssnow%wb_hys = 0.99*soil%ssat_vec ssnow%smp_hys = -0.99*soil%sucs_vec - ssnow%ssat_hys = gw_params%ssat_wet_factor*soil%ssat_vec +!ccc To come back to it ssnow%ssat_hys = gw_params%ssat_wet_factor*soil%ssat_vec ssnow%watr_hys = soil%watr ssnow%hys_fac = 1.0 diff --git a/src/offline/cable_mpicommon.F90 b/src/offline/cable_mpicommon.F90 index 5070ca57d..103f6fcd9 100644 --- a/src/offline/cable_mpicommon.F90 +++ b/src/offline/cable_mpicommon.F90 @@ -29,7 +29,8 @@ MODULE cable_mpicommon ! base number of input fields: must correspond to CALLS to ! MPI_address (field ) in *_mpimaster/ *_mpiworker - INTEGER, PARAMETER :: nparam = 340 +! INTEGER, PARAMETER :: nparam = 341 ! replaced by below by ccc - GW + INTEGER, PARAMETER :: nparam =351 ! MPI: extra params sent only if nsoilparmnew is true INTEGER, PARAMETER :: nsoilnew = 1 @@ -77,7 +78,8 @@ MODULE cable_mpicommon !INTEGER, PARAMETER :: nmat = 29 ! MPI: CABLE_r491, after following up with Bernard on the new variables ! vh sli nmat + 4 36 -> 40 - INTEGER, PARAMETER :: nmat = 39 +! INTEGER, PARAMETER :: nmat = 40 ! replaced by below by rk4417 - phase2 + INTEGER, PARAMETER :: nmat = 45 !hysteresis 41 ! MPI: number of contig vector parts / worker (results) !INTEGER, PARAMETER :: nvec = 149 @@ -94,7 +96,8 @@ MODULE cable_mpicommon ! vh sli nvec + 6 162 -> 168 ! INTEGER, PARAMETER :: nvec = 172! 168 ! INH REV_CORR +3 (SSEB +2 will be needed) - INTEGER, PARAMETER :: nvec = 175 +! INTEGER, PARAMETER :: nvec = 175 ! replaced by below by rk4417 - phase2 + INTEGER, PARAMETER :: nvec = 176! 176!175 ! MPI: number of final casa result matrices and vectors to receive ! by the master for casa_poolout and casa_fluxout diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index a4529f9b7..f7051ee47 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -684,7 +684,9 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) canopy%oldcansto=canopy%cansto ! Zero out lai where there is no vegetation acc. to veg. index - WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0. +! WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0. +! line above replaced by below - rk4417 - phase2 + WHERE ( veg%iveg(:) .GE. 14 ) veg%vlai = 0. ! MMY change from iveg%vlai to veg%vlai @@ -758,7 +760,9 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) !$ ENDIF ! Zero out lai where there is no vegetation acc. to veg. index - WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0. +! WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0. +! line above replaced by below - rk4417 - phase2 + WHERE ( veg%iveg(:) .GE. 14 ) veg%vlai = 0. ! MMY change from iveg%vlai to veg%vlai !$ ! At first time step of year, set tile area according to updated LU areas !$ IF (ktau == 1 .and. CABLE_USER%POPLUC) THEN @@ -951,7 +955,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) met%ofsd = met%fsd(:,1) + met%fsd(:,2) canopy%oldcansto=canopy%cansto - IF ( (TRIM(cable_user%MetType) .EQ. "gswp") .OR. (TRIM(cable_user%MetType) .EQ. "gswp3") ) & + IF ( (TRIM(cable_user%MetType) .EQ. "gswp") .OR. (TRIM(cable_user%MetType) .EQ. "gswp3") & + .OR. (TRIM(cable_user%MetType) .EQ. "prin") ) & !MMY CALL close_met_file IF (icycle>0 .AND. cable_user%CALL_POP) THEN @@ -1046,7 +1051,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) IF ( (.NOT. CASAONLY) .AND. spinConv ) THEN SELECT CASE (TRIM(cable_user%MetType)) - CASE ('plum', 'cru', 'gswp', 'gswp3') + CASE ('plum', 'cru', 'gswp', 'gswp3', 'prin') CALL write_output( dels, ktau_tot, met, canopy, casaflux, casapool, casamet, & ssnow, rad, bal, air, soil, veg, CSBOLTZ, CEMLEAF, CEMSOIL ) CASE DEFAULT @@ -1113,10 +1118,20 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) IF( INT( ktau_tot/kend ) > 1 ) THEN ! evaluate spinup + ! =================== MMY_phase2 + ! commented out this region in favor of the one below - rk4417 + ! ===================== IF (ANY( ABS(ssnow%wb-soilMtemp)>delsoilM).OR. & ANY(ABS(ssnow%tgg-soilTtemp)>delsoilT) .OR. & MAXVAL(ABS(ssnow%GWwb-GWtemp),dim=1)>delgwM) THEN + ! =================== MMY_phase2 uncomment ===================== +! IF( (ANY( ABS(ssnow%wb-soilMtemp)>delsoilM).OR. & +! ANY( ABS(ssnow%tgg-soilTtemp)>delsoilT) .or. & +! maxval(ABS(ssnow%GWwb-GWtemp),dim=1) > delgwM) .and. & +! ( (int(ktau_tot/kend) .lt. cable_user%max_spins) .and.& +! (cable_user%max_spins .gt. 0) ) ) THEN + ! No complete convergence yet ! PRINT *, 'ssnow%wb : ', ssnow%wb ! PRINT *, 'soilMtemp: ', soilMtemp @@ -1323,7 +1338,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) IF ( TRIM(cable_user%MetType) .NE. "gswp" .AND. & TRIM(cable_user%MetType) .NE. "gswp3" .AND. & TRIM(cable_user%MetType) .NE. "plum" .AND. & - TRIM(cable_user%MetType) .NE. "cru") CALL close_met_file + TRIM(cable_user%MetType) .NE. "cru" .AND. & + TRIM(cable_user%MetType) .NE. "prin" ) CALL close_met_file ! Close log file CLOSE(logn) @@ -2536,6 +2552,10 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Get_address (canopy%fwsoil(off), displs(bidx), ierr) blen(bidx) = r2len + bidx = bidx + 1 ! block inserted by rk4417 - phase2 + CALL MPI_Get_address (canopy%sublayer_dz(off), displs(bidx), ierr) + blen(bidx) = r2len + bidx = bidx + 1 CALL MPI_Get_address (canopy%gswx(off,1), displs(bidx), ierr) CALL MPI_Type_create_hvector (mf, r1len, r1stride, MPI_BYTE, & @@ -3124,19 +3144,36 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & & types(bidx), ierr) blen(bidx) = 1 + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%ssat_hys(off,1), displs(bidx), ierr) + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & types(bidx), ierr) + blen(bidx) = 1 bidx = bidx + 1 - CALL MPI_Get_address (soil%smpc_vec(off,1), displs(bidx), ierr) + CALL MPI_Get_address (ssnow%watr_hys(off,1), displs(bidx), ierr) CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & - & types(bidx), ierr) + & types(bidx), ierr) blen(bidx) = 1 bidx = bidx + 1 - CALL MPI_Get_address (soil%wbc_vec(off,1), displs(bidx), ierr) + CALL MPI_Get_address (ssnow%smp_hys(off,1), displs(bidx), ierr) CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & - & types(bidx), ierr) + & types(bidx), ierr) blen(bidx) = 1 + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%wb_hys(off,1), displs(bidx), ierr) + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & types(bidx), ierr) + blen(bidx) = 1 + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%hys_fac(off,1), displs(bidx), ierr) + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & types(bidx), ierr) + blen(bidx) = 1 + !1D bidx = bidx + 1 @@ -3159,14 +3196,18 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Get_address (soil%GWwatr(off), displs(bidx), ierr) blen(bidx) = r2len - bidx = bidx + 1 - CALL MPI_Get_address (soil%GWz(off), displs(bidx), ierr) - blen(bidx) = r2len + ! bidx = bidx + 1 + ! CALL MPI_Get_address (soil%GWz(off), displs(bidx), ierr) + ! blen(bidx) = r2len bidx = bidx + 1 CALL MPI_Get_address (soil%GWdz(off), displs(bidx), ierr) blen(bidx) = r2len + bidx = bidx + 1 ! inserted by rk4417 - phase2 + CALL MPI_Get_address (soil%elev(off), displs(bidx), ierr) + blen(bidx) = r2len + bidx = bidx + 1 CALL MPI_Get_address (soil%slope(off), displs(bidx), ierr) blen(bidx) = r2len @@ -3175,6 +3216,43 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Get_address (soil%slope_std(off), displs(bidx), ierr) blen(bidx) = r2len + +! block below inserted by rk4417 - phase2 + bidx = bidx + 1 + CALL MPI_Get_address (soil%drain_dens(off), displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%hkrz(off), displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%zdepth(off), displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%srf_frac_ma(off), displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%edepth_ma(off), displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%qhz_max(off), displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%qhz_efold(off), displs(bidx), ierr) + blen(bidx) = r2len +! end of block - rk4417 - phase2 + bidx = bidx + 1 CALL MPI_Get_address (ssnow%GWwb(off), displs(bidx), ierr) blen(bidx) = r2len @@ -3182,6 +3260,7 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& ! MPI: sanity check IF (bidx /= ntyp) THEN WRITE (*,*) 'master: invalid number of param_t fields ',bidx,', fix it!' + WRITE (*,*) 'local counbt bidx is ',bidx,' while ntyp is ',ntyp ! inserted by rk4417 - phase2 CALL MPI_Abort (comm, 1, ierr) END IF @@ -4404,6 +4483,7 @@ SUBROUTINE master_casa_params (comm,casabiome,casapool,casaflux,casamet,& ! MPI: sanity check IF (bidx /= ntyp) THEN WRITE (*,*) 'master: invalid number of casa_t param fields ',bidx,', fix it!' + WRITE (*,*) 'local counbt bidx is ',bidx,' while ntyp is ',ntyp ! inserted by rk4417 - phase2 CALL MPI_Abort (comm, 1, ierr) END IF @@ -4862,7 +4942,52 @@ SUBROUTINE master_outtypes (comm,met,canopy,ssnow,rad,bal,air,soil,veg) & mat_t(midx, rank), ierr) CALL MPI_Type_commit (mat_t(midx, rank), ierr) midx = midx + 1 + +! block below inserted by rk4417 - phase2 ! REAL(r_2) + CALL MPI_Get_address (ssnow%smp(off,1), maddr(midx), ierr) ! 12 + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & mat_t(midx, rank), ierr) + CALL MPI_Type_commit (mat_t(midx, rank), ierr) + midx = midx + 1 + + ! REAL(r_2) + CALL MPI_Get_address (ssnow%wb_hys(off,1), maddr(midx), ierr) ! 12 + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & mat_t(midx, rank), ierr) + CALL MPI_Type_commit (mat_t(midx, rank), ierr) + midx = midx + 1 + + ! REAL(r_2) + CALL MPI_Get_address (ssnow%smp_hys(off,1), maddr(midx), ierr) ! 12 + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & mat_t(midx, rank), ierr) + CALL MPI_Type_commit (mat_t(midx, rank), ierr) + midx = midx + 1 + + ! REAL(r_2) + CALL MPI_Get_address (ssnow%ssat_hys(off,1), maddr(midx), ierr) ! 12 + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & mat_t(midx, rank), ierr) + CALL MPI_Type_commit (mat_t(midx, rank), ierr) + midx = midx + 1 + + ! REAL(r_2) + CALL MPI_Get_address (ssnow%watr_hys(off,1), maddr(midx), ierr) ! 12 + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & mat_t(midx, rank), ierr) + CALL MPI_Type_commit (mat_t(midx, rank), ierr) + midx = midx + 1 + + ! REAL(r_2) + CALL MPI_Get_address (ssnow%hys_fac(off,1), maddr(midx), ierr) ! 12 + CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & + & mat_t(midx, rank), ierr) + CALL MPI_Type_commit (mat_t(midx, rank), ierr) + midx = midx + 1 +! end of block - rk4417 - phase2 + + ! REAL(r_1) CALL MPI_Get_address (ssnow%evapfbl(off,1), maddr(midx), ierr) ! 12 CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, & & mat_t(midx, rank), ierr) @@ -5372,6 +5497,11 @@ SUBROUTINE master_outtypes (comm,met,canopy,ssnow,rad,bal,air,soil,veg) CALL MPI_Get_address (canopy%fwsoil(off), vaddr(vidx), ierr) ! 59 blen(vidx) = cnt * extr2 + vidx = vidx + 1 ! block inserted by rk4417 - phase2 + ! REAL(r_2) + CALL MPI_Get_address (canopy%sublayer_dz(off), vaddr(vidx), ierr) ! 59 + blen(vidx) = cnt * extr2 + ! MPI: 2D vars moved above ! rwater ! evapfbl diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index 8358b21f8..1f30e5d5b 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -599,20 +599,6 @@ SUBROUTINE mpidrv_worker (comm) ENDIF - - - - - - - - - - - - - - IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)).AND. & CABLE_USER%CALL_POP) THEN @@ -1770,6 +1756,11 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Get_address (canopy%fwsoil, displs(bidx), ierr) blen(bidx) = r2len +! block below inserted by rk4417 - phase2 + bidx = bidx + 1 + CALL MPI_Get_address (canopy%sublayer_dz, displs(bidx), ierr) + blen(bidx) = r2len + bidx = bidx + 1 CALL MPI_Get_address (canopy%gswx, displs(bidx), ierr) blen(bidx) = mf * r1len @@ -2279,14 +2270,39 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& bidx = bidx + 1 CALL MPI_Get_address (soil%rhosoil_vec, displs(bidx), ierr) blen(bidx) = ms * r2len - + bidx = bidx + 1 - CALL MPI_Get_address (soil%smpc_vec, displs(bidx), ierr) + CALL MPI_Get_address (ssnow%ssat_hys, displs(bidx), ierr) blen(bidx) = ms * r2len - + + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%watr_hys, displs(bidx), ierr) + blen(bidx) = ms * r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%smp_hys, displs(bidx), ierr) + blen(bidx) = ms * r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%wb_hys, displs(bidx), ierr) + blen(bidx) = ms * r2len + + bidx = bidx + 1 - CALL MPI_Get_address (soil%wbc_vec, displs(bidx), ierr) + CALL MPI_Get_address (ssnow%hys_fac, displs(bidx), ierr) blen(bidx) = ms * r2len + ! end of block - rk4417 - phase2 + + ! bidx = bidx + 1 + ! CALL MPI_Get_address (soil%smpc_vec, displs(bidx), ierr) + ! blen(bidx) = ms * r2len + + ! bidx = bidx + 1 + ! CALL MPI_Get_address (soil%wbc_vec, displs(bidx), ierr) + ! blen(bidx) = ms * r2len !1d @@ -2310,14 +2326,20 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Get_address (soil%GWwatr, displs(bidx), ierr) blen(bidx) = r2len - bidx = bidx + 1 - CALL MPI_Get_address (soil%GWz, displs(bidx), ierr) - blen(bidx) = r2len +! commented out by rk4417 - phase2 +! bidx = bidx + 1 +! CALL MPI_Get_address (soil%GWz, displs(bidx), ierr) +! blen(bidx) = r2len bidx = bidx + 1 CALL MPI_Get_address (soil%GWdz, displs(bidx), ierr) blen(bidx) = r2len +! block below inserted by rk4417 - phase2 + bidx = bidx + 1 + CALL MPI_Get_address (soil%elev, displs(bidx), ierr) + blen(bidx) = r2len + bidx = bidx + 1 CALL MPI_Get_address (soil%slope, displs(bidx), ierr) blen(bidx) = r2len @@ -2326,6 +2348,42 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& CALL MPI_Get_address (soil%slope_std, displs(bidx), ierr) blen(bidx) = r2len +! block below inserted by rk4417 - phase2 + bidx = bidx + 1 + CALL MPI_Get_address (soil%drain_dens, displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%hkrz, displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%zdepth, displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%srf_frac_ma, displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%edepth_ma, displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%qhz_max, displs(bidx), ierr) + blen(bidx) = r2len + + + bidx = bidx + 1 + CALL MPI_Get_address (soil%qhz_efold, displs(bidx), ierr) + blen(bidx) = r2len +! end of block - rk4417 - phase2 + bidx = bidx + 1 CALL MPI_Get_address (ssnow%GWwb, displs(bidx), ierr) blen(bidx) = r2len @@ -2333,6 +2391,7 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,& ! MPI: sanity check IF (bidx /= ntyp) THEN WRITE (*,*) 'worker ',rank,' invalid number of param_t fields',bidx,', fix it!' + WRITE (*,*) 'worker ',rank,' invalid number ntyp is ',ntyp,', fix it!' ! inserted by rk4417 - phase2 CALL MPI_Abort (comm, 1, ierr) END IF @@ -3782,6 +3841,32 @@ SUBROUTINE worker_outtype (comm,met,canopy,ssnow,rad,bal,air,soil,veg) CALL MPI_Get_address (ssnow%evapfbl(off,1), displs(bidx), ierr) blocks(bidx) = r2len * ms +! block below inserted by rk4417 - phase2 + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%smp(off,1), displs(bidx), ierr) + blocks(bidx) = r2len * ms + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%wb_hys(off,1), displs(bidx), ierr) + blocks(bidx) = r2len * ms + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%smp_hys(off,1), displs(bidx), ierr) + blocks(bidx) = r2len * ms + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%ssat_hys(off,1), displs(bidx), ierr) + blocks(bidx) = r2len * ms + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%watr_hys(off,1), displs(bidx), ierr) + blocks(bidx) = r2len * ms + + bidx = bidx + 1 + CALL MPI_Get_address (ssnow%hys_fac(off,1), displs(bidx), ierr) + blocks(bidx) = r2len * ms +! end of block - rk4417 - phase2 + !midx = midx + 1 ! REAL(r_1) !CALL MPI_Get_address (ssnow%wbfice(off,1), maddr(midx), ierr) ! 13 @@ -4541,6 +4626,11 @@ SUBROUTINE worker_outtype (comm,met,canopy,ssnow,rad,bal,air,soil,veg) CALL MPI_Get_address (canopy%fwsoil(off), displs(bidx), ierr) blocks(bidx) = r2len +! block below inserted by rk4417 - phase2 + bidx = bidx + 1 + CALL MPI_Get_address (canopy%sublayer_dz(off), displs(bidx), ierr) + blocks(bidx) = r2len + ! MPI: 2D vars moved above ! rwater ! evapfbl diff --git a/src/offline/cable_output.F90 b/src/offline/cable_output.F90 index f8b5600a6..b32f75a8e 100644 --- a/src/offline/cable_output.F90 +++ b/src/offline/cable_output.F90 @@ -1800,24 +1800,29 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) out_settings%dimswitch = "real" IF (output%params .AND. cable_user%gw_model) THEN CALL check_and_write(opid%SatFracmax, & - 'SatFracmax', SPREAD(REAL(gw_params%MaxSatFraction,4),1,mp), & + 'SatFracmax', SPREAD(REAL(SQRT(gw_params%MaxSatFraction),4),1,mp), & ranges%gw_default, patchout%SatFracmax, out_settings) CALL check_and_write(opid%Qhmax, & - 'Qhmax', SPREAD(REAL(gw_params%MaxHorzDrainRate, 4),1,mp), & + 'Qhmax', REAL(soil%qhz_max, 4), & ranges%gw_default, patchout%Qhmax, out_settings) - CALL check_and_write(opid%QhmaxEfold, & - 'QhmaxEfold', SPREAD(REAL(gw_params%EfoldHorzDrainRate, 4),1,mp), & - ranges%gw_default, patchout%QhmaxEfold, out_settings) + ! part below commented out by rk4417 - phase2 + ! CALL check_and_write(opid%QhmaxEfold, & + ! 'QhmaxEfold', SPREAD(REAL(gw_params%EfoldHorzDrainRate, 4),1,mp), & + ! ranges%gw_default, patchout%QhmaxEfold, out_settings) CALL check_and_write(opid%HKefold, & - 'HKefold', SPREAD(REAL(gw_params%hkrz, 4),1,mp), & + 'HKefold', REAL(soil%hkrz, 4), & ranges%gw_default, patchout%HKefold, out_settings) + ! note specially for ranges%gw_default that I have instead + ! (/0.0,1000000.0/) - rk4417 - phase2 CALL check_and_write(opid%HKdepth, & - 'HKdepth', SPREAD(REAL(gw_params%zdepth, 4),1,mp), & + 'HKdepth', REAL(soil%zdepth, 4), & ranges%gw_default, patchout%HKdepth, out_settings) + ! note specially for ranges%gw_default that I have instead + ! (/0.0,1000000.0/) - rk4417 - phase2 END IF END SUBROUTINE open_output_file @@ -2146,6 +2151,29 @@ SUBROUTINE write_output(dels, ktau, met, canopy, casaflux, casapool, casamet, ss END IF END IF +! block below inserted by rk4417 - phase2 +! ccc - This should only be a condition on output%SMP. Need to implement it in the initialisation. + IF((output%soil .OR. output%SMP)) THEN + out_settings%dimswitch = 'soil' + CALL generate_out_write_acc(ovid%SMP, 'SMP', out%SMP, REAL(ssnow%smp, 4),& + (/-1.0e36,1.0e36/), patchout%SMP, out_settings) + END IF + + IF(gw_params%bc_hysteresis) THEN + out_settings%dimswitch = 'soil' + CALL generate_out_write_acc(ovid%SMP_hys, 'SMP_hys', out%SMP_hys, REAL(ssnow%smp_hys, 4),& + (/-1.0e36,1.0e36/), patchout%SMP_hys, out_settings) + CALL generate_out_write_acc(ovid%WB_hys, 'WB_hys', out%WB_hys, REAL(ssnow%wb_hys, 4),& + (/-1.0e36,1.0e36/), patchout%WB_hys, out_settings) + CALL generate_out_write_acc(ovid%SSAT_hys, 'SSAT_hys', out%SSAT_hys, REAL(ssnow%ssat_hys, 4),& + (/-1.0e36,1.0e36/), patchout%SSAT_hys, out_settings) + CALL generate_out_write_acc(ovid%WATR_hys, 'WATR_hys', out%WATR_hys, REAL(ssnow%watr_hys, 4),& + (/-1.0e36,1.0e36/), patchout%WATR_hys, out_settings) + CALL generate_out_write_acc(ovid%hys_fac, 'hys_fac', out%hys_fac, REAL(ssnow%hys_fac, 4),& + (/-1.0e36,1.0e36/), patchout%hys_fac, out_settings) + END IF +! end of block - rk4417 - phase2 + IF (output%Qrecharge) THEN ! recharge rate CALL generate_out_write_acc(ovid%Qrecharge, 'Qrecharge', out%Qrecharge, REAL(ssnow%Qrecharge, 4), ranges%Qrecharge, patchout%Qrecharge, out_settings) @@ -2680,8 +2708,9 @@ SUBROUTINE create_restart(logn, dels, ktau, soil, veg, ssnow, canstoID, albsoilsnID, gammzzID, tggsnID, sghfluxID, & ghfluxID, runoffID, rnof1ID, rnof2ID, gaID, dgdtgID, & fevID, fesID, fhsID, wbtot0ID, osnowd0ID, cplantID, & - csoilID, tradID, albedoID, gwID + csoilID, tradID, albedoID, gwID, subdzID ! added subdzID - rk4417 - phase2 INTEGER :: h0ID, snowliqID, SID, TsurfaceID, scondsID, nsnowID, TsoilID + INTEGER :: hys(6) ! inserted by rk4417 - phase2 CHARACTER(LEN=10) :: todaydate, nowtime ! used to timestamp netcdf file ! CHARACTER :: FRST_OUT*100, CYEAR*4 CHARACTER :: FRST_OUT*200, CYEAR*4 @@ -2964,6 +2993,31 @@ SUBROUTINE create_restart(logn, dels, ktau, soil, veg, ssnow, !$ CALL define_ovar(ncid_restart, rpid%rs20, 'rs20', '-', & !$ 'Soil respiration coefficient at 20C', & !$ .TRUE., 'real', 0, 0, 0, mpID, dummy, .TRUE.) + +! inserted block below - rk4417 - phase2 + CALL define_ovar(ncid_restart, rpid%froot, 'froot', '-', & + 'Fraction of roots in each soil layer', & + .TRUE., soilID, 'soil', 0, 0, 0, mpID, dummy, .TRUE.) + CALL define_ovar(ncid_restart, rpid%bch, 'bch', '-', & + 'Parameter b, Campbell eqn 1985', & + .TRUE., 'real', 0, 0, 0, mpID, dummy, .TRUE.) + CALL define_ovar(ncid_restart, rpid%hyds, 'hyds', 'mm/s', & ! MMY m/s->mm/s & + 'Hydraulic conductivity @ saturation', & + .TRUE., 'real', 0, 0, 0, mpID, dummy, .TRUE.) + CALL define_ovar(ncid_restart, rpid%sucs, 'sucs', 'mm', & ! MMY m->mm & + 'Suction @ saturation', .TRUE., & + 'real', 0, 0, 0, mpID, dummy, .TRUE.) + CALL define_ovar(ncid_restart, rpid%css, 'css', 'J/kg/C', & + 'Heat capacity of soil minerals', & + .TRUE., 'real', 0, 0, 0, mpID, dummy, .TRUE.) + CALL define_ovar(ncid_restart, rpid%rhosoil, 'rhosoil', 'kg/m^3', & + 'Density of soil minerals', & + .TRUE., 'real', 0, 0, 0, mpID, dummy, .TRUE.) + CALL define_ovar(ncid_restart, rpid%rs20, 'rs20', '-', & + 'Soil respiration coefficient at 20C', & + .TRUE., 'real', 0, 0, 0, mpID, dummy, .TRUE.) +! end of block - rk4417 - phase2 + CALL define_ovar(ncid_restart, rpid%albsoil, 'albsoil', '-', & 'Soil reflectance', .TRUE., & radID, 'radiation', 0, 0, 0, mpID, dummy, .TRUE.) @@ -3097,6 +3151,27 @@ SUBROUTINE create_restart(logn, dels, ktau, soil, veg, ssnow, .TRUE.,'real',0,0,0,mpID,dummy,.TRUE.) END IF ! SLI soil model + +! inserted block below - rk4417 - phase2 + if (cable_user%gw_model) then + CALL define_ovar(ncid_restart,hys(1),'wb_hys','-',& + 'water (volumetric) at dry/wet switch', & + .TRUE.,soilID,'soil',0,0,0,mpID,dummy,.TRUE.) + CALL define_ovar(ncid_restart,hys(2),'smp_hys','-',& + 'smp [mm] at dry/wet switch', & + .TRUE.,soilID,'soil',0,0,0,mpID,dummy,.TRUE.) + CALL define_ovar(ncid_restart,hys(3),'ssat_hys','-',& + 'ssat water (volumetric) from hyst', & + .TRUE.,soilID,'soil',0,0,0,mpID,dummy,.TRUE.) + CALL define_ovar(ncid_restart,hys(4),'watr_hys','-',& + 'ssat water (volumetric) from hyst', & + .TRUE.,soilID,'soil',0,0,0,mpID,dummy,.TRUE.) + CALL define_ovar(ncid_restart,hys(5),'hys_fac','-',& + 'water (volumetric) at dry/wet switch', & + .TRUE.,soilID,'soil',0,0,0,mpID,dummy,.TRUE.) + end if +! end of block - rk4417 - phase2 + ! Write global attributes for file: CALL DATE_AND_TIME(todaydate, nowtime) todaydate = todaydate(1:4)//'/'//todaydate(5:6)//'/'//todaydate(7:8) @@ -3385,6 +3460,21 @@ SUBROUTINE create_restart(logn, dels, ktau, soil, veg, ssnow, ranges%default_s, patchout_var, out_settings) END IF +! inserted IF block below - rk4417 - phase2 + IF (cable_user%gw_model) THEN + out_settings%dimswitch = "soil" + CALL check_and_write(hys(1),'wb_hys',REAL(ssnow%wb_hys,4), & + (/0.0,1.0/),patchout_var,out_settings) + CALL check_and_write(hys(2),'smp_hys',REAL(ssnow%smp_hys,4), & + (/-1.0e10,1.0e10/),patchout_var,out_settings) + CALL check_and_write(hys(3),'ssat_hys',REAL(ssnow%ssat_hys,4), & + (/0.0,1.0/),patchout_var,out_settings) + CALL check_and_write(hys(4),'watr_hys',REAL(ssnow%watr_hys,4), & + (/0.0,1.0/),patchout_var,out_settings) + CALL check_and_write(hys(5),'hys_fac',REAL(ssnow%hys_fac,4), & + (/0.0,1.0/),patchout_var,out_settings) + END IF + ! Close restart file ok = NF90_CLOSE(ncid_restart) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 44ebc6b9d..31829f011 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -59,7 +59,7 @@ MODULE cable_param_module USE phenvariable USE cable_abort_module USE cable_IO_vars_module - USE cable_common_module, ONLY: cable_user, gw_params + USE cable_common_module, ONLY: cable_user, gw_params, filename ! added 'filename' for now - rk4417 - phase2 USE cable_pft_params_mod USE cable_soil_params_mod USE CABLE_LUC_EXPT, ONLY: LUC_EXPT, LUC_EXPT_TYPE, LUC_EXPT_SET_TILES @@ -67,7 +67,7 @@ MODULE cable_param_module PRIVATE PUBLIC get_default_params, write_default_params, derived_parameters, & check_parameter_values, report_parameters, parID_type, & - write_cnp_params, consistency_ice_veg_soil + write_cnp_params, consistency_ice_veg_soil, GWspatialParameters ! GWspatialParameters added by rk4417 - phase2 INTEGER :: patches_in_parfile=4 ! # patches in default global parameter ! file @@ -105,22 +105,22 @@ MODULE cable_param_module REAL, DIMENSION(:, :), ALLOCATABLE :: insilt REAL, DIMENSION(:, :), ALLOCATABLE :: insand - !MD temp vars for reading in aquifer properties - LOGICAL :: found_explicit_gw_parameters - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWbch - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWssat - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWhyds - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWsucs - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWrhosoil - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWclay - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWsilt - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWsand - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWWatr - REAL, DIMENSION(:, :), ALLOCATABLE :: inWatr - REAL, DIMENSION(:, :), ALLOCATABLE :: inSlope - REAL, DIMENSION(:, :), ALLOCATABLE :: inGWdz - REAL, DIMENSION(:, :), ALLOCATABLE :: inSlopeSTD - REAL, DIMENSION(:, :), ALLOCATABLE :: inORG +! !MD temp vars for reading in aquifer properties ! block commented out by rk4417 - phase2 +! LOGICAL :: found_explicit_gw_parameters +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWbch +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWssat +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWhyds +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWsucs +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWrhosoil +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWclay +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWsilt +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWsand +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWWatr +! REAL, DIMENSION(:, :), ALLOCATABLE :: inWatr +! REAL, DIMENSION(:, :), ALLOCATABLE :: inSlope +! REAL, DIMENSION(:, :), ALLOCATABLE :: inGWdz +! REAL, DIMENSION(:, :), ALLOCATABLE :: inSlopeSTD +! REAL, DIMENSION(:, :), ALLOCATABLE :: inORG ! vars intro for Ticket #27 INTEGER, DIMENSION(:, :), ALLOCATABLE :: inSoilColor @@ -156,6 +156,12 @@ SUBROUTINE get_default_params(logn, vegparmnew, LUC_EXPT) INTEGER :: nlon INTEGER :: nlat +! Two calls below inserted by rk4417 - phase2 + ! Get parameter values for all default veg and soil types: + !CALL get_type_parameters(logn, vegparmnew, classification) + CALL cable_pft_params() + CALL cable_soil_params() + WRITE(logn,*) ' Reading grid info from ', TRIM(filename%type) WRITE(logn,*) ' And assigning C4 fraction according to veg classification.' WRITE(logn,*) @@ -178,16 +184,24 @@ SUBROUTINE get_default_params(logn, vegparmnew, LUC_EXPT) CALL spatialSoil(nlon, nlat, logn) ENDIF - ! Get parameter values for all default veg and soil types: - !CALL get_type_parameters(logn, vegparmnew, classification) - CALL cable_pft_params() - CALL cable_soil_params() +! Block below moved to the top - rk4417 - phase2 +! ! Get parameter values for all default veg and soil types: +! !CALL get_type_parameters(logn, vegparmnew, classification) +! CALL cable_pft_params() +! CALL cable_soil_params() + ! include prescribed soil colour in determining albedo - Ticket #27 IF (calcsoilalbedo) THEN CALL read_soilcolor(logn) END IF +! line below inserted by rk4417 - phase2 + IF (cable_user%force_npatches_as .GT. 0) npatch=cable_user%force_npatches_as +! MMY@13April force_npatches_as seems a new flag read from cable.nml +! the default value is -1 set in TYPE kbl_user_switches in cable_common +! TYPE kbl_user_switches was moved from cable_common_module to cable_runtime_opts_mod - rk4417 - phase2 + ! count to obtain 'landpt', 'max_vegpatches' and 'mp' CALL countPatch(nlon, nlat, npatch) @@ -249,7 +263,10 @@ SUBROUTINE read_gridinfo(nlon, nlat, npatch) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error inquiring y dimension.') ok = NF90_INQUIRE_DIMENSION(ncid, yID, LEN=nlat) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error getting y dimension.') - IF(.NOT. exists%patch) THEN + IF(.NOT. exists%patch) THEN ! MMY@13April here is right, when npatch doesn't exist in met file (i.e. exists%patch=False), + ! read npatch from gridinfo file. Note that it is to say, + ! npatch in meteorology file has a priority over in gridinfo + ok = NF90_INQ_DIMID(ncid, 'patch', pID) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error inquiring patch dimension.') ok = NF90_INQUIRE_DIMENSION(ncid, pID, LEN=npatch) @@ -473,14 +490,14 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) INTEGER, INTENT(IN) :: logn ! log file unit number ! local variables - INTEGER :: ncid, ok, ii, jj, kk, ok2, ncid_elev + INTEGER :: ncid, ok, ii, jj, kk, ok2 !, ncid_elev ! is not used - rk4417 - phase2 INTEGER :: xID, yID, fieldID INTEGER :: xlon, xlat - REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: indummy +! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: indummy ! appears only in commented code - rk4417 - phase2 REAL, DIMENSION(:,:), ALLOCATABLE :: sfact, dummy2 REAL, DIMENSION(:,:), ALLOCATABLE :: in2alb - ok = NF90_OPEN(filename%type, 0, ncid) +! ok = NF90_OPEN(filename%type, 0, ncid) ! moved below - rk4417 - phase2 ALLOCATE( in2alb(nlon, nlat) ) ! local ALLOCATE( dummy2(nlon, nlat) ) ! local @@ -498,16 +515,9 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ALLOCATE( insilt(nlon, nlat) ) ALLOCATE( insand(nlon, nlat) ) - !MD Aquifer properties - ALLOCATE( inGWssat(nlon, nlat) ) - ALLOCATE( inGWbch(nlon, nlat) ) - ALLOCATE( inGWhyds(nlon, nlat) ) - ALLOCATE( inGWsucs(nlon, nlat) ) - ALLOCATE( inGWrhosoil(nlon, nlat) ) - ALLOCATE( inGWWatr(nlon, nlat) ) - ALLOCATE( inWatr(nlon, nlat) ) - ALLOCATE( inORG(nlon, nlat) ) - + ! MMY start reading from gridinfo file ! 2 lines inserted by rk4417 - phase2 + ok = NF90_OPEN(filename%type, 0, ncid) + IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error MMY finding gridinfo file.') ! MMY ! 1 ok = NF90_INQ_VARID(ncid, 'swilt', fieldID) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error finding variable swilt.') @@ -574,90 +584,17 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ok = NF90_GET_VAR(ncid, fieldID, in2alb) IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error reading variable UM albedo') - !MD try to read aquifer properties from the file - ! if they don't exist set aquifer properties to the same as the soil - ok = NF90_INQ_VARID(ncid, 'Watr', fieldID) - WRITE(*,*) NF90_NOERR - ok2= ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inWatr) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inWatr(:,:) = 0.05 - END IF - - found_explicit_gw_parameters = .TRUE. - - ok = NF90_INQ_VARID(ncid, 'GWssat', fieldID) - WRITE(*,*) NF90_NOERR - ok2= ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inGWssat) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inGWssat(:,:) = inssat(:,:) - found_explicit_gw_parameters = .FALSE. - END IF - - ok = NF90_INQ_VARID(ncid, 'GWWatr', fieldID) - ok2 = ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inGWssat) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inGWWatr(:,:) = 0.05 - END IF - - ok = NF90_INQ_VARID(ncid, 'GWsucs', fieldID) - ok2 = ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inGWsucs) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inGWsucs(:,:) = ABS(insucs(:,:)) * 1000.0 - found_explicit_gw_parameters = .FALSE. - END IF - - ok = NF90_INQ_VARID(ncid, 'GWbch', fieldID) - ok2 = ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inGWbch) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inGWbch(:,:) = inbch(:,:) - found_explicit_gw_parameters = .FALSE. - END IF - - ok = NF90_INQ_VARID(ncid, 'GWhyds', fieldID) - ok2 = ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inGWhyds) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inGWhyds(:,:) = inhyds(:,:)*1000.0 - found_explicit_gw_parameters = .FALSE. - END IF - - ok = NF90_INQ_VARID(ncid, 'GWrhosoil', fieldID) - ok2 = ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inGWrhosoil) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inGWrhosoil(:,:) = inrhosoil(:,:) - END IF - - ok = NF90_INQ_VARID(ncid, 'organic', fieldID) - ok2 = ok - IF (ok .EQ. NF90_NOERR) THEN - ok2 = NF90_GET_VAR(ncid, fieldID, inORG) - WRITE(logn,*) 'READ FORG FROM THE DATA FILE, yeidling ' - WRITE(logn,*) 'A maximum value of ',MAXVAL(inORG),' and min val of',MINVAL(inORG) - END IF - IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN - inORG(:,:) = 0.0 - WRITE(logn,*) 'COULD NOT READ FORG FROM THR SRF FILE setting to 0.0' - END IF + ! ok = NF90_INQ_VARID(ncid, 'organic', fieldID) + ! ok2 = ok + ! IF (ok .EQ. NF90_NOERR) THEN + ! ok2 = NF90_GET_VAR(ncid, fieldID, inORG) + ! WRITE(logn,*) 'READ FORG FROM THE DATA FILE, yeidling ' + ! WRITE(logn,*) 'A maximum value of ',MAXVAL(inORG),' and min val of',MINVAL(inORG) + ! END IF + ! IF ((ok2 .NE. NF90_NOERR) .OR. (ok .NE. NF90_NOERR)) THEN + ! inORG(:,:) = 0.0 + ! WRITE(logn,*) 'COULD NOT READ FORG FROM THR SRF FILE setting to 0.0' + ! END IF ! Use this code if need to process original UM file soil fields into CABLE @@ -734,7 +671,7 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) ! CALL NSflip(nlon,nlat,in2alb) ok = NF90_CLOSE(ncid) - IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error closing IGBP soil map.') + IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error closing gridinfo file.') ! Code if using UM soil file ! unit change and glacial-point check were done in preprocessing @@ -775,49 +712,49 @@ SUBROUTINE spatialSoil(nlon, nlat, logn) inALB(:, :, 1, 1) = sfact(:, :) * dummy2(:, :) - ALLOCATE(inSlope(nlon,nlat),stat=ok) - IF (ok .NE. 0) CALL nc_abort(ok, 'Error allocating inSlope ') - inSlope(:,:) = 0.0 + ! ALLOCATE(inSlope(nlon,nlat),stat=ok) + ! IF (ok .NE. 0) CALL nc_abort(ok, 'Error allocating inSlope ') + ! inSlope(:,:) = 0.0 - ALLOCATE(inSlopeSTD(nlon,nlat),stat=ok) - IF (ok .NE. 0) CALL nc_abort(ok, 'Error allocating inSlopeSTD ') - inSlopeSTD(:,:) = 0.0 + ! ALLOCATE(inSlopeSTD(nlon,nlat),stat=ok) + ! IF (ok .NE. 0) CALL nc_abort(ok, 'Error allocating inSlopeSTD ') + ! inSlopeSTD(:,:) = 0.0 - ALLOCATE(inGWdz(nlon,nlat),stat=ok) - IF (ok .NE. 0) CALL nc_abort(ok, 'Error allocating inGWdz ') - inGWdz(:,:) = 20.0 + ! ALLOCATE(inGWdz(nlon,nlat),stat=ok) + ! IF (ok .NE. 0) CALL nc_abort(ok, 'Error allocating inGWdz ') + ! inGWdz(:,:) = 20.0 - IF (cable_user%GW_MODEL) THEN - ok = NF90_OPEN(TRIM(filename%gw_elev),NF90_NOWRITE,ncid_elev) - IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error opening GW elev param file.') + ! IF (cable_user%GW_MODEL) THEN + ! ok = NF90_OPEN(TRIM(filename%gw_elev),NF90_NOWRITE,ncid_elev) + ! IF (ok /= NF90_NOERR) CALL nc_abort(ok, 'Error opening GW elev param file.') - ok = NF90_INQ_VARID(ncid_elev, 'slope', fieldID) - IF (ok /= NF90_NOERR) WRITE(logn,*) 'Error finding variable slope' - ok = NF90_GET_VAR(ncid_elev, fieldID, inSlope) - IF (ok /= NF90_NOERR) THEN - inSlope = 0.0 - WRITE(logn, *) 'Could not read slope data for SSGW, set to 0.0' - END IF + ! ok = NF90_INQ_VARID(ncid_elev, 'slope', fieldID) + ! IF (ok /= NF90_NOERR) WRITE(logn,*) 'Error finding variable slope' + ! ok = NF90_GET_VAR(ncid_elev, fieldID, inSlope) + ! IF (ok /= NF90_NOERR) THEN + ! inSlope = 0.0 + ! WRITE(logn, *) 'Could not read slope data for SSGW, set to 0.0' + ! END IF - ok = NF90_INQ_VARID(ncid_elev, 'slope_std', fieldID) !slope_std - IF (ok /= NF90_NOERR) WRITE(logn,*) 'Error finding variable slope std' - ok = NF90_GET_VAR(ncid_elev, fieldID, inSlopeSTD) - IF (ok /= NF90_NOERR) THEN - inSlopeSTD = 0.0 - WRITE(logn, *) 'Could not read slope stddev data for SSGW, set to 0.0' - END IF + ! ok = NF90_INQ_VARID(ncid_elev, 'slope_std', fieldID) !slope_std + ! IF (ok /= NF90_NOERR) WRITE(logn,*) 'Error finding variable slope std' + ! ok = NF90_GET_VAR(ncid_elev, fieldID, inSlopeSTD) + ! IF (ok /= NF90_NOERR) THEN + ! inSlopeSTD = 0.0 + ! WRITE(logn, *) 'Could not read slope stddev data for SSGW, set to 0.0' + ! END IF - ok = NF90_INQ_VARID(ncid_elev, 'dtb', fieldID) - IF (ok /= NF90_NOERR) WRITE(logn,*) 'Error finding variable dtb' - ok = NF90_GET_VAR(ncid_elev, fieldID, inGWdz) - IF (ok /= NF90_NOERR) THEN - inGWdz = 20.0 - WRITE(logn, *) 'Could not read dtb data for SSGW, set to 0.0' - END IF + ! ok = NF90_INQ_VARID(ncid_elev, 'dtb', fieldID) + ! IF (ok /= NF90_NOERR) WRITE(logn,*) 'Error finding variable dtb' + ! ok = NF90_GET_VAR(ncid_elev, fieldID, inGWdz) + ! IF (ok /= NF90_NOERR) THEN + ! inGWdz = 20.0 + ! WRITE(logn, *) 'Could not read dtb data for SSGW, set to 0.0' + ! END IF - ok = NF90_CLOSE(ncid_elev) + ! ok = NF90_CLOSE(ncid_elev) - ENDIF !running gw model + ! ENDIF !running gw model DEALLOCATE(in2alb, sfact, dummy2) @@ -1064,7 +1001,9 @@ SUBROUTINE countPatch(nlon, nlat, npatch) STOP END IF ! CLN added for npatches - ELSE IF ( npatch .GT. 1 ) THEN +! ELSE IF ( npatch .GT. 1 ) THEN +! line above replaced by below - rk4417 - phase2 + ELSE IF ( npatch .GT. 1 .AND. cable_user%force_npatches_as .LE. 0 ) THEN landpt(kk)%nap = 0 DO tt = 1, npatch @@ -1080,6 +1019,20 @@ SUBROUTINE countPatch(nlon, nlat, npatch) PRINT *, 'vegtype_metfile = ', vegtype_metfile(kk,:) STOP END IF +! inserted ELSE IF block below - rk4417 - phase2 + ! CLN added for npatches + ELSE IF ( npatch .GT. 1 .AND. cable_user%force_npatches_as .EQ. npatch) THEN + !force the number of tiles per cell + !will read in the tile info from the GWspatial subroutine + landpt(kk)%nap = npatch + ncount = ncount + landpt(kk)%nap + landpt(kk)%cend = ncount + IF (landpt(kk)%cend < landpt(kk)%cstart) THEN + PRINT *, 'Land point ', kk, ' does not have veg type!' + PRINT *, 'landpt%cstart, cend = ', landpt(kk)%cstart, landpt(kk)%cend + PRINT *, 'vegtype_metfile = ', vegtype_metfile(kk,:) + STOP + END IF ELSE ! assume nmetpatches to be 1 IF (nmetpatches == 1) THEN @@ -1130,7 +1083,8 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & vegparmnew, month, TFRZ, LUC_EXPT) ! Initialize many canopy_type, soil_snow_type, soil_parameter_type and ! roughness_type variables; - ! Calculate 'froot' from 'rootbeta' parameter; + ! Calculate 'froot' from 'rootbeta' parameter. Does not assign values to + ! the soil hydraulic parameters. ! Assign values from input file to their proper variables in soil_snow_type, ! soil_parameter_type, veg_parameter_type and patch_type; ! Prescribe parameters for each point based on its veg/soil type. @@ -1151,6 +1105,12 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & USE cable_common_module, ONLY : calcsoilalbedo,cable_user + USE cable_pft_params_mod, ONLY : vegin !added by rk4417 - phase2 + USE cable_soil_params_mod, ONLY : soilin !added by rk4417 - phase2 + +! Note that subroutine init_veg_from_vegin has been moved from cable_common.F90 +! to cable_parameters.F90 (this file) and no longer needs to be imported - rk4417 - phase2 + IMPLICIT NONE INTEGER, INTENT(IN) :: logn ! log file unit number INTEGER, INTENT(IN) :: month ! month of year @@ -1214,7 +1174,6 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & canopy%fe = 0.0 ! sensible heat flux !mrd ssnow%qrecharge = 0.0 - ssnow%GWwb = -1.0 ssnow%wtd = 1.0 canopy%sublayer_dz = 0.001 !could go into restart to ensure starting/stopping runs gives identical results !however the impact is negligible @@ -1247,7 +1206,9 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & END SELECT - +! line below added by rk4417 - phase2 + soil%zse_vec = REAL(SPREAD(soil%zse,1,mp),r_2) ! MMY@13April, this line is needed since zec_vec is defined in UM/init/cable_um_init_subrs.F90, + ! but it isn't used in offline CABLE !ELSE ! ! parameters that are not spatially dependent @@ -1335,6 +1296,7 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & ! Set initial snow depth and snow-free soil albedo + ! MMY Note @Oct2022 Please notice ssnow%albsoilsn is set here rather than in GWspatialParameters DO is = 1, landpt(e)%cend - landpt(e)%cstart + 1 ! each patch DO ir = 1, nrb IF (CABLE_USER%POPLUC) THEN !vh! use same soilalbedo for all land-use tiles @@ -1347,6 +1309,7 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & ENDIF END DO ! total depth, change from m to mm !see Ticket #57 + ! Ticket #57: changing the scaling factor from 1000 (water density) to 140 (old ice density) ! MMY @Oct2022 ssnow%snowd(landpt(e)%cstart + is - 1) & = inSND(landpt(e)%ilon, landpt(e)%ilat, is, month) * 140.0 END DO @@ -1357,101 +1320,104 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & inLAI(landpt(e)%ilon,landpt(e)%ilat,is) END DO - ! Set IGBP soil texture values, Q.Zhang @ 12/20/2010. - IF (soilparmnew) THEN - - soil%swilt(landpt(e)%cstart:landpt(e)%cend) = & - inswilt(landpt(e)%ilon, landpt(e)%ilat) - soil%sfc(landpt(e)%cstart:landpt(e)%cend) = & - insfc(landpt(e)%ilon, landpt(e)%ilat) - soil%ssat(landpt(e)%cstart:landpt(e)%cend) = & - inssat(landpt(e)%ilon, landpt(e)%ilat) - soil%bch(landpt(e)%cstart:landpt(e)%cend) = & - inbch(landpt(e)%ilon, landpt(e)%ilat) - soil%hyds(landpt(e)%cstart:landpt(e)%cend) = & - inhyds(landpt(e)%ilon, landpt(e)%ilat) - soil%sucs(landpt(e)%cstart:landpt(e)%cend) = & - -1.* ABS(insucs(landpt(e)%ilon, landpt(e)%ilat)) !ensure negative - soil%rhosoil(landpt(e)%cstart:landpt(e)%cend) = & - inrhosoil(landpt(e)%ilon, landpt(e)%ilat) - soil%css(landpt(e)%cstart:landpt(e)%cend) = & - incss(landpt(e)%ilon, landpt(e)%ilat) - soil%cnsd(landpt(e)%cstart:landpt(e)%cend) = & - incnsd(landpt(e)%ilon, landpt(e)%ilat) - - !possibly heterogeneous soil properties - ! Set all heterogeneous soil properties to their equivalent uniform values. - ! These values can be overridden later on by input values in files. - ! `smpc_vec` and `wbc_vec` are only used in the ground water and have no - ! equivalent in the uniform parameters (non _vec). So we do not initialise - ! these variables when we are not using the ground water scheme. - DO klev=1,ms - - soil%clay_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(inclay(landpt(e)%ilon, landpt(e)%ilat),r_2) - - soil%sand_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(insand(landpt(e)%ilon, landpt(e)%ilat),r_2) - - soil%silt_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(insilt(landpt(e)%ilon, landpt(e)%ilat),r_2) - - soil%rhosoil_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(inrhosoil(landpt(e)%ilon, landpt(e)%ilat),r_2) - - soil%org_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(inORG(landpt(e)%ilon, landpt(e)%ilat),r_2) - - soil%watr(landpt(e)%cstart:landpt(e)%cend,klev) = & - REAL(inWatr(landpt(e)%ilon, landpt(e)%ilat),r_2) - - soil%css_vec(landpt(e)%cstart:landpt(e)%cend, klev) = & - REAL(incss(landpt(e)%ilon, landpt(e)%ilat), r_2) - - soil%cnsd_vec(landpt(e)%cstart:landpt(e)%cend, klev) = & - REAL(incnsd(landpt(e)%ilon, landpt(e)%ilat), r_2) + ! ! Set IGBP soil texture values, Q.Zhang @ 12/20/2010. + ! IF (soilparmnew) THEN - END DO + ! soil%swilt(landpt(e)%cstart:landpt(e)%cend) = & + ! inswilt(landpt(e)%ilon, landpt(e)%ilat) + ! soil%sfc(landpt(e)%cstart:landpt(e)%cend) = & + ! insfc(landpt(e)%ilon, landpt(e)%ilat) + ! soil%ssat(landpt(e)%cstart:landpt(e)%cend) = & + ! inssat(landpt(e)%ilon, landpt(e)%ilat) + ! soil%bch(landpt(e)%cstart:landpt(e)%cend) = & + ! inbch(landpt(e)%ilon, landpt(e)%ilat) + ! soil%hyds(landpt(e)%cstart:landpt(e)%cend) = & + ! inhyds(landpt(e)%ilon, landpt(e)%ilat) + ! soil%sucs(landpt(e)%cstart:landpt(e)%cend) = & + ! -1.* ABS(insucs(landpt(e)%ilon, landpt(e)%ilat)) !ensure negative + ! soil%rhosoil(landpt(e)%cstart:landpt(e)%cend) = & + ! inrhosoil(landpt(e)%ilon, landpt(e)%ilat) + ! soil%css(landpt(e)%cstart:landpt(e)%cend) = & + ! incss(landpt(e)%ilon, landpt(e)%ilat) + ! soil%cnsd(landpt(e)%cstart:landpt(e)%cend) = & + ! incnsd(landpt(e)%ilon, landpt(e)%ilat) - !Aquifer properties same as bottom soil layer for now - soil%GWsucs_vec(landpt(e)%cstart:landpt(e)%cend) = & - REAL(inGWsucs(landpt(e)%ilon, landpt(e)%ilat),r_2) + ! !possibly heterogeneous soil properties + ! ! Set all heterogeneous soil properties to their equivalent uniform values. + ! ! These values can be overridden later on by input values in files. + ! ! `smpc_vec` and `wbc_vec` are only used in the ground water and have no + ! ! equivalent in the uniform parameters (non _vec). So we do not initialise + ! ! these variables when we are not using the ground water scheme. + ! DO klev=1,ms - soil%GWhyds_vec(landpt(e)%cstart:landpt(e)%cend) = & - REAL(inGWhyds(landpt(e)%ilon, landpt(e)%ilat),r_2) + ! soil%clay_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(inclay(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%GWbch_vec(landpt(e)%cstart:landpt(e)%cend) = & - REAL(inGWbch(landpt(e)%ilon, landpt(e)%ilat),r_2) + ! soil%sand_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(insand(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%GWrhosoil_vec(landpt(e)%cstart:landpt(e)%cend) = & - REAL(inGWrhosoil(landpt(e)%ilon, landpt(e)%ilat),r_2) + ! soil%silt_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(insilt(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%GWssat_vec(landpt(e)%cstart:landpt(e)%cend) = & - REAL(inGWssat(landpt(e)%ilon, landpt(e)%ilat),r_2) + ! soil%rhosoil_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(inrhosoil(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%GWwatr(landpt(e)%cstart:landpt(e)%cend) = & - soil%watr(landpt(e)%cstart:landpt(e)%cend,ms) + ! soil%org_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(inORG(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%slope(landpt(e)%cstart:landpt(e)%cend) = & - MIN(MAX(inSlope(landpt(e)%ilon,landpt(e)%ilat),1e-8),0.95) + ! soil%watr(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(inWatr(landpt(e)%ilon, landpt(e)%ilat),r_2) - soil%slope_std(landpt(e)%cstart:landpt(e)%cend) = & - MIN(MAX(inSlopeSTD(landpt(e)%ilon,landpt(e)%ilat),1e-8),0.95) + ! soil%zse_vec(landpt(e)%cstart:landpt(e)%cend,klev) = & + ! REAL(soil%zse(landpt(e)%cstart:landpt(e)%cend), r_2) - soil%GWdz(landpt(e)%cstart:landpt(e)%cend) = & - inGWdz(landpt(e)%ilon,landpt(e)%ilat) + ! soil%css_vec(landpt(e)%cstart:landpt(e)%cend, klev) = & + ! REAL(incss(landpt(e)%ilon, landpt(e)%ilat), r_2) - ! vh ! - soil%silt(landpt(e)%cstart:landpt(e)%cend) = & - insilt(landpt(e)%ilon, landpt(e)%ilat) + ! soil%cnsd_vec(landpt(e)%cstart:landpt(e)%cend, klev) = & + ! REAL(incnsd(landpt(e)%ilon, landpt(e)%ilat), r_2) - soil%sand(landpt(e)%cstart:landpt(e)%cend) = & - insand(landpt(e)%ilon, landpt(e)%ilat) + ! END DO - soil%clay(landpt(e)%cstart:landpt(e)%cend) = & - inclay(landpt(e)%ilon, landpt(e)%ilat) + ! !Aquifer properties same as bottom soil layer for now + ! soil%GWsucs_vec(landpt(e)%cstart:landpt(e)%cend) = & + ! REAL(inGWsucs(landpt(e)%ilon, landpt(e)%ilat),r_2) - ENDIF + ! soil%GWhyds_vec(landpt(e)%cstart:landpt(e)%cend) = & + ! REAL(inGWhyds(landpt(e)%ilon, landpt(e)%ilat),r_2) + + ! soil%GWbch_vec(landpt(e)%cstart:landpt(e)%cend) = & + ! REAL(inGWbch(landpt(e)%ilon, landpt(e)%ilat),r_2) + + ! soil%GWrhosoil_vec(landpt(e)%cstart:landpt(e)%cend) = & + ! REAL(inGWrhosoil(landpt(e)%ilon, landpt(e)%ilat),r_2) + + ! soil%GWssat_vec(landpt(e)%cstart:landpt(e)%cend) = & + ! REAL(inGWssat(landpt(e)%ilon, landpt(e)%ilat),r_2) + + ! soil%GWwatr(landpt(e)%cstart:landpt(e)%cend) = & + ! soil%watr(landpt(e)%cstart:landpt(e)%cend,ms) + + ! soil%slope(landpt(e)%cstart:landpt(e)%cend) = & + ! MIN(MAX(inSlope(landpt(e)%ilon,landpt(e)%ilat),1e-8),0.95) + + ! soil%slope_std(landpt(e)%cstart:landpt(e)%cend) = & + ! MIN(MAX(inSlopeSTD(landpt(e)%ilon,landpt(e)%ilat),1e-8),0.95) + + ! soil%GWdz(landpt(e)%cstart:landpt(e)%cend) = & + ! inGWdz(landpt(e)%ilon,landpt(e)%ilat) + + ! ! vh ! + ! soil%silt(landpt(e)%cstart:landpt(e)%cend) = & + ! insilt(landpt(e)%ilon, landpt(e)%ilat) + + ! soil%sand(landpt(e)%cstart:landpt(e)%cend) = & + ! insand(landpt(e)%ilon, landpt(e)%ilat) + + ! soil%clay(landpt(e)%cstart:landpt(e)%cend) = & + ! inclay(landpt(e)%ilon, landpt(e)%ilat) + + ! ENDIF ! vars intro for Ticket #27 IF (calcsoilalbedo) THEN @@ -1580,18 +1546,20 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & ! vegin%gswmin, vegin%conkc0,vegin%conko0,vegin%ekc,vegin%eko ) ! ! vegf_temp,urbanf_temp,lakef_temp,icef_temp, & - IF (ALLOCATED(inGWsucs )) DEALLOCATE(inGWsucs) - IF (ALLOCATED(inGWhyds )) DEALLOCATE(inGWhyds) - IF (ALLOCATED(inGWbch )) DEALLOCATE(inGWbch) - IF (ALLOCATED(inGWsilt )) DEALLOCATE(inGWsilt) - IF (ALLOCATED(inGWsand )) DEALLOCATE(inGWsand) - IF (ALLOCATED(inGWclay )) DEALLOCATE(inGWclay) - IF (ALLOCATED(inGWssat )) DEALLOCATE(inGWssat) - IF (ALLOCATED(inGWWatr )) DEALLOCATE(inGWWatr) - IF (ALLOCATED(inWatr )) DEALLOCATE(inWatr) - IF (ALLOCATED(inSlope )) DEALLOCATE(inSlope) - IF (ALLOCATED(inSlopeSTD)) DEALLOCATE(inSlopeSTD) - IF (ALLOCATED(inORG )) DEALLOCATE(inORG) +! Block below commented out by rk4417 - phase2 +! IF (ALLOCATED(inGWsucs )) DEALLOCATE(inGWsucs) +! IF (ALLOCATED(inGWhyds )) DEALLOCATE(inGWhyds) +! IF (ALLOCATED(inGWbch )) DEALLOCATE(inGWbch) +! IF (ALLOCATED(inGWsilt )) DEALLOCATE(inGWsilt) +! IF (ALLOCATED(inGWsand )) DEALLOCATE(inGWsand) +! IF (ALLOCATED(inGWclay )) DEALLOCATE(inGWclay) +! IF (ALLOCATED(inGWssat )) DEALLOCATE(inGWssat) +! IF (ALLOCATED(inGWWatr )) DEALLOCATE(inGWWatr) +! IF (ALLOCATED(inWatr )) DEALLOCATE(inWatr) +! IF (ALLOCATED(inSlope )) DEALLOCATE(inSlope) +! IF (ALLOCATED(inSlopeSTD)) DEALLOCATE(inSlopeSTD) +! IF (ALLOCATED(inORG )) DEALLOCATE(inORG) + ! if using old format veg_parm input file, need to define veg%deciduous ! BP dec 2007 ! IF (.NOT. vegparmnew) THEN @@ -1646,6 +1614,13 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & ssnow%wbliq = ssnow%wb - ssnow%wbice ssnow%GWwb = 0.9*soil%ssat +! 5 lines below inserted by rk4417 - phase2 + ssnow%wb_hys = -1.0e+36 + ssnow%hys_fac = 1.0 + ssnow%watr_hys = soil%watr + ssnow%ssat_hys = soil%ssat_vec + ssnow%smp_hys = -1.0e+36 + !IF(hide%Ticket49Bug5) THEN ! vh_js ! neeed to remove this if to enable the code below @@ -1678,16 +1653,18 @@ SUBROUTINE write_default_params(met, air, ssnow, veg, bgc, & veg%disturbance_intensity = 0. ENDIF +! MMY note @Oct2022 Important differ to MMY's trunk: assume depth to bedrock = 4.6m + aquifer soil%GWdz = MAX(1.0,MIN(20.0,soil%GWdz - SUM(soil%zse,dim=1))) - !set vectorized versions as same as defaut for now - soil%swilt_vec(:,:) = REAL(SPREAD(soil%swilt(:),2,ms),r_2) - soil%sfc_vec(:,:) = REAL(SPREAD(soil%sfc(:),2,ms),r_2) - soil%sucs_vec(:,:) = REAL(SPREAD(soil%sucs(:),2,ms),r_2) - soil%bch_vec(:,:) = REAL(SPREAD(soil%bch(:),2,ms),r_2) - soil%ssat_vec(:,:) = REAL(SPREAD(soil%ssat(:),2,ms),r_2) - soil%hyds_vec(:,:) = REAL(SPREAD(soil%hyds(:),2,ms),r_2) - soil%zse_vec(:,:) = REAL(SPREAD(soil%zse(:),1,mp),r_2) +! Block below commented out by rk4417 - phase2 +! ____________ MMY @Oct2022 comment out ________ + !set vectorized versions as same as default for now +! soil%swilt_vec(:,:) = REAL(SPREAD(soil%swilt(:),2,ms),r_2) +! soil%sfc_vec(:,:) = REAL(SPREAD(soil%sfc(:),2,ms),r_2) +! soil%sucs_vec(:,:) = REAL(SPREAD(soil%sucs(:),2,ms),r_2) +! soil%bch_vec(:,:) = REAL(SPREAD(soil%bch(:),2,ms),r_2) +! soil%ssat_vec(:,:) = REAL(SPREAD(soil%ssat(:),2,ms),r_2) +! soil%hyds_vec(:,:) = REAL(SPREAD(soil%hyds(:),2,ms),r_2) END SUBROUTINE write_default_params !============================================================================= @@ -2158,6 +2135,11 @@ SUBROUTINE derived_parameters(soil, sum_flux, bal, ssnow, veg, rough) end do if ((gw_params%MaxSatFraction .lt. -9999.9) .and. (mp .eq. 1)) soil%slope(:) = 0.01 + ELSE !not gw model + + !These are not used when gw_model == false + soil%watr = 0._r_2 + soil%GWwatr = 0._r_2 END IF @@ -2366,10 +2348,30 @@ SUBROUTINE check_parameter_values(soil, veg, ssnow) IF(ANY(soil%isoilm(landpt(i)%cstart:(landpt(i)%cstart + landpt(i)%nap & - 1)) < 1 ) .OR. ANY(soil%isoilm(landpt(i)%cstart:(landpt(i)%cstart & + landpt(i)%nap - 1)) > mstype)) THEN - WRITE(*,*) 'SUBROUTINE load_parameters:' - WRITE(*,*) 'Land point number:',i + WRITE(*,*) 'SUBROUTINE load_parameters: soil < 1 or > ',mstype + DO j=landpt(i)%cstart,landpt(i)%cend + IF (soil%isoilm(j) .LT. 1 .OR. soil%isoilm(j) .GT. mstype) THEN + WRITE(*,*) 'AT ',i,j,' indices ' + WRITE(*,*) '~~~~~~~~~~~~~~lat=>',& + patch(landpt(i)%cstart:landpt(i)%cend)%latitude + WRITE(*,*) '~~~~~~~~~~~~~~lon=>',& + patch(landpt(i)%cstart:landpt(i)%cend)%longitude + WRITE(*,*) '~~~~~~~~~~~~~~soil%isoilm',soil%isoilm(j) + END IF + ENDDO CALL abort('Unknown soil type! Aborting.') END IF + +! IF block above replaces one below - rk4417 - phase2 + +! IF(ANY(soil%isoilm(landpt(i)%cstart:(landpt(i)%cstart + landpt(i)%nap & +! - 1)) < 1 ) .OR. ANY(soil%isoilm(landpt(i)%cstart:(landpt(i)%cstart & +! + landpt(i)%nap - 1)) > mstype)) THEN +! WRITE(*,*) 'SUBROUTINE load_parameters:' +! WRITE(*,*) 'Land point number:',i +! CALL abort('Unknown soil type! Aborting.') +! END IF + ! Check patch fractions sum to 1 in each grid cell: IF((SUM(patch(landpt(i)%cstart:landpt(i)%cend)%frac) - 1.0) & > 1.0E-6) THEN diff --git a/src/offline/cable_pft_params.F90 b/src/offline/cable_pft_params.F90 index 92d09d971..3303249c7 100644 --- a/src/offline/cable_pft_params.F90 +++ b/src/offline/cable_pft_params.F90 @@ -42,6 +42,10 @@ MODULE cable_pft_params_mod g1(ntype_max), & zr(ntype_max), & clitt(ntype_max), & +! froot1 - froot6 do not serve any purpose +! their values inside offline/pft_params.nml are not used +! veg%froot is initialized instead from formula inside offline/cable_parameters.F90 +! rk4417 - phase2 froot1(ntype_max), & froot2(ntype_max), & froot3(ntype_max), & @@ -60,10 +64,10 @@ MODULE cable_pft_params_mod ratecp3(ntype_max), & refl1(ntype_max), & refl2(ntype_max), & - refl3(ntype_max), & + refl3(ntype_max), & ! not used - rk4417 - phase2 taul1(ntype_max), & taul2(ntype_max), & - taul3(ntype_max), & + taul3(ntype_max), & ! not used - rk4417 - phase2 dleaf(ntype_max), & lai(ntype_max) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 366aa58de..97a007ed1 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -803,10 +803,20 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) YYYY.EQ. CABLE_USER%YearEnd ) THEN ! evaluate spinup + ! =================== MMY_phase2 + ! commented out this region in favor of the one below - rk4417 + ! ===================== IF( ANY( ABS(ssnow%wb-soilMtemp)>delsoilM).OR. & ANY( ABS(ssnow%tgg-soilTtemp)>delsoilT) .OR. & MAXVAL(ABS(ssnow%GWwb-GWtemp),dim=1) > delgwM) THEN + ! =================== MMY_phase2 uncomment ===================== +! IF( (ANY( ABS(ssnow%wb-soilMtemp)>delsoilM).OR. & +! ANY( ABS(ssnow%tgg-soilTtemp)>delsoilT) .or. & +! maxval(ABS(ssnow%GWwb-GWtemp),dim=1) > delgwM) .and. & +! ( (int(ktau_tot/kend) .lt. cable_user%max_spins) .and.& +! (cable_user%max_spins .gt. 0) ) ) THEN + ! No complete convergence yet PRINT *, 'ssnow%wb : ', ssnow%wb PRINT *, 'soilMtemp: ', soilMtemp diff --git a/src/offline/cbl_model_driver_offline.F90 b/src/offline/cbl_model_driver_offline.F90 index 763e68075..1604c4537 100644 --- a/src/offline/cbl_model_driver_offline.F90 +++ b/src/offline/cbl_model_driver_offline.F90 @@ -188,8 +188,23 @@ SUBROUTINE cbm( ktau,dels, air, bgc, canopy, met, !update the various biophysics state variables ssnow%owetfac = ssnow%wetfac + +IF(cable_user%SOIL_STRUC=='default') THEN -CALL soil_snow(dels, soil, ssnow, canopy, met, bal,veg) + IF (cable_user%gw_model) THEN + CALL soil_snow_gw(dels, soil, ssnow, canopy, met, bal,veg) + ELSE + CALL soil_snow(dels, soil, ssnow, canopy, met, bal,veg) + ENDIF + +ELSEIF (cable_user%SOIL_STRUC=='sli') THEN + + IF (cable_user%test_new_gw) & + CALL sli_hydrology(dels,ssnow,soil,veg,canopy) + + CALL sli_main(ktau,dels,veg,soil,ssnow,met,canopy,air,rad,0) + +ENDIF !#539 - move snow_aging now after soil_snow - uses this timestep snow amount CALL snow_aging(ssnow%snage,mp,dels,ssnow%snowd,ssnow%osnowd,ssnow%tggsn(:,1),& diff --git a/src/science/canopy/cable_canopy.F90 b/src/science/canopy/cable_canopy.F90 index c7f1f4cb5..153490554 100644 --- a/src/science/canopy/cable_canopy.F90 +++ b/src/science/canopy/cable_canopy.F90 @@ -556,7 +556,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima ELSE - write(6,*) "SLI is not an option right now" - commented out by rk4417 - phase2 + write(6,*) "SLI is not an option right now" ! SLI SEB to get canopy%fhs, canopy%fess, canopy%ga ! (Based on old Tsoil, new canopy%tv, new canopy%fns) !H!CALL sli_main(1,dels,veg,soil,ssnow,met,canopy,air,rad,1) @@ -638,7 +638,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima canopy%ga = canopy%fns-canopy%fhs-canopy%fes ! *ssnow%cls ELSE - write(6,*) "SLI is not an option right now" - commented out by rk4417 - phase2 + write(6,*) "SLI is not an option right now" ! SLI SEB to get canopy%fhs, canopy%fess, canopy%ga ! (Based on old Tsoil, new canopy%tv, new canopy%fns) From 67f74b3315152472596e1c99cced74091f4b22a7 Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 29 Jul 2025 11:50:15 +1000 Subject: [PATCH 28/28] Changes to compile the library for coupled models --- CMakeLists.txt | 57 +++++--- build.bash | 2 +- src/coupled/esm/LAI_canopy_height_cbl.F90 | 156 ++++++++++++++++++++++ src/coupled/esm/cable_surface_types.F90 | 4 +- 4 files changed, 197 insertions(+), 22 deletions(-) create mode 100644 src/coupled/esm/LAI_canopy_height_cbl.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index a186c1f27..a67431ce1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,11 +36,11 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") "$<$:${CABLE_GNU_Fortran_FLAGS_DEBUG}>" ) endif() - +set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/include") if(CABLE_LIBRARY) add_library( - cable_science + cable STATIC src/science/albedo/cbl_albedo.F90 src/science/albedo/cbl_snow_albedo.F90 @@ -70,8 +70,8 @@ if(CABLE_LIBRARY) src/science/casa-cnp/casa_rplant.F90 src/science/casa-cnp/casa_sumcflux.F90 src/science/casa-cnp/casa_variable.F90 - src/science/gw_hydro/cable_gw_hydro.F90 - src/science/gw_hydro/cable_psm.F90 + #src/science/gw_hydro/cable_gw_hydro.F90 + #src/science/gw_hydro/cable_psm.F90 src/science/landuse/landuse3.F90 src/science/landuse/landuse_constant.F90 src/science/misc/cable_air.F90 @@ -79,9 +79,9 @@ if(CABLE_LIBRARY) src/science/misc/cable_climate.F90 src/science/pop/pop_constants.F90 src/science/pop/pop_def.F90 - src/science/pop/POP.F90 - src/science/pop/pop_io.F90 - src/science/pop/POPLUC.F90 + #src/science/pop/POP.F90 + #src/science/pop/pop_io.F90 + #src/science/pop/POPLUC.F90 src/science/pop/pop_types.F90 src/science/radiation/cbl_init_radiation.F90 src/science/radiation/cbl_radiation.F90 @@ -90,11 +90,11 @@ if(CABLE_LIBRARY) src/science/radiation/cbl_spitter.F90 src/science/roughness/cable_roughness.F90 src/science/roughness/roughnessHGT_effLAI_cbl.F90 - src/science/sli/cable_sli_main.F90 - src/science/sli/cable_sli_numbers.F90 - src/science/sli/cable_sli_roots.F90 - src/science/sli/cable_sli_solve.F90 - src/science/sli/cable_sli_utils.F90 + #src/science/sli/cable_sli_main.F90 + #src/science/sli/cable_sli_numbers.F90 + #src/science/sli/cable_sli_roots.F90 + #src/science/sli/cable_sli_solve.F90 + #src/science/sli/cable_sli_utils.F90 src/science/soilsnow/cbl_conductivity.F90 src/science/soilsnow/cbl_GW.F90 src/science/soilsnow/cbl_hyd_redistrib.F90 @@ -122,14 +122,33 @@ if(CABLE_LIBRARY) src/params/cable_maths_constants_mod.F90 src/util/cable_runtime_opts_mod.F90 src/util/cable_common.F90 - src/coupled/esm16/casa_offline_inout.F90 - src/coupled/esm16/casa_ncdf.F90 - src/coupled/esm16/cable_iovars.F90 - src/coupled/esm16/cable_surface_types.F90 - src/coupled/esm16/cable_define_types.F90 - src/coupled/esm16/cable_phenology.F90 - src/coupled/esm16/cable_LUC_EXPT.F90 + src/coupled/esm/cable_iovars.F90 + src/coupled/esm/cable_surface_types.F90 + src/coupled/esm/cable_define_types.F90 + src/coupled/esm/cable_soil_params_mod.F90 + src/coupled/esm/cable_pft_params_mod.F90 + src/coupled/esm/casa_landuse.F90 + src/coupled/esm/LAI_canopy_height_cbl.F90 + #src/shared/casa_offline_inout.F90 + #src/shared/casa_ncdf.F90 + #src/shared/cable_phenology.F90 + #src/shared/cable_LUC_EXPT.F90 + src/coupled/shared/cable_canopy_type_mod.F90 + src/coupled/shared/cable_soilsnow_type_mod.F90 + src/coupled/shared/cable_soil_type_mod.F90 + src/coupled/shared/cable_veg_type_mod.F90 + src/util/cable_climate_type_mod.F90 + src/util/masks_cbl.F90 ) + target_compile_definitions(cable PRIVATE UM_CBL) + target_include_directories(cable PUBLIC "$") + + install(TARGETS cable + EXPORT CABLEFortranTargets + LIBRARY DESTINATION lib + INCLUDES DESTINATION include) + + install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY}/ TYPE INCLUDE) else() diff --git a/build.bash b/build.bash index c698d596f..ee86a1b1e 100755 --- a/build.bash +++ b/build.bash @@ -62,7 +62,7 @@ while [ ${#} -gt 0 ]; do cmake_args+=(-DCABLE_MPI="ON") ;; -l|--library) - build_args+=(--target cable_science) + build_args+=(--target cable) cmake_args+=(-DCABLE_LIBRARY="ON") do_install=0 # Disable installation when only building the science library ;; diff --git a/src/coupled/esm/LAI_canopy_height_cbl.F90 b/src/coupled/esm/LAI_canopy_height_cbl.F90 new file mode 100644 index 000000000..864bf626e --- /dev/null +++ b/src/coupled/esm/LAI_canopy_height_cbl.F90 @@ -0,0 +1,156 @@ +!******************************COPYRIGHT******************************************** +! (c) CSIRO 2022. +! All rights reserved. +! +! This routine has been licensed to the other JULES partners for use and +! distribution under the JULES collaboration agreement, subject to the terms and +! conditions set out therein. +! +! [Met Office Ref SC0237] +!******************************COPYRIGHT******************************************** +MODULE cbl_LAI_canopy_height_mod + +!----------------------------------------------------------------------------- +! Description: +! Restricts the range of canopy height and LAI inherited from JULES/UM +! spatial maps +! +! This MODULE is USEd by: +! cable_land_albedo_mod_cbl.F90 +! +! This MODULE contains 1 public Subroutine: +! limit_HGT_LAI +! +! Code Owner: Please refer to ModuleLeaders.txt +! This file belongs in CABLE SCIENCE +! +! Code Description: +! Language: Fortran 90. +! This code is written to JULES coding standards v1. +!----------------------------------------------------------------------------- + +IMPLICIT NONE +PUBLIC :: limit_HGT_LAI +PRIVATE + +CONTAINS + +SUBROUTINE Limit_HGT_LAI( HGT_pft_temp, LAI_pft_cbl, HGT_pft_cbl, mp, land_pts,& + ntiles, npft, tile_pts, tile_index, tile_frac, & + L_tile_pts, LAI_pft, HGT_pft, CLAI_thresh ) + +!*## Purpose +! +! This SUBROUTINE checks that input values of leaf area and canopy height +! lie between set values (defined by PFT) and overwrites if necessary. +! These limits are needed to ensure that the canopy roughness properties of +! the land point are well behaved (see [[ruff_resist]]). +! +! The mapping from grid cells (land_pts,ntiles) to the (mp) structure +! is also undertaken. +! +! **This SUBROUTINE is active when CABLE is run within a coupled model +! (i.e. ACCESS)** +! +!## Method +! +! The cell-level input (land_pts,ntile) variables for leaf area `LAI_pft` and +! canopy height `HGT_pft` are checked for +! +! - whether the land_pt has a non-zero fraction of that tile (if zero +! fraction then LAI and canopy height are set to zero) +! - whether the LAI lies above a minimum value `CLAI_thresh` +! - whether the canopy height lies above a minimum value (which is PFT dependent) +! - whether non-vegetated tiles have a non-zero canopy height. +! +! Outputs are `LAI_pft_cbl` and `HGT_pft_cbl` +! +! **WARNINGS** +! +! - INTENT statements need to be added to the argument lists +! - hardwired indexing is used throughout. This SUBROUTINE assumes that +! 1. non-vegetated tiles take indexes 14 and onwards +! 2. tall vegetation occupies tile indexes 1-4 +! 3. other vegetation (shrubs/grasses/crops/wetlands) occupy tile indexes +! 5-13 +! - the limits on canopy height are hardwired +! 1. canopy height for tall vegetation is limited to be greater than 1m +! 2. canopy height for other vegetation is limited to be greater than 0.1m + +USE cable_surface_types_mod, ONLY: shrub_cable, evergreen_broadleaf, & + deciduous_broadleaf, aust_temperate, & + aust_tropical + +IMPLICIT NONE + +INTEGER, INTENT(IN) :: land_pts, ntiles, npft !! # land points here + !! max # of tiles, PFTs +!! CABLE 1-D vector of length "mp" active tiles +INTEGER, INTENT(IN) :: mp !! total number of tiles + +!! range limited LAI/ canopy height +REAL, INTENT(OUT) :: HGT_pft_temp(land_pts,ntiles) !! UM dimensions (m) + +!! PACKed to CABLE 1-D vector of length mp +REAL, INTENT(OUT) :: LAI_pft_cbl(mp) !! veg%vlai +REAL, INTENT(OUT) :: HGT_pft_cbl(mp) !! veg%hc + +!! From UM spatial maps via restart or ancillary. +REAL, INTENT(IN) :: tile_frac(land_pts,ntiles) !! frac. veg/non-veg type +INTEGER, INTENT(IN) :: tile_pts(ntiles) !! # tiles per PFT type +INTEGER, INTENT(IN) :: tile_index(land_pts,ntiles) !! land_pt index per tile +!! updates monthly +REAL, INTENT(IN) :: LAI_pft(land_pts, npft) !! LAI (m\(^2\)m\(^{-2}\)) +REAL, INTENT(IN) :: HGT_pft(land_pts, npft) !! canopy height (m) +!! logical mask. TRUE if tilefrac > 0 0 (OR some threshold) +LOGICAL, INTENT(IN) :: L_tile_pts(land_pts,ntiles) + +!! scalar constant threshold for "cell" to be considred vegetated +REAL, INTENT(IN) :: Clai_thresh ! minimum LAI threshold + +!local vars +INTEGER :: i,j, n +REAL :: LAI_pft_temp(land_pts,ntiles) ! needed to filter spatail map + +!Retain init where tile_frac=0 +LAI_pft_temp(:,:) = 0.0 +HGT_pft_temp(:,:) = 0.0 + +!Comensurate with restricted ordering of PFTs in UM. Not generically applicable +DO n=1,ntiles + DO j=1,tile_pts(n) + + i = tile_index(j,n) ! landpt index + + IF( tile_frac(i,n) .GT. 0.0 ) THEN + + IF( n < shrub_cable ) THEN ! trees + + LAI_pft_temp(i,n) = MAX( 0.99*CLAI_thresh, LAI_pft(i,n) ) + HGT_pft_temp(i,n) = MAX( 1.0, HGT_pft(i,n) ) + + ELSE IF( n >= shrub_cable .AND. n < aust_temperate ) THEN ! shrubs/grass + + LAI_pft_temp(i,n) = MAX( 0.99*CLAI_thresh, LAI_pft(i,n) ) + HGT_pft_temp(i,n) = MAX( 0.1, HGT_pft(i,n) ) + + ELSE IF( n == aust_temperate .OR. n == aust_tropical ) THEN ! Aust. trees + + LAI_pft_temp(i,n) = MAX( CLAI_thresh, LAI_pft(i,n) ) + HGT_pft_temp(i,n) = MAX( 1.0, HGT_pft(i,n) ) + + ENDIF + + ENDIF + + ENDDO +ENDDO + +!surface_type = PACK(surface_type_temp, um1%L_TILE_PTS) +LAI_pft_cbl = PACK(LAI_pft_temp, l_tile_pts) +HGT_pft_cbl = PACK(HGT_pft_temp, l_tile_pts) + +END SUBROUTINE limit_HGT_LAI + +END MODULE cbl_LAI_canopy_height_mod + diff --git a/src/coupled/esm/cable_surface_types.F90 b/src/coupled/esm/cable_surface_types.F90 index f68d5bf36..4a0f2244b 100644 --- a/src/coupled/esm/cable_surface_types.F90 +++ b/src/coupled/esm/cable_surface_types.F90 @@ -24,8 +24,8 @@ MODULE cable_surface_types_mod INTEGER, PARAMETER :: c3_cropland = 9 INTEGER, PARAMETER :: c4_cropland = 10 INTEGER, PARAMETER :: wetland = 11 - INTEGER, PARAMETER :: empty1 = 12 - INTEGER, PARAMETER :: empty2 = 13 + INTEGER, PARAMETER :: aust_temperate = 12 + INTEGER, PARAMETER :: aust_tropical = 13 INTEGER, PARAMETER :: barren_cable = 14 INTEGER, PARAMETER :: urban_cable = 15 INTEGER, PARAMETER :: lakes_cable = 16