== Shushi page == Bugs in MERGE-OCN (r2558) 1. corrected in slowproc.f90 line 4972, soiltext ---> soilbulk line 4975, soiltext ---> soilph 2. corrected bug as ticket#109 3. corrected bug for "carbon"initialization, actually, this variable should contain a element dimension in the calculation, at least icarbon and initrogen. So it is better to rename this variable as ele_pools... in stomate_io.f90 {{{ 232 !!!!pssdebugok#1+++ 233 REAL(r_std), DIMENSION(npts,ncarb,nvm,nelements), INTENT(out) :: carbon !! som pool: active, slow, or passive, surface 234 !! for carbon and nitrogen (gC/m**2) 235 !!!!pssdebugok#1--- 236 ! REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(out) :: carbon !! carbon pool: active, slow, or passive, 237 ! !! (gC/m**2) }}} 4. leaching only have iammonium and initrate, other gas (inox, initrous) for leaching should be zero; in stomate_soilcarbon.f90 debug: initiate leaching(:,:,:) = 0 {{{ 654 !!!!pssdebug++++ 655 leaching(:,:,:) = 0 656 !!!!pssdebug---- }}} For denitrification calculation, there is no leaching(:,m,inox) and leaching(:,m,initrous) calculation before. I think we do not need to minus leaching for this two species here. {{{ 1079 !!!!pssdebugok+ ! 1080 ! NO consumption 1081 denitrification(:,m,i_nh4_to_no) = MIN(soil_n_min(:,m,inox), & 1082 anvf(:,m) * ( mu_no(:) / Y_no + M_no * soil_n_min(:,m,inox) / sum_n(:)) * bact(:,m) * & 1083 24. * dt ) 1084 ! 1085 ! N2O consumption 1086 denitrification(:,m,i_nh4_to_n2o) = MIN(soil_n_min(:,m,initrous), & 1087 anvf(:,m) * ( mu_n2o(:) / Y_n2o + M_n2o * soil_n_min(:,m,initrous) / sum_n(:)) * bact(:,m) * & 1088 24. * dt ) 1089 !!!!pssdebugok- 1090 !!!!pssdebug ! NO consumption 1091 !!!!pssdebug denitrification(:,m,i_nh4_to_no) = MIN(soil_n_min(:,m,inox)-leaching(:,m,inox), & 1092 !!!!pssdebug anvf(:,m) * ( mu_no(:) / Y_no + M_no * soil_n_min(:,m,inox) / sum_n(:)) * bact(:,m) * & 1093 !!!!pssdebug 24. * dt ) 1094 !!!!pssdebug ! 1095 !!!!pssdebug ! N2O consumption 1096 !!!!pssdebug denitrification(:,m,i_nh4_to_n2o) = MIN(soil_n_min(:,m,initrous)-leaching(:,m,initrous), & 1097 !!!!pssdebug anvf(:,m) * ( mu_n2o(:) / Y_n2o + M_n2o * soil_n_min(:,m,initrous) / sum_n(:)) * bact(:,m) * & 1098 !!!!pssdebug 24. * dt ) }}} 5. bug corrected: update nitrogen pool of circ_class_biomass when carbon pool is updated, otherwise CN ratio of each pool will increase with timestep in src_stomate/stomate_growth_fun_all.f90 {{{ 2049 !!!!pssdebug_ok add circ_class_biomass initrogen to keep CN ratio, but need to check nitrogen mass balance later!!!! 2050 !!!!CAUTION1: nitgrogen maybe imbalance, should substracb from ilibile N pool 2051 !!!!CAUTION2: fruit CN ratio is assumed as same as root, this may be modified in the later version 2052 circ_class_biomass(ipts,j,:,ileaf,initrogen) = circ_class_biomass(ipts,j,:,ileaf,initrogen) + & 2053 ( Cl_inc(:) * total_inc / circ_class_n(ipts,j,:) )/cn_leaf(ipts,j) 2054 circ_class_biomass(ipts,j,:,isapabove,initrogen) = circ_class_biomass(ipts,j,:,isapabove,initrogen) + & 2055 ( Cs_inc(:) * alloc_sap_above * total_inc / circ_class_n(ipts,j,:) )/cn_leaf(ipts,j)*fcn_wood(j) 2056 circ_class_biomass(ipts,j,:,isapbelow,initrogen) = circ_class_biomass(ipts,j,:,isapbelow,initrogen) + & 2057 ( Cs_inc(:) * (un - alloc_sap_above) * total_inc / circ_class_n(ipts,j,:) )/cn_leaf(ipts,j)*fcn_wood(j) 2058 circ_class_biomass(ipts,j,:,iroot,initrogen) = circ_class_biomass(ipts,j,:,iroot,initrogen) + & 2059 ( Cr_inc(:) * total_inc / circ_class_n(ipts,j,:) )/cn_leaf(ipts,j)*fcn_root(j) 2060 circ_class_biomass(ipts,j,:,ifruit,initrogen) = circ_class_biomass(ipts,j,:,ifruit,initrogen) + & 2061 ( Cf_inc(:) * total_inc / circ_class_n(ipts,j,:) )/cn_leaf(ipts,j)*fcn_root(j) 2062 !!!!pssdebug_ok- }}} 6. nitrogen demand for new biomass problem, I think in the code, the N demand for the new biomass is not equal to bm_alloc. Two problems: 1) the growth_resp is already subtracted before, and do not need to add (1-frac_growthresp(j)) to calculate N demand; 2) bm_alloc_tot should be equal to sum of bm_alloc = bm_alloc_tot(ipts,j) - residual(ipts,j). Here I think it is better to use the demand as the total new biomass allocated in one timestep. >>DSG: I can confirm, (1-frac_growthresp(j)) must not be subtracted. in src_stomate/stomate_growth_fun_all.f90 {{{ 2857 ! only check if there is biomass growth 2858 IF ( costf.GT.min_stomate ) THEN 2859 2860 ! fraction of labile N allocatable for growth 2861 IF ( control%ok_functional_allocation) THEN 2862 !no growth respiration calculated here! 2863 n_avail = & !MIN( & 2864 MAX(biomass(ipts,j,ilabile,initrogen)*0.9-min_stomate,0.0)!,& 2865 ELSE 2866 n_avail = n_alloc_tot(ipts,j) 2867 ENDIF 2868 2869 ! carbon growth possible given nitrogen availability and current nitrogen concentration 2870 ! bm_supply_n = n_avail / costf / (1.-frac_growthresp(j)) * cn_leaf(ipts,j) 2871 !!!!pssdebug, I think here we DO NOT need to substract growth_resp in this block, because growth_resp is already substracted before 2872 !(1.-frac_growthresp(j)) 2873 bm_supply_n = n_avail / costf * cn_leaf(ipts,j) 2874 !!!!pssdebug---- 2875 ! elasticity of leaf nitrogen concentration 2876 deltacnmax=exp(-(1./(cn_leaf(ipts,j)*0.5*(1./cn_leaf_max(j)+1./cn_leaf_min(j))))**8) 2877 2878 IF ( bm_alloc_tot(ipts,j) .GT. bm_supply_n ) THEN 2879 ! case of not enough nitrogen to sustain intended growth, 2880 ! reduce carbon allocation to meet nitrogen availability 2881 ! taking account of the maximal change of nitrogen concentrations 2882 2883 ! delta of nitrogen concetrations in response to nitrogen deficit 2884 IF( (control%ok_functional_allocation.AND.lab_fac(ipts,j).GT.0.3) .OR. & 2885 (.NOT.control%ok_functional_allocation.AND.use_reserve(ipts,j).GT.min_stomate))THEN 2886 deltacn=1.0 2887 ELSE 2888 deltacnmax=Dmax * (1.-deltacnmax) 2889 ! deltacn = n_avail / ( bm_alloc_tot(ipts,j) * (1.-frac_growthresp(j)) * costf * 1./cn_leaf(ipts,j) ) 2890 !!!!pssdebug++++ 2891 deltacn = n_avail / ( bm_alloc_tot(ipts,j) * costf * 1./cn_leaf(ipts,j) ) 2892 !!!!pssdebug---- 2893 deltacn=MIN(MAX(deltacn,1.0-deltacnmax),1.0) 2894 ENDIF 2895 2896 ! nitrogen demand given possible nitrogen concentration change 2897 n_alloc_tot(ipts,j) = MIN( n_avail , & 2898 ! bm_alloc_tot(ipts,j) * (1.-frac_growthresp(j)) * costf * 1./cn_leaf(ipts,j)*deltacn ) 2899 !!!!pssdebug++++ 2900 (bm_alloc_tot(ipts,j) - residual(ipts,j)) * costf * 1./cn_leaf(ipts,j)*deltacn ) 2901 !!!!pssdebug---- 2902 ! if not succesful, reduce growth 2903 ! constrain carbon used for growth dependent on available nitrogen 2904 ! under the assumption that the f_allocs are piecewise linear with 2905 ! bm_alloc_tot, which is first-order correct 2906 bm_alloc_tot(ipts,j) = MIN( bm_alloc_tot(ipts,j) , & 2907 ! n_alloc_tot(ipts,j) / costf / (1.-frac_growthresp(j)) / (1./cn_leaf(ipts,j) * deltacn) ) 2908 !!!!pssdebug++++ 2909 n_alloc_tot(ipts,j) / costf / (1./cn_leaf(ipts,j) * deltacn) ) 2910 !!!!pssdebug---- 2911 2912 ELSE 2913 ! sufficient nitrogen, increaase of leaf nitrogen concentration dependent on 2914 ! distance to maximal leaf nitrogen concentration 2915 ! cannot change leaf C:N in bud burst period 2916 ! nitrogen constrained such that nitrogen concentration can only increase 1% per day 2917 IF( (control%ok_functional_allocation.AND.lab_fac(ipts,j).GT.0.3) .OR. & 2918 (.NOT.control%ok_functional_allocation.AND.use_reserve(ipts,j).GT.min_stomate))THEN 2919 deltacn=1.0 2920 ELSE 2921 deltacnmax=Dmax * deltacnmax !(1.-deltacnmax) 2922 deltacn = n_avail / & 2923 !!!!pssdebug ( bm_alloc_tot(ipts,j) * (1.-frac_growthresp(j)) * costf * 1./cn_leaf(ipts,j) ) 2924 ( bm_alloc_tot(ipts,j) * costf * 1./cn_leaf(ipts,j) ) 2925 deltacn=MIN(MAX(deltacn,1.0),1.+deltacnmax) 2926 ENDIF 2927 n_alloc_tot(ipts,j) = MIN( n_avail , & 2928 !!!!pssdebug bm_alloc_tot(ipts,j) * (1.-frac_growthresp(j)) * costf * 1./cn_leaf(ipts,j)*deltacn ) 2929 (bm_alloc_tot(ipts,j) - residual(ipts,j)) * costf * 1./cn_leaf(ipts,j)*deltacn ) 2930 2931 ENDIF 2932 2933 ENDIF }}} 7. bug corrected: avoid negative values for FixNH4 in src_stomate/stomate_soilcarbon.f90 {{{ 838 WHERE(soil_n_min(:,m,iammonium).GT.min_stomate) 839 FixNH4(:) = a_FixNH4 + b_FixNH4 * & 840 log10(soil_n_min(:,m,iammonium) / ( dpu_max * Bd(:) ) ) & 841 * MIN(clay(:),clay_max) / clay_max 842 ELSEWHERE 843 FixNH4(:) = 0.0 844 ENDWHERE 845 FixNH4(:) = MAX( zero, FixNH4(:)) }}} 8. The unit of history write from nitrogen_dynamics subroutine is per timestep, not per day 9. I plot parameter sensitivity tests for most of the equations in nitrogen_dynamics subroutine. It seems OK, but need further careful test in the future.