Here are the modified subroutines within JULES v2.1, which we had permission to provide in isolation. Permission and license-information to use the whole JULES model is available from http://jules- lsm.github.io/access_req/JULES_access.html , from where the complete original JULES v2.1 model may be requested by email. Others are allowed to freely use our modifications, giving appropriate credit, but the JULES conditions referred to above still apply. Ideally it should be a case of copying and pasting, overwriting, these subroutines into the locations of the same-name subroutines in JULES v2.1, except for 'smc_sph' which needs a new file creating and registering or copied into the same file as ‘smc_ext’ (not overwriting in this last case). Problems may be encountered to do with local-variable-scoping and local variable would need to be added to some subroutines to pass parameters between the modified functions within the existing JULES structure. Additionally, your JULES config file will need a new PFT type creating (call it ‘sphag’ for instance – this is what I have assumed it’s called in the attached code) with the parameters as given in the paper Coppell, Gloor & Holden. Four subroutines are modified – ‘physiol’, ‘smc_sph’ (new subroutine), ‘leaf_limits’ and ‘leaf’. ‘Physiol’ calculates model gridbox fluxes, including, specific to this case, evapotranspiration and water transport, and calculating the new desiccation-stress function as an output parameter (called soil moisture availability factor in the code – fsmc) from Coppell, Gloor & Holden Table 2 point 1. ‘physiol’ directly calls ‘smc_sph’ ‘leaf_limits’ simulates leaf resistance and the C3 photosynthesis processes. It takes fsmc (from ‘physiol’ above) as one of its input parameters. It incorporates the changes in Table 2 point 3 for the photosynthesis mechanism in Sphagnum. ‘Leaf’ calculates final net leaf photosynthesis using Cox’s smoothed- minimum function, and also leaf surface conductance, which for this sphagnum model is hardcoded as a single value. ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! Description: ! Subroutine to calculate gridbox mean values of surface conductance ! and carbon fluxes. Also returns net primary productivity, leaf ! turnover and wood respiration of each plant functional type for ! driving TRIFFID. !********************************************************************** SUBROUTINE physiol (row_length,rows,land_pts,land_index & , nshyd,ntiles,tile_pts,tile_index,l_aggregate & , dim_cs1, dim_cs2 & , co2,co2_3d,co2_dim_len & , co2_dim_row,l_co2_interactive & , l_triffid, l_q10 & , can_model,cs,frac,ht,ipar,lai,pstar,q1 & , sthu,tsoil,tstar_tile & , v_crit,v_sat,v_wilt,wind,z0_tile,z1 & , canhc_tile,vfrac_tile,emis_tile,emis_soil & , flake,g_leaf,gs,gs_tile & , gpp,gpp_ft,npp,npp_ft,resp_p,resp_p_ft & , resp_s,resp_w_ft,smct,wt_ext_tile,fsmc,wt_ext & , ra,albsoil,cos_zenith_angle & , can_rad_mod,ilayers) USE c_densty USE nstypes, ONLY : & ! imported scalars with intent(in) lake,npft,ntype,soil USE nvegparm USE soil_param, ONLY : dzsoil USE surf_param, ONLY : diff_frac USE pftparm IMPLICIT NONE INTEGER & row_length & ! IN Number of points on a row ,rows & ! IN Number of rows in a theta field ,land_pts & ! IN Number of land points to be ! ! processed. ,land_index(land_pts) & ! IN Index of land points on the ! ! P-grid. ,co2_dim_len & ! IN Length of a CO2 field row. ,co2_dim_row & ! IN Number of CO2 field rows. ,nshyd & ! IN Number of soil moisture ! ! levels. ,ntiles & ! IN Number of surface tiles. ,tile_pts(ntype) & ! IN Number of land points which ! ! include the nth surface type. ,tile_index(land_pts,ntype) & ! IN Indices of land points which ! ! include the nth surface type. ,can_model & ! IN Swith for thermal vegetation ! ! canopy ,dim_cs1, dim_cs2 ! soil carbon dimensions LOGICAL & l_co2_interactive & ! switch for 3D CO2 field , l_aggregate & ! IN Logical to set aggregate ! surface scheme , l_triffid & ! TRUE if using TRIFFID , l_q10 ! TRUE if using Q10 for soil resp INTEGER & can_rad_mod & ! !Switch for canopy radiation model ,ilayers ! !No of layers in canopy radiation model REAL & alb_type_dummy(land_pts,ntype,4) & ! ! WORK Dummy argument for albedo ! ! subroutine ,fapar_dir(land_pts,npft,ilayers) & ! ! WORK Profile of absorbed PAR - ! ! Direct beam ,fapar_dif(land_pts,npft,ilayers) & ! ! WORK Profile of absorbed PAR - ! ! Diffuse beam ,fapar_dir_tot(land_pts,npft) & ! ! WORK Total absorbed PAR - ! ! Direct beam ,fapar_dif_tot(land_pts,npft) & ! ! WORK Total absorbed PAR - ! ! Diffuse beam ,faparv(land_pts,ilayers) ! WORK Profile of absorbed PAR. REAL & co2 & ! IN Atmospheric CO2 concentration ,co2_3d(co2_dim_len,co2_dim_row) & ! ! IN 3D atmos CO2 concentration ! ! (kg CO2/kg air). ,cs(land_pts,dim_cs1) & ! IN Soil carbon (kg C/m2). ,veg_frac(dim_cs2) & ! WORK vegetated fraction of gridbox ,frac(land_pts,ntype) & ! IN Surface type fractions. ,ht(land_pts,npft) & ! IN Canopy height (m). ,ipar(row_length,rows) & ! IN Incident PAR (W/m2). ,lai(land_pts,npft) & ! IN Leaf area index. ,pstar(row_length,rows) & ! IN Surface pressure (Pa). ,q1(row_length,rows) & ! IN Specific humidity at level 1 ! ! (kg H2O/kg air). ,sthu(land_pts,nshyd) & ! IN Soil moisture content in each ! ! layer as a fraction of saturation ,tsoil(land_pts) & ! IN Soil temperature (K). ,tstar_tile(land_pts,ntiles) & ! IN Tile surface temperatures (K). ,v_crit(land_pts,nshyd) & ! IN Volumetric soil moisture ! ! concentration above which ! ! stomata are not sensitive ! ! to soil water (m3 H2O/m3 soil). ,v_sat(land_pts,nshyd) & ! IN Volumetric soil moisture ! ! concentration at saturation ! ! (m3 H2O/m3 soil). ,v_wilt(land_pts,nshyd) & ! IN Volumetric soil moisture ! ! concentration below which ! ! stomata close (m3 H2O/m3 soil). ,wind(row_length,rows) & ! IN Windspeed (m/s). ,z0_tile(land_pts,ntiles) & ! IN Tile roughness lengths (m). ,z1(row_length,rows) & ! IN Windspeed reference height(m). ,gs(land_pts) & ! INOUT Gridbox mean surface ! ! conductance (m/s). ,wt_ext(land_pts,nshyd) & ! OUT Gridbox-mean WT_EXT. ! NB This is only non-zero if ! l_aggregate=TRUE. ,ra(land_pts) & ! OUT Aerodynamic resistance (s/m). ,albsoil(land_pts) & ! ! Soil albedo. ,cos_zenith_angle(row_length, rows) ! ! Cosine of the zenith angle REAL & canhc_tile(land_pts,ntiles) & ! OUT Areal heat capacity of canopy ! ! for land tiles (J/K/m2). ,flake(land_pts,ntiles) & ! OUT Lake fraction. ,g_leaf(land_pts,npft) & ! OUT Leaf turnover rate (/360days). ,gs_tile(land_pts,ntiles) & ! OUT Surface conductance for ! ! land tiles (m/s). ,gpp(land_pts) & ! OUT Gridbox mean gross primary ! ! productivity (kg C/m2/s). ,gpp_ft(land_pts,npft) & ! OUT Gross primary productivity ! ! (kg C/m2/s). ,npp(land_pts) & ! OUT Gridbox mean net primary ! ! productivity (kg C/m2/s). ,npp_ft(land_pts,npft) & ! OUT Net primary productivity ! ! (kg C/m2/s). ,resp_p(land_pts) & ! OUT Gridbox mean plant respiration ! ! (kg C/m2/s). ,resp_p_ft(land_pts,npft) & ! OUT Plant respiration (kg C/m2/s). ,resp_s(land_pts,dim_cs1) & ! OUT Soil respiration (kg C/m2/s). ,resp_w_ft(land_pts,npft) & ! OUT Wood maintenance respiration ! ! (kg C/m2/s). ,smct(land_pts) & ! OUT Available moisture in the ! ! soil profile (mm). ,vfrac_tile(land_pts,ntiles) & ! OUT Fractional canopy coverage for ! ! land tiles. ,emis_tile(land_pts,ntiles) & ! OUT Emmissivity for land tiles ,emis_soil(land_pts) & ! OUT Emmissivity of underlying soil ,wt_ext_tile(land_pts,nshyd,ntiles) & ! ! OUT Fraction of evapotranspiration ! ! which is extracted from each ! ! soil layer by each tile. ,fsmc(land_pts,npft) ! OUT Moisture availability factor. ! External routines called :- EXTERNAL root_frac,smc_ext,raero,sf_stom,soil_evap, & leaf_lit,cancap,microbe REAL & canhc(land_pts) & ! WORK Canopy heat capacity (J/K/m2). ,ch_type(land_pts,ntype) & ! WORK CANHC for surface types. ,f_root(nshyd) & ! WORK Fraction of roots in each soil ! ! layer. ,gsoil(land_pts) & ! WORK Bare soil conductance. ,gs_type(land_pts,ntype) & ! WORK Conductance for surface types. ,pstar_land(land_pts) & ! WORK Surface pressure (Pa). ,rib(row_length,rows) & ! WORK Bulk Richardson Number. ,tstar(land_pts) & ! WORK Surface temperature (K). ,vfrac(land_pts) & ! WORK Fractional canopy coverage. ,vf_type(land_pts,ntype) & ! WORK VFRAC for surface types. ,wt_ext_type(land_pts,nshyd,ntype) & ! ! WORK WT_EXT for surface types. ,fsoil(land_pts,npft) & ! WORK Fraction of ground below canopy ! ! contributing to evaporation. ,fsoil_tot(land_pts) & ! WORK Total fraction of soil ! ! contributing to evaporation ! ! ( = bare soil fraction ! ! + fraction of ground below canopy ,z0(land_pts) & ! WORK Roughness length (m). ,WTP_RIUTTA !WORK - RIUTTA WATER TABLE FOR FIXED LAB TYPE SPHAG CALIBRATION INTEGER & i,j,k,l,m,n ! Loop indices CHARACTER*50 char_debug !----------------------------------------------------------------------- ! Initialisations !----------------------------------------------------------------------- f_root(:)=0.0 gpp(:)=0.0 npp(:)=0.0 resp_p(:)=0.0 smct(:)=0.0 ra(:)=0.0 canhc(:)=0.0 vfrac(:)=0.0 veg_frac(:)=0.0 wt_ext(:,:)=0.0 rib(:,:)=0.0 g_leaf(:,:)=0.0 gpp_ft(:,:)=0.0 npp_ft(:,:)=0.0 resp_p_ft(:,:)=0.0 resp_w_ft(:,:)=0.0 fsmc(:,:)=0.0 resp_s(:,:)=0.0 faparv(:,:)=0.0 wt_ext_type(:,:,:)=0.0 alb_type_dummy(:,:,:)=0.0 fapar_dir(:,:,:)=0.0 fapar_dif(:,:,:)=0.0 !****************************************************************** ! Sphagnum hard-coded test variables for running laboratory-equivalent ! decoupling from climate signal !****************************************************************** !ipar(:,:)=440.0 !Riutta sphag ps set PAR ipar(:,:)=330.0 !set for 'full daylight' Strack 1250 um m-2 s-1 !tstar(:)=273.15+13 !Riutta respiration set temp !tstar(:)=273.15+22.81 !Topt for sphag ps from Riutta - 22.81 (3.05 S.E.); sedges - 29.27 (3.34); sedges+shrubs - 24.88 (0.85) !tstar(:)=273.15+12.6 !Strack mean daily summer temp tstar(:)=273.15+25 !w&f STANDARDISED TEMP !tstar_tile(:,:)=273.15+12.6 !copied from tstar, so amend when tstar amended tstar_tile(:,:)=273.15+25 !copied from tstar, so amend when tstar amended !WTP_RIUTTA=0.14 !Riutta respiration set WL WTP_RIUTTA=0.0854 !WL opt for sphag ps for Riutta - 0.0854(0.0137); sedges - 0.16 (0.0116); sedges and shrubs 0.2910 (0.0948) DO l=1,land_pts j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length pstar_land(l) = pstar(i,j) END DO DO n=1,ntype DO l=1,land_pts gs_type(l,n)=gs(l) END DO END DO gs(:)=0.0 DO l=1,land_pts fsoil_tot(l) = frac(l,soil) gsoil(l) = 0. IF (v_crit(l,1) > 0.) & gsoil(l) = gs_nvg(soil - npft) * & (sthu(l,1) * v_sat(l,1) / v_crit(l,1))**2 END DO !----------------------------------------------------------------------- ! Calculate light absorption by the plant canopy !----------------------------------------------------------------------- IF ( can_rad_mod /= 1 ) THEN ! DEPENDS ON: albpft CALL albpft (row_length * rows,land_pts, & land_index,tile_index,tile_pts,ilayers, & albsoil,cos_zenith_angle,lai,alb_type_dummy, & fapar_dir,fapar_dif) END IF !----------------------------------------------------------------------- ! Loop over Plant Functional Types to calculate the available moisture ! and the values of canopy conductance, the carbon fluxes and the leaf ! turnover rate !----------------------------------------------------------------------- DO n=1,npft IF ( l_aggregate ) THEN DO l=1,land_pts tstar(l) = tstar_tile(l,1) z0(l) = z0_tile(l,1) END DO ELSE DO l=1,land_pts tstar(l) = tstar_tile(l,n) z0(l) = z0_tile(l,n) END DO END IF ! DEPENDS ON: root_frac CALL root_frac(nshyd,dzsoil,rootd_ft(n),f_root) !************************************************** ! Sphagnum PFT !************************************************** IF (pftName(N).EQ.'sphag') THEN CALL smc_sph (land_pts,nshyd,tile_pts(n),tile_index(:,n) & , f_root ,sthu,v_crit,v_sat,v_wilt & , wt_ext_type(:,:,n),fsmc(:,n), WTP_RIUTTA) ELSE ! all the other PFT types ! DEPENDS ON: smc_ext CALL smc_ext (land_pts,nshyd,tile_pts(n),tile_index(:,n) & , f_root,sthu,v_crit,v_sat,v_wilt & , wt_ext_type(:,:,n),fsmc(:,n)) END IF ! DEPENDS ON: raero CALL raero (row_length,rows,land_pts,land_index & , tile_pts(n),tile_index(:,n) & , rib,wind,z0,z0,z1,ra) !----------------------------------------------------------------------- ! Calculate light absorption by the plant canopy !----------------------------------------------------------------------- IF (can_rad_mod /= 1) THEN DO l=1,land_pts fapar_dif_tot(l,n)=0.0 fapar_dir_tot(l,n)=0.0 END DO DO m=1,ilayers DO k=1,tile_pts(n) l = tile_index(k,n) i = land_index(l) fapar_dir_tot(l,n) = fapar_dir_tot(l,n) & + fapar_dir(l,n,m) fapar_dif_tot(l,n) = fapar_dif_tot(l,n) & + fapar_dif(l,n,m) faparv(l,m) = (1-diff_frac(i)) * fapar_dir(l,n,m) & + diff_frac(i) * fapar_dif(l,n,m) END DO END DO END IF ! !write(char_debug, '(7I)') n ! !write(6,*) 'n(pft#)= '//char_debug ! DEPENDS ON: sf_stom CALL sf_stom (row_length,rows,land_pts,land_index & , tile_pts(n),tile_index(:,n),n & , co2,co2_3d,co2_dim_len & , co2_dim_row,l_co2_interactive & , fsmc(:,n),ht(:,n),ipar,lai(:,n),pstar_land & , q1,ra,tstar & , can_rad_mod,ilayers,faparv & , gpp_ft(:,n),npp_ft(:,n),resp_p_ft(:,n) & , resp_w_ft(:,n),gs_type(:,n)) ! DEPENDS ON: soil_evap CALL soil_evap (land_pts,nshyd,tile_pts(n),tile_index(:,n) & , gsoil,lai(:,n),gs_type(:,n),wt_ext_type(:,:,n) & , fsoil(:,n)) ! DEPENDS ON: leaf_lit CALL leaf_lit (land_pts,tile_pts(n),tile_index(:,n) & , n,fsmc(:,n),tstar,g_leaf(:,n)) ! DEPENDS ON: cancap CALL cancap (land_pts,tile_pts(n),tile_index(:,n),can_model,n & , ht(:,n),lai(:,n),ch_type(:,n),vf_type(:,n)) DO l=1,land_pts fsoil_tot(l) = fsoil_tot(l) + frac(l,n)*fsoil(l,n) END DO END DO !---------------------------------------------------------------------- ! Non-vegetated surface types !---------------------------------------------------------------------- DO n=npft+1,ntype DO m=1,tile_pts(n) l=tile_index(m,n) gs_type(l,n) = gs_nvg(n-npft) END DO END DO ! Copy soil conductance and add bare soil fraction to extraction from ! surface layer n = soil DO m=1,tile_pts(n) l=tile_index(m,n) gs_type(l,n) = gsoil(l) wt_ext_type(l,1,n) = 1. END DO !---------------------------------------------------------------------- ! Canopy heat capacity and coverage for non-vegetated surfaces !---------------------------------------------------------------------- DO n=npft+1,ntype DO m=1,tile_pts(n) l=tile_index(m,n) ch_type(l,n) = ch_nvg(n-npft) vf_type(l,n) = vf_nvg(n-npft) END DO END DO !---------------------------------------------------------------------- ! Surface emissivity !---------------------------------------------------------------------- IF ( l_aggregate ) THEN emis_tile(:,:) = 0.0 DO n=1,npft DO m=1,tile_pts(n) l = tile_index(m,n) emis_tile(l,1) = emis_tile(l,1) + frac(l,n) * emis_pft(n) END DO END DO DO n=npft+1,ntype DO m=1,tile_pts(n) l = tile_index(m,n) emis_tile(l,1) = emis_tile(l,1) & + frac(l,n) * emis_nvg(n - npft) END DO END DO ELSE DO n=1,npft DO m=1,tile_pts(n) l = tile_index(m,n) emis_tile(l,n) = emis_pft(n) END DO END DO DO n=npft+1,ntype DO m=1,tile_pts(n) l = tile_index(m,n) emis_tile(l,n) = emis_nvg(n - npft) END DO END DO END IF DO l=1,land_pts emis_soil(l) = emis_nvg(soil - npft) END DO !---------------------------------------------------------------------- ! Calculate the rate of soil respiration !---------------------------------------------------------------------- ! set VEG_FRAC according to whether it is full or dummy field IF (l_triffid) THEN DO l=1,land_pts veg_frac(l) = SUM(frac(l,1:npft)) END DO END IF ! DEPENDS ON: microbe CALL microbe (land_pts,dim_cs1,dim_cs2,l_triffid,l_q10,cs, & sthu(:,1),v_sat(:,1),v_wilt(:,1),tsoil,resp_s, & veg_frac) !---------------------------------------------------------------------- ! Form gridbox mean values !---------------------------------------------------------------------- DO n=1,ntype DO m=1,tile_pts(n) l=tile_index(m,n) gs(l) = gs(l) + frac(l,n)*gs_type(l,n) END DO END DO IF ( l_aggregate ) THEN DO n=1,ntype DO m=1,tile_pts(n) l=tile_index(m,n) canhc(l) = canhc(l) + frac(l,n)*ch_type(l,n) vfrac(l) = vfrac(l) + frac(l,n)*vf_type(l,n) DO k=1,nshyd wt_ext(l,k) = wt_ext(l,k) + frac(l,n)*wt_ext_type(l,k,n) END DO END DO END DO DO l=1,land_pts IF ( lake > 0 ) THEN flake(l,1) = frac(l,lake) ELSE flake(l,1) = 0. END IF gs_tile(l,1) = 0. IF (flake(l,1) < 1.) & gs_tile(l,1) = gs(l) / (1. - flake(l,1)) canhc_tile(l,1) = canhc(l) vfrac_tile(l,1) = vfrac(l) DO k=1,nshyd wt_ext_tile(l,k,1) = wt_ext(l,k) END DO END DO ELSE gs_tile(:,:)=0.0 canhc_tile(:,:)=0.0 vfrac_tile(:,:)=0.0 DO n=1,ntype DO m=1,tile_pts(n) l=tile_index(m,n) flake(l,n) = 0. gs_tile(l,n) = gs_type(l,n) canhc_tile(l,n) = ch_type(l,n) vfrac_tile(l,n) = vf_type(l,n) DO k=1,nshyd wt_ext(l,k) = wt_ext(l,k) + frac(l,n)*wt_ext_type(l,k,n) wt_ext_tile(l,k,n) = wt_ext_type(l,k,n) END DO END DO END DO IF ( lake > 0 ) THEN n = lake ! Lake tile DO m=1,tile_pts(n) l=tile_index(m,n) flake(l,n) = 1. END DO END IF END IF DO n=1,npft DO m=1,tile_pts(n) l=tile_index(m,n) gpp(l) = gpp(l) + frac(l,n) * gpp_ft(l,n) npp(l) = npp(l) + frac(l,n) * npp_ft(l,n) resp_p(l) = resp_p(l) + frac(l,n) * resp_p_ft(l,n) END DO END DO !---------------------------------------------------------------------- ! Diagnose the available moisture in the soil profile !---------------------------------------------------------------------- ! Available water for plant transpiration DO n=1,nshyd DO l=1,land_pts smct(l) = smct(l) + MAX( 0. , & wt_ext(l,n)*rho_water*dzsoil(n)* & (sthu(l,n)*v_sat(l,n)-v_wilt(l,n))) END DO END DO ! Add available water for evaporation from bare soil DO l=1,land_pts smct(l) = (1.0-fsoil_tot(l))*smct(l) + & fsoil_tot(l)*rho_water*dzsoil(1)* & MAX(0.0,sthu(l,1))*v_sat(l,1) END DO RETURN END SUBROUTINE physiol ? (smc_sph – called for Sphagnum in place of smc_ext – add this subroutine into the same file in addition to smc_ext or create a new one) ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! SUBROUTINE SMC_SPH----------------------------------------------- ! Description: ! Calculates the soil moisture availability factor and ! the fraction of the transpiration which is extracted from each ! soil layer. ! Documentation : UM Documentation Paper 25 !--------------------------------------------------------------------- ! Subroutine Interface: SUBROUTINE smc_sph (npnts,nshyd,tile_pts,tile_index & , f_root,sthu,v_crit,v_sat,v_wilt & , wt_ext,fsmc,wtp) IMPLICIT NONE ! Subroutine arguments: ! Scalar arguments with intent(IN) : INTEGER & npnts & ! IN Number of gridpoints. ,nshyd & ! IN Number of soil moisture layers. ,tile_pts & ! IN Number of points containing the ! ! given surface type. ,tile_index(npnts) ! IN Indices on the land grid of the ! ! points containing the given ! ! surface type. ! Array arguments with intent(IN) : REAL & f_root(nshyd) & ! IN Fraction of roots in each soil ! ! layer. ,sthu(npnts,nshyd) & ! IN Unfrozen soil moisture content of ! ! each layer as a fraction of ! ! saturation. ,v_crit(npnts,nshyd) & ! IN Volumetric soil moisture ! ! concentration above which ! ! evapotranspiration is not sensitive ! ! to soil water (m3 H2O/m3 soil). ,v_sat(npnts,nshyd) & ! IN Volumetric soil moisture ! ! concentration at saturation ! ! (m3 H2O/m3 soil). ,v_wilt(npnts,nshyd) ! IN Volumetric soil moisture ! ! concentration below which ! ! stomata close (m3 H2O/m3 soil). ! Array arguments with intent(INOUT) : REAL & wt_ext(npnts,nshyd) ! OUT Cummulative fraction of transpiration ! ! extracted from each soil layer ! ! (kg/m2/s). ! Array arguments with intent(OUT) : REAL & fsmc(npnts) ! OUT Soil moisture availability ! ! factor. ! Local scalars: INTEGER & i,j,n ! WORK Loop counters ! Local arrays: REAL & fsmc_l(npnts,nshyd) & ,WTP ! DEV (RC), WORK water table position, m ! WORK Soil moisture availability ! ! factor for each soil layer. !---------------------------------------------------------------------- ! Initialisations !---------------------------------------------------------------------- DO i=1,npnts fsmc(i)=0.0 END DO WTP=0.2 !---------------------------------------------------------------------- ! WTP= <> !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! Sphagnum - Applies equation in right-hand-column for point 1 in Table2 ! from Coppell, Gloor & Holden ! Empirical relation for passive water extraction in sphagnum-peat/soil column !---------------------------------------------------------------------- DO J=1,tile_pts I=tile_index(J) IF (WTP .LE. 0.10 ) THEN fsmc(I)=1.0 ELSE fsmc(I)=EXP(-1*((WTP*100)-10)/16) END IF END DO !---------------------------------------------------------------------- ! Sphagnum - evapotranspiration is very simple here, just based on ! unmodified very shallow rootd fraction, so effectively just ! near-surface evaporation only, since Sphagnum has no active transport !---------------------------------------------------------------------- DO n=1,nshyd DO j=1,tile_pts i=tile_index(j) IF (fsmc(i) > 0.0) & wt_ext(i,n) = f_root(n) END DO END DO RETURN END SUBROUTINE smc_sph ? ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! Purpose: ! Calculates the leaf resistance and net photosynthesis using: ! (i) Collatz et al. (1992) C3 photosynthesis model ! (ii) Jacobs (1994) CI/CA closure. !********************************************************************** SUBROUTINE leaf_limits(land_field,veg_pts,veg_index,ft & , nl,dq,apar,tl,ca,oa,pstar,fsmc & , clos_pts,open_pts,clos_index,open_index & , ci,rd,wcarb,wexpt,wlite) USE c_0_dg_c, ONLY : zerodegc USE pftparm USE surf_param USE c_rmol !added for sphagnum ! Need to work out how to set this ! REAL & ! & FD(NPFT) & ! ! Dark respiration coefficient. ! &,NEFF(NPFT) ! Constant relating VCMAX and leaf N ! DATA FD / 0.015, 0.015, 0.015, 0.025, 0.015 / ! DATA NEFF / 0.8E-3 , 0.8E-3 , 0.8E-3 , 0.4E-3, 0.8E-3 / IMPLICIT NONE ! (mol/sec) / (watts) conversion for PAR: REAL, PARAMETER :: conpar = 2.19e5 INTEGER & land_field & ! IN Total number of land points. ,veg_pts & ! IN Number of vegetated points. ,veg_index(land_field) & ! IN Index of vegetated points ! on the land grid. ,ft ! IN Plant functional type. INTEGER & clos_index(land_field) & ! OUT Index of land points ! with closed stomata. ,clos_pts & ! OUT Number of land points ! with closed stomata. ,open_index(land_field) & ! OUT Index of land points ! with open stomata. ,open_pts ! OUT Number of land points ! ! with open stomata. REAL & nl(land_field) & ! IN Leaf nitrogen ! ! concentration (kg N/kg C). ,dq(land_field) & ! IN Canopy level specific humidity ! ! deficit (kg H2O/kg air). ,apar(land_field) & ! IN Absorbed PAR (W/m2) ,tl(land_field) & ! IN Leaf temperature (K). ,ca(land_field) & ! IN Canopy CO2 pressure (Pa). ,oa(land_field) & ! IN Atmospheric O2 pressure (Pa). ,pstar(land_field) & ! IN Atmospheric pressure (Pa). ,fsmc(land_field) & ! IN Soil water factor. ,ci(land_field) & ! OUT Internal CO2 pressure (Pa). ,rd(land_field) & ! OUT Dark respiration (mol CO2/m2/s). ,wcarb(land_field) & ! OUT Carboxylation, ... ,wlite(land_field) & ! ... Light, and ... ,wexpt(land_field) & ! ... export limited gross ... ! ! ... photosynthetic rates ... ! ! ... (mol CO2/m2/s). ,ccarbi(land_field) & ! WORK Internal sphagnum CO2 pressure limited by Carboxylation (Pa). ,clitei(land_field) & ! WORK Internal sphagnum CO2 pressure limited by Light(Pa).! ,acr(land_field) & ! WORK Absorbed PAR ! ! (mol photons/m2/s). ,ccp(land_field) & ! WORK Photorespiratory compensatory ! ! point (mol/m3). ,denom(land_field) & ! WORK Denominator in equation for VCM ,kc(land_field) & ! WORK Michaelis constant for CO2 (Pa) ,ko(land_field) & ! WORK Michaelis constant for O2 (Pa). ,qtenf(land_field) & ! WORK Q10 function. ,tau(land_field) & ! WORK CO2/O2 specificity ratio. ,tdegc(land_field) & ! WORK Leaf temperature (deg C). ,vcm(land_field) & ! WORK Maximum rate of carboxylation ! ! of Rubisco (mol CO2/m2/s). ,vcmax(land_field) & ! WORK Maximum rate of carboxylation ! ! of Rubisco - without the ! ! temperature factor ! ! (mol CO2/m2/s). ,conv(land_field) & ! WORK Factor for converting mol/m3 ! ! into Pa (J/mol). !,VMMAX & ! !maximum Vmax for mosses ,gl(land_field) & ! WORK Leaf conductance for H2O (m/s). ,anetl(land_field) & ! WORK Net leaf photosynthesis ! (mol CO2/m2/s/LAI). ,alpha_rc & ! WORK - sphag ,K_rc & ! WORK - sphag ,gamma_rc & ! WORK - sphag ,Ac & ! WORK - sphag ,Bc & ! WORK - sphag ,Cc & ! WORK - sphag ,Al & ! WORK - sphag ,Bl & ! WORK - sphag ,Cl & ! WORK - sphag ,glco2(land_field) INTEGER & j,l ! WORK Loop counters. CHARACTER*50 char_debug !---------------------------------------------------------------------- ! Initialise counters !---------------------------------------------------------------------- clos_pts = 0 open_pts = 0 DO j=1,veg_pts l = veg_index(j) !---------------------------------------------------------------------- ! Calculate the points with closed stomata !---------------------------------------------------------------------- ! all stom points OPEN for sphag IF ( & ((pftName(ft).NE.'sphag') & .AND. & ( fsmc(l)==0.0 .OR. dq(l) >= dqcrit(ft) & .OR. apar(l)==0.0 ) ) THEN clos_pts = clos_pts + 1 clos_index(clos_pts) = j ELSE open_pts = open_pts + 1 open_index(open_pts) = j END IF END DO !---------------------------------------------------------------------- ! Calculate the photosynthetic parameters !---------------------------------------------------------------------- !DIR$ IVDEP DO j=1,open_pts l = veg_index(open_index(j)) vcmax(l) = neff(ft) * nl(l) tdegc(l) = tl(l) - zerodegc ! TAU is the Rubisco specificity for CO2 relative to O2. The numbers ! in this equation are from Cox, HCTN 24, "Description ... Vegetation ! Model", equation 53. tau(l) = 2600.0 * (0.57 ** (0.1 * (tdegc(l) - 25.0))) ccp(l) = 0.5 * oa(l) / tau(l) * REAL(c3(ft)) qtenf(l) = vcmax(l) * (q10_leaf ** (0.1 * (tdegc(l) - 25.0))) denom(l) = (1 + EXP (0.3 * (tdegc(l) - tupp(ft)))) & * (1 + EXP (0.3 * (tlow(ft) - tdegc(l)))) vcm(l) = qtenf(l) / denom(l) ! Cox, HCTN 24, equation 49. rd(l) = fd(ft) * vcm(l) !---------------------------------------------------------------------- ! Calculate the factor for converting mol/m3 into Pa (J/m3). !---------------------------------------------------------------------- conv(l) = rmol * tl(l) !---------------------------------------------------------------------- ! Calculate the internal CO2 pressure (Jacobs, 1994). !---------------------------------------------------------------------- IF (pftName(ft).EQ.'sphag') THEN glco2(l) = 0.0237*0.07 ! simple fixed value for sphagnum END IF gl(l) = ratio * glco2(l) !ratio=1.6 anyway ! then calc simultanoeus equ from all these and then do minimum as per usual below... alpha_rc = gl(l) / (1.6 *rmol* tl(l)* fsmc(l)) K_rc = kc(l)*(1+(oa(l)/ko(l))) gamma_rc= alpha(ft) * acr(l) Ac = alpha_rc/vcm(l) Bc = ((alpha_rc/vcm(l))*(K_rc-ca(l)))-rd(l)+1 Cc = -1*( (K_rc*ca(l)*alpha_rc/vcm(l)) + (rd(l)*K_rc) + ccp(l) ) Al= alpha_rc/gamma_rc Bl= ((alpha_rc/gamma_rc)*((2*ccp(l))-ca(l)))-rd(l)+1 Cl= -1*ccp(l)*( (2*ca(l)*alpha_rc/gamma_rc ) + ( 2*rd(l) ) + 1 ) ccarbi(l) = ( -Bc+sqrt(Bc**2-(4*Ac*Cc)) )/(2*Ac) ! WORK Internal sphagnum CO2 pressure limited by Carboxylation (Pa). clitei(l) = ( -Bl+sqrt(Bl**2-(4*Al*Cl)) )/(2*Al) ! WORK Internal sphagnum CO2 pressure limited by Light(Pa). ELSE !non-sphag ci(l) = (ca(l) - ccp(l)) * f0(ft) & * (1 - dq(l) / dqcrit(ft)) + ccp(l) ENDIF !---------------------------------------------------------------------- ! Convert absorbed PAR into mol PAR photons/m2/s !---------------------------------------------------------------------- acr(l) = apar(l) / conpar END DO !! REPEAT FOR CLOSED STOMATA (not sphagnum) !DIR$ IVDEP DO j=1,clos_pts l = veg_index(clos_index(j)) !---------------------------------------------------------------------- ! Calculate the factor for converting mol/m3 into Pa (J/m3). !---------------------------------------------------------------------- conv(l) = rmol * tl(l) vcmax(l) = neff(ft) * nl(l) tdegc(l) = tl(l) - zerodegc ! TAU is the Rubisco specificity for CO2 relative to O2. The numbers ! in this equation are from Cox, HCTN 24, "Description ... Vegetation ! Model", equation 53. tau(l) = 2600.0 * (0.57 ** (0.1 * (tdegc(l) - 25.0))) ccp(l) = 0.5 * oa(l) / tau(l) * REAL( c3(ft) ) qtenf(l) = vcmax(l) * (q10_leaf ** (0.1 * (tdegc(l) - 25.0))) denom(l) = (1 + EXP (0.3 * (tdegc(l) - tupp(ft)))) & * (1 + EXP (0.3 * (tlow(ft) - tdegc(l)))) vcm(l) = qtenf(l) / denom(l) ! Cox, HCTN 24, equation 49. rd(l) = fd(ft) * vcm(l) !---------------------------------------------------------------------- ! Calculate the internal CO2 pressure (Jacobs, 1994). !---------------------------------------------------------------------- ci(l) = (ca(l) - ccp(l)) * f0(ft) & * (1 - dq(l) / dqcrit(ft)) + ccp(l) END DO !---------------------------------------------------------------------- ! Calculate the gross photosynthesis for RuBP-Carboxylase, Light and ! Export limited photosynthesis (Collatz et al., 1992). !---------------------------------------------------------------------- !DIR$ IVDEP IF (c3(ft) == 1) THEN DO j=1,open_pts l = veg_index(open_index(j)) ! The numbers ! in these 2 equations are from Cox, HCTN 24, "Description ... Vegetation ! Model", equations 54 and 55. kc(l) = 30.0 * (2.1 ** (0.1 * (tdegc(l) - 25.0))) ko(l) = 30000.0 * (1.2 ** (0.1 * (tdegc(l) - 25.0))) IF (pftName(ft).EQ.'sphag') & THEN wcarb(l) = vcm(l) * (ccarbi(l) - ccp(l)) & / (ccarbi(l) + kc(l) * (1. + oa(l) / ko(l))) wlite(l) = alpha(ft) * acr(l) * (clitei(l) - ccp(l)) & / (clitei(l) + 2 * ccp(l)) wlite(l) = MAX(wlite(l), TINY(1.0e0)) wexpt(l) = fwe_c3 * vcm(l) ELSE wcarb(l) = vcm(l) * (ci(l) - ccp(l)) & / (ci(l) + kc(l) * (1. + oa(l) / ko(l))) wlite(l) = alpha(ft) * acr(l) * (ci(l) - ccp(l)) & / (ci(l) + 2 * ccp(l)) wlite(l) = MAX(wlite(l), TINY(1.0e0)) wexpt(l) = fwe_c3 * vcm(l) END IF END DO ELSE DO j=1,open_pts l = veg_index(open_index(j)) wcarb(l) = vcm(l) wlite(l) = alpha(ft) * acr(l) wlite(l) = MAX(wlite(l), TINY(1.0e0)) wexpt(l) = fwe_c4 * vcm(l) * ci(l) / pstar(l) END DO END IF RETURN END SUBROUTINE leaf_limits ? ? ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! Purpose: ! Subroutine to calculate gridbox mean values of surface conductance ! and carbon fluxes. Also returns net primary productivity, leaf ! turnover and wood respiration of each plant functional type for ! driving TRIFFID. ! Calculates the leaf resistance and net photosynthesis using: ! (i) Collatz et al. (1992) C3 photosynthesis model, and the ! Collatz et al. (1991) C4 photosynthesis model. ! (ii) Jacobs (1994) CI/CA closure. ! Written by Peter Cox (February 1996) ! Adapted for MOSES II tile model by Richard Essery (July 1997) ! Coded into UMVN 6.2 by Pete Falloon (July 2006) !********************************************************************** SUBROUTINE leaf (land_field,veg_index,ft & , clos_pts,open_pts,clos_index,open_index & , fsmc,tl,ca,ci,rd,wcarb,wexpt,wlite & , gl,al) USE pftparm USE c_rmol USE surf_param IMPLICIT NONE INTEGER & land_field & ! IN Total number of land points. ,veg_index(land_field) & ! IN Index of vegetated points ! on the land grid. ,ft & ! IN Plant functional type. ,clos_index(land_field) & ! IN Index of land points ! with closed stomata. ,clos_pts & ! IN Number of land points ! with closed stomata. ,open_index(land_field) & ! IN Index of land points ! with open stomata. ,open_pts ! IN Number of land points ! with open stomata. REAL & fsmc(land_field) & ! IN Soil water factor. ,tl(land_field) & ! IN Leaf temperature (K). ,rd(land_field) & ! IN Dark respiration (mol CO2/m2/s). ,ca(land_field) & ! IN Canopy CO2 pressure (Pa). ,ci(land_field) & ! IN Internal CO2 pressure (Pa). ,wcarb(land_field) & ,wlite(land_field) & ,wexpt(land_field) & ! IN Carboxylation, ! Light, and ! export limited gross ! photosynthetic rates ! (mol CO2/m2/s). ,gl(land_field) & ! OUT Leaf conductance for H2O (m/s). ,al(land_field) & ! OUT Net Leaf photosynthesis ! ! (mol CO2/m2/s). ,b1(land_field) & ,b2(land_field) & ,b3(land_field) & ! WORK Coefficients of the quadratic. ,conv(land_field) & ! WORK Factor for converting mol/m3 ! ! into Pa (J/mol). ,glco2(land_field) & ! WORK Leaf conductance for CO2 (m/s). ,wl(land_field) & ! WORK Gross leaf photosynthesis ! ! (mol CO2/m2/s). ,wp(land_field) & ! WORK Smoothed minimum of ! ! Carboxylation and Light ! ! limited gross photosynthesis ! ! (mol CO2/m2/s). ,Vm & ! molar volume = RT/P INTEGER & j,l ! WORK Loop counters. CHARACTER*50 char_debug !---------------------------------------------------------------------- ! Calculate the co-limited rate of gross photosynthesis !---------------------------------------------------------------------- !DIR$ IVDEP DO j=1,open_pts l = veg_index(open_index(j)) b1(l) = beta1 b2(l) = - (wcarb(l) + wlite(l)) b3(l) = wcarb(l) * wlite(l) wp(l) = -b2(l)/(2*b1(l)) & - SQRT(b2(l)*b2(l)/(4*b1(l)*b1(l)) - b3(l)/b1(l)) END DO !DIR$ IVDEP DO j=1,open_pts l = veg_index(open_index(j)) b1(l) = beta2 b2(l) = - (wp(l) + wexpt(l)) b3(l) = wp(l) * wexpt(l) wl(l) = -b2(l)/(2*b1(l)) & - SQRT(b2(l)*b2(l)/(4*b1(l)*b1(l)) - b3(l)/b1(l)) END DO !---------------------------------------------------------------------- ! Carry out calculations for points with open stomata !---------------------------------------------------------------------- !DIR$ IVDEP DO j=1,open_pts l = veg_index(open_index(j)) !---------------------------------------------------------------------- ! Calculate the net rate of photosynthesis !---------------------------------------------------------------------- al(l) = (wl(l) - rd(l)) * fsmc(l) !---------------------------------------------------------------------- ! Calculate the factor for converting mol/m3 into Pa (J/m3). !---------------------------------------------------------------------- conv(l) = rmol * tl(l) !---------------------------------------------------------------------- ! Diagnose the leaf conductance !---------------------------------------------------------------------- IF (pftName(ft).EQ.'sphag') THEN glco2(l) = 0.07*0.0237 !same as calc also used in leaf_limits subroutine ! FIXED VALUE GLCO2 ELSE ! not sphagnum glco2(l) = (al(l) * conv(l)) / (ca(l) - ci(l)) END IF gl(l) = ratio * glco2(l) END DO IF (pftName(ft).NE.'sphag') THEN !---------------------------------------------------------------------- ! Close stomata at points with negative or zero net photosynthesis ! or where the leaf resistance exceeds its maximum value. !---------------------------------------------------------------------- !DIR$ IVDEP DO j=1,open_pts l = veg_index(open_index(j)) IF (gl(l).le.glmin(ft) .OR. al(l).le.0.0) THEN gl(l) = glmin(ft) glco2(l) = gl(l) / ratio al(l) = -rd(l)*fsmc(l) END IF END DO !---------------------------------------------------------------------- ! Define fluxes and conductances for points with closed stomata !---------------------------------------------------------------------- !DIR$ IVDEP DO j=1,clos_pts l = veg_index(clos_index(j)) gl(l) = glmin(ft) glco2(l) = gl(l) / ratio al(l) = -rd(l) *fsmc(l) END DO END IF RETURN END SUBROUTINE leaf