diff --git a/CMakeLists.txt b/CMakeLists.txt index 693712860..93d150c94 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,7 +84,7 @@ add_library( src/science/canopy/cbl_dryLeaf.F90 src/science/canopy/cbl_friction_vel.F90 src/science/canopy/cbl_fwsoil.F90 - src/science/canopy/cbl_init_wetfac_mod.f90 + src/science/canopy/cbl_init_wetfac_mod.F90 src/science/canopy/cbl_latent_heat.F90 src/science/canopy/cbl_photosynthesis.F90 src/science/canopy/cbl_pot_evap_snow.F90 diff --git a/src/science/canopy/cable_canopy.F90 b/src/science/canopy/cable_canopy.F90 index 7927c3da7..465a1d1e8 100644 --- a/src/science/canopy/cable_canopy.F90 +++ b/src/science/canopy/cable_canopy.F90 @@ -6,24 +6,22 @@ MODULE cable_canopy_module PRIVATE CONTAINS - + SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,climate, sunlit_veg_mask, reducedLAIdue2snow ) - USE cable_def_types_mod - USE cbl_radiation_module, ONLY : radiation - USE cable_air_module - USE cable_common_module - USE cable_roughness_module - -USE cbl_friction_vel_module, ONLY : comp_friction_vel, psim, psis -USE cbl_pot_evap_snow_module, ONLY : Penman_Monteith, Humidity_deficit_method -USE cbl_qsat_module, ONLY : qsatfjh, qsatfjh2 -USE cbl_zetar_module, ONLY : update_zetar -USE cable_latent_heat_module, ONLY : latent_heat_flux -USE cable_wetleaf_module, ONLY : wetleaf -USE cbl_dryLeaf_module, ONLY : dryLeaf +! subrs +USE cbl_radiation_module, ONLY : radiation +USE cbl_friction_vel_module, ONLY : comp_friction_vel, psim, psis +USE cbl_pot_evap_snow_module, ONLY : Penman_Monteith, Humidity_deficit_method +USE cbl_qsat_module, ONLY : qsatfjh, qsatfjh2 +USE cbl_zetar_module, ONLY : update_zetar +USE cable_latent_heat_module, ONLY : latent_heat_flux +USE cable_wetleaf_module, ONLY : wetleaf +USE cbl_dryLeaf_module, ONLY : dryLeaf USE cable_within_canopy_module, ONLY : within_canopy USE cbl_SurfaceWetness_module, ONLY : Surf_wetness_fact +! data +USE cable_common_module, ONLY: cable_runtime, cable_user ! physical constants USE cable_phys_constants_mod, ONLY : CTFRZ => TFRZ USE cable_phys_constants_mod, ONLY : CRMAIR => RMAIR @@ -65,6 +63,12 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima USE cable_math_constants_mod, ONLY : CPI_C => PI USE cable_other_constants_mod, ONLY : CLAI_THRESH => LAI_THRESH +USE grid_constants_mod_cbl, ONLY: ICE_SoilType + + USE cable_def_types_mod + USE cable_air_module + USE cable_roughness_module + IMPLICIT NONE TYPE (balances_type), INTENT(INOUT) :: bal @@ -80,7 +84,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima TYPE (veg_parameter_type), INTENT(INOUT) :: veg REAL :: reducedLAIdue2snow(mp) -logical :: sunlit_veg_mask(mp) +LOGICAL :: sunlit_veg_mask(mp) REAL, INTENT(IN) :: dels ! integration time setp (s) INTEGER :: & iter, & ! iteration # @@ -151,16 +155,21 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima INTEGER :: j - INTEGER, SAVE :: call_number =0 - - ! END header - - call_number = call_number + 1 - - !H! vNot sure that this is appropriate for JULES standalone - HaC - !H!IF( .NOT. cable_runtime%um) & - canopy%cansto = canopy%oldcansto - +! local vars in calc of Surface conductance +REAL :: LAI_min(mp) +REAL :: Rel_bigLeafLAI_sun(mp) +REAL :: Rel_bigLeafLAI_shd(mp) +REAL :: minCanopyCond(mp) +REAL :: canopy_conductance(mp) +REAL :: Rel_moisture(mp) +REAL :: soil_conductance(mp) +REAL :: Surf_conductance(mp) +! END header + + ! Not sure that this is appropriate for JULES standalone - HaC either + IF( .NOT. cable_runtime%um ) THEN + canopy%cansto = canopy%oldcansto + ENDIF ALLOCATE( cansat(mp), gbhu(mp,mf)) ALLOCATE( dsx(mp), fwsoil(mp), tlfx(mp), tlfy(mp) ) ALLOCATE( ecy(mp), hcy(mp), rny(mp)) @@ -394,23 +403,17 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima hcy = 0.0 ! init current estimate lat heat ecy = rny - hcy ! init current estimate lat heat - CALL dryLeaf( dels, rad, rough, air, met, & - veg, canopy, soil, ssnow, dsx, & - fwsoil, tlfx, tlfy, ecy, hcy, & - rny, gbhu, gbhf, csx, cansat, & - ghwet, iter,climate, sum_rad_gradis, sum_rad_rniso ) + CALL dryLeaf( dels, rad, rough, air, met, veg, canopy, soil, ssnow, & + dsx, fwsoil, tlfx, tlfy, ecy, hcy, rny, gbhu, gbhf, csx, & + cansat, ghwet, iter, climate, sum_rad_gradis, & + sum_rad_rniso ) - - CALL wetLeaf( dels, & - cansat, tlfy, & - gbhu, gbhf, ghwet, & - mp, CLAI_thresh, CCAPP, CRmair, & - reducedLAIdue2snow, sum_rad_rniso, sum_rad_gradis, & - canopy%fevw, canopy%fevw_pot, canopy%fhvw, & - canopy%fwet, canopy%cansto, air%rlam, air%dsatdk, & + CALL wetLeaf( dels, cansat, tlfy, gbhu, gbhf, ghwet, mp, CLAI_thresh, & + CCAPP, CRmair, reducedLAIdue2snow, sum_rad_rniso, & + sum_rad_gradis, canopy%fevw, canopy%fevw_pot, canopy%fhvw, & + canopy%fwet, canopy%cansto, air%rlam, air%dsatdk, & met%tvair, met%tk, met%dva, air%psyc ) - ! Calculate latent heat from vegetation: ! Calculate sensible heat from vegetation: ! Calculate net rad absorbed by canopy: @@ -688,16 +691,36 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima canopy%cduv = canopy%us * canopy%us / (MAX(met%ua,CUMIN))**2 - !---diagnostic purposes - canopy%gswx_T = rad%fvlai(:,1)/MAX( CLAI_THRESH, canopy%vlaiw(:) ) & - * canopy%gswx(:,1) + rad%fvlai(:,2) / MAX(CLAI_THRESH, & - canopy%vlaiw(:))*canopy%gswx(:,2) +! Evaluate Total Surface conductance !---diagnostic purposes + +!vlaiw is already limted? +LAI_min(:) = MAX( CLAI_THRESH, canopy%vlaiw(:) ) +Rel_bigLeafLAI_sun(:) = rad%fvlai(:,1) / LAI_min(:) +Rel_bigLeafLAI_shd(:) = rad%fvlai(:,2) / LAI_min(:) - ! The surface conductance below is required by dust scheme; it is composed from canopy and soil conductances - canopy%gswx_T = (1.-rad%transd)*MAX(1.e-06,canopy%gswx_T ) + & !contribution from canopy conductance - rad%transd*(.01*ssnow%wb(:,1)/soil%sfc)**2 ! + soil conductance; this part is done as in Moses - WHERE ( soil%isoilm == 9 ) canopy%gswx_T = 1.e6 ! this is a value taken from Moses for ice points +!total stomatal conductance +canopy_conductance(:) = Rel_bigLeafLAI_sun(:) * canopy%gswx(:,1) & + + Rel_bigLeafLAI_shd(:) * canopy%gswx(:,2) +! Surface conductance required by dust scheme +! Canopy conductance +minCanopyCond(:) = MAX( 1.e-06, canopy_conductance(:) ) +canopy_conductance(:) = (1.-rad%transd(:))* minCanopyCond(:) + +! Soil conductance - following MOSES method +Rel_moisture(:) = ssnow%wb(:,1) / soil%sfc(:) +soil_conductance(:) = rad%transd(:) * ( 0.01*Rel_moisture )**2 + +! Combined Surface conductance +Surf_conductance(:) = canopy_conductance(:) + soil_conductance(:) + +! this is a value taken from MOSES for ice points +WHERE ( soil%isoilm == ICE_SoilType ) + Surf_conductance = 1.e6 +END WHERE + +canopy%gswx_T(:) = Surf_conductance(:) !fill CABLE type for now + canopy%cdtq = canopy%cduv * & ( LOG( rough%zref_uv / rough%z0m) - & psim( canopy%zetar(:,NITER) * rough%zref_uv/rough%zref_tq, mp, CPI_C ) & diff --git a/src/science/canopy/cbl_SurfaceWetness.F90 b/src/science/canopy/cbl_SurfaceWetness.F90 index bbc530312..83cf4731a 100644 --- a/src/science/canopy/cbl_SurfaceWetness.F90 +++ b/src/science/canopy/cbl_SurfaceWetness.F90 @@ -66,13 +66,17 @@ SUBROUTINE Surf_wetness_fact( cansat, canopy, ssnow,veg, met, soil, dels ) !call saturated_fraction(ssnow,soil,veg) ssnow%satfrac(:) = 1.0e-8 ssnow%rh_srf(:) = 1.0 - - CALL initialize_wetfac(mp, ssnow%wetfac, soil%swilt, soil%sfc, ssnow%wb, ssnow%wbice, 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) + + !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) diff --git a/src/science/canopy/cbl_dryLeaf.F90 b/src/science/canopy/cbl_dryLeaf.F90 index c9d5dea8c..f71f299fe 100644 --- a/src/science/canopy/cbl_dryLeaf.F90 +++ b/src/science/canopy/cbl_dryLeaf.F90 @@ -18,6 +18,11 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & 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 @@ -26,6 +31,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & 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 @@ -35,6 +41,8 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & 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 TYPE (air_type), INTENT(INOUT) :: air @@ -145,7 +153,9 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & ! END header ALLOCATE( gswmin(mp,mf )) - + + gs_coeff(:,:) = 0.0 ! CABLE_ISSUE39 + ! Soil water limitation on stomatal conductance: IF( iter ==1) THEN IF ((cable_user%soil_struc=='default').AND.(cable_user%FWSOIL_SWITCH.NE.'Haverd2013')) THEN @@ -319,36 +329,38 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & !1.2877 + (0.0116 × Vcmax,a25) – (0.0334 × TWQ) !Shrubs: Rdark,a25 = 1.5758 + (0.0116 × Vcmax,a25) – (0.0334 × TWQ) - IF (veg%iveg(i).EQ.2 .OR. veg%iveg(i).EQ. 4 ) THEN ! broadleaf forest - - rdx(i,1) = 0.60*(1.2818e-6+0.0116*veg%vcmax(i)- & - 0.0334*climate%qtemp_max_last_year(i)*1e-6) - rdx(i,2) = rdx(i,1) - - ELSEIF (veg%iveg(i).EQ.1 .OR. veg%iveg(i).EQ. 3 ) THEN ! needleleaf forest - rdx(i,1) = 1.0*(1.2877e-6+0.0116*veg%vcmax(i)- & - 0.0334*climate%qtemp_max_last_year(i)*1e-6) - rdx(i,2) = rdx(i,1) - - ELSEIF (veg%iveg(i).EQ.6 .OR. veg%iveg(i).EQ.8 .OR. & - veg%iveg(i).EQ. 9 ) THEN ! C3 grass, tundra, crop - rdx(i,1) = 0.60*(1.6737e-6+0.0116*veg%vcmax(i)- & - 0.0334*climate%qtemp_max_last_year(i)*1e-6) - rdx(i,2) = rdx(i,1) - - ELSE ! shrubs and other (C4 grass and crop) - rdx(i,1) = 0.60*(1.5758e-6+0.0116*veg%vcmax(i)- & - 0.0334*climate%qtemp_max_last_year(i)*1e-6) - rdx(i,2) = rdx(i,1) + ! broadleaf forest + IF ( veg%iveg(i) .EQ. evergreen_broadleaf .OR. & + veg%iveg(i) .EQ. deciduous_broadleaf ) THEN + + rdx(i,:) = 0.60 * ( 1.2818e-6 + 0.0116 * veg%vcmax(i) - & + 0.0334 * climate%qtemp_max_last_year(i) * 1.0e-6 ) + + ! needleleaf forest + ELSEIF ( veg%iveg(i) .EQ. evergreen_needleleaf .OR. & + veg%iveg(i) .EQ. deciduous_needleleaf ) THEN + + rdx(i,:) = 1.0 * ( 1.2877e-6 + 0.0116 * veg%vcmax(i) - & + 0.0334 * climate%qtemp_max_last_year(i) * 1.0e-6 ) + + ! C3 grass, tundra , C3 crop + ELSEIF ( veg%iveg(i) .EQ. c3_grassland .OR. & + veg%iveg(i) .EQ. tundra .OR. & + veg%iveg(i) .EQ. c3_cropland ) THEN + + rdx(i,:) = 0.60 * ( 1.6737e-6 + 0.0116 * veg%vcmax(i) - & + 0.0334 * climate%qtemp_max_last_year(i) * 1e-6 ) + + ! shrub & other (C4 grass and C4 crop) (wetlands, nveg) TBC + ELSE + rdx(i,:) = 0.60 * ( 1.5758e-6 + 0.0116 * veg%vcmax(i) - & + 0.0334 * climate%qtemp_max_last_year(i) * 1.0e-6 ) ENDIF - ! modify for leaf area and instanteous temperature response (Rd25 -> Rd) rdx(i,1) = rdx(i,1) * xrdt(tlfx(i)) * rad%scalex(i,1) rdx(i,2) = rdx(i,2) * xrdt(tlfx(i)) * rad%scalex(i,2) - - ! reduction of daytime leaf dark-respiration to account for !photo-inhibition !Mercado, L. M., Huntingford, C., Gash, J. H. C., Cox, P. M., @@ -369,15 +381,8 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & rdx(i,2) = rdx(i,2) * & (0.5 - 0.05*LOG(jtomol*1.0e6*rad%qcan(i,1,2))) -!$ xleuning(i,1) = ( fwsoil(i) / ( csx(i,1) - co2cp3 ) ) & -!$ * ( veg%a1gs(i) / ( 1.0 + dsx(i)/veg%d0gs(i))) -!$ xleuning(i,2) = ( fwsoil(i) / ( csx(i,2) - co2cp3 ) ) & -!$ * ( veg%a1gs(i) / ( 1.0 + dsx(i)/veg%d0gs(i))) - ELSE !cable_user%call_climate -!$!Vanessa:note there is no xleuning to go into photosynthesis etc anymore -!$ gs_coeff = xleuning rdx(i,1) = (veg%cfrd(i)*vcmxt3(i,1) + veg%cfrd(i)*vcmxt4(i,1)) rdx(i,2) = (veg%cfrd(i)*vcmxt3(i,2) + veg%cfrd(i)*vcmxt4(i,2)) @@ -494,12 +499,12 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & REAL(veg%froot(i,:),r_2),& soil%ssat_vec(i,:), & soil%swilt_vec(i,:), & - MAX(REAL(canopy%fevc(i)/air%rlam(i)/1000_r_2,r_2),0.0_r_2), & + MAX(REAL(canopy%fevc(i)/air%rlam(i)/density_liq,r_2),0.0_r_2), & REAL(veg%gamma(i),r_2), & REAL(soil%zse,r_2), REAL(dels,r_2), REAL(veg%zr(i),r_2)) fwsoil(i) = canopy%fwsoil(i) - ssnow%evapfbl(i,:) = ssnow%rex(i,:)*dels*1000_r_2 ! mm water & + ssnow%evapfbl(i,:) = ssnow%rex(i,:)*dels*density_liq ! mm water & !(root water extraction) per time step ELSE @@ -513,7 +518,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met, & 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) * 1000.0 ) + soil%zse(kk) * density_liq ) ENDDO IF (cable_user%soil_struc=='default') THEN diff --git a/src/science/canopy/cbl_init_wetfac_mod.f90 b/src/science/canopy/cbl_init_wetfac_mod.F90 similarity index 100% rename from src/science/canopy/cbl_init_wetfac_mod.f90 rename to src/science/canopy/cbl_init_wetfac_mod.F90 diff --git a/src/science/canopy/cbl_latent_heat.F90 b/src/science/canopy/cbl_latent_heat.F90 index 040828830..aed415afa 100644 --- a/src/science/canopy/cbl_latent_heat.F90 +++ b/src/science/canopy/cbl_latent_heat.F90 @@ -52,7 +52,8 @@ SUBROUTINE Latent_heat_flux( mp, CTFRZ, dels, soil_zse, soil_swilt, & ! over snow. And the value for points with new snow could be smaller. -USE cable_def_types_mod, ONLY : r_2 +USE cable_def_types_mod, ONLY : r_2 + IMPLICIT NONE INTEGER :: mp diff --git a/src/science/canopy/cbl_photosynthesis.F90 b/src/science/canopy/cbl_photosynthesis.F90 index b9ff57450..8d4caddfe 100644 --- a/src/science/canopy/cbl_photosynthesis.F90 +++ b/src/science/canopy/cbl_photosynthesis.F90 @@ -44,6 +44,11 @@ SUBROUTINE photosynthesis( csxz, cx1z, cx2z, gswminz, & INTEGER :: i,j + anrubpz(:,:) = 0.0 + ansinkz(:,:) = 0.0 + anxz(:,:) = 0.0 + anrubiscoz(:,:) = 0.0 + DO i=1,mp IF (SUM(vlaiz(i,:)) .GT. CLAI_THRESH) THEN