diff --git a/couplage/CARAIB/ver01_Iv_couplage/carbon_npp_cal.f b/couplage/CARAIB/ver01_Iv_couplage/carbon_npp_cal.f
index ade0cadee5e30e2ff6f13b511524fd29e7371a06..ee0746a0fa529b220676fb84d4ffeeae5fe8ade8 100644
--- a/couplage/CARAIB/ver01_Iv_couplage/carbon_npp_cal.f
+++ b/couplage/CARAIB/ver01_Iv_couplage/carbon_npp_cal.f
@@ -85,7 +85,7 @@ c
integer , intent(in) :: ngt,ny0,iread,itmt
integer ih,ipool,iday2,iday_sow,iharv
real*4 , intent(in) :: tmt_it
- real*4 faparl,tot_frac,faparh
+ real*4 faparl,tot_frac,faparh,rbmax
c------------------------------------------------
c second argument of year_iteration :
@@ -207,33 +207,33 @@ c=======================================================================
c yearly gpp and fractionation calculation
c=======================================================================
do ip = 1,npft
- iharv = 0
- if(frac(ip).gt.1.e-7) then
- do iday = 1, nd
+ iharv = 0
+ if(frac(ip).gt.1.e-7) then
+ do iday = 1, nd
if(ip.le.npft0.or.ip.ge.npft0+ncrop+1) then
- ygppf(ip) = ygppf(ip) + xgpp(ip,iday)
- ygppf_ini(ip,ngt) = ygppf(ip)
- else
- if(iharv.eq.0) then
- ygppf(ip) = ygppf(ip) + xgpp(ip,iday)
- if(maturity(ip).eq.-999) ygppf(ip) =0.
- ygppf_ini(ip,ngt) = ygppf(ip)
- if(iphase(ip,iday).eq.2) then
- ygppf_ini(ip,ngt) = 0.
- iharv = 1
- endif
- else
- ygppf_ini(ip,ngt) = ygppf_ini(ip,ngt) + xgpp(ip,iday)
- endif
- endif
+ ygppf(ip) = ygppf(ip) + xgpp(ip,iday)
+ ygppf_ini(ip,ngt) = ygppf(ip)
+ else
+ if(iharv.eq.0) then
+ ygppf(ip) = ygppf(ip) + xgpp(ip,iday)
+ if(maturity(ip).eq.-999) ygppf(ip) =0.
+ ygppf_ini(ip,ngt) = ygppf(ip)
+ if(iphase(ip,iday).eq.2) then
+ ygppf_ini(ip,ngt) = 0.
+ iharv = 1
+ endif
+ else
+ ygppf_ini(ip,ngt) = ygppf_ini(ip,ngt) + xgpp(ip,iday)
+ endif
+ endif
yfractf(ip) = yfractf(ip) + xfract(ip,iday)
- enddo
- if (ygppf(ip).eq.0.) then
+ enddo
+ if (ygppf(ip).eq.0.) then
yfractf(ip) = 0.
- else
+ else
yfractf(ip) = yfractf(ip) / ygppf(ip)
- endif
- endif
+ endif
+ endif
end do
@@ -262,62 +262,67 @@ c grid point gpp, npp, lai, biomass, 13C fract and fire emission
c=======================================================================
do iday = 1, nd
- iday2 = iday-1
- if(iday .eq. 1) iday2 = nd
- do ip = 1, npft
- xdgpp(iday) = xdgpp(iday) + frac(ip) * xgpp(ip,iday)
- xdnpp(iday) = xdnpp(iday) + frac(ip) * xnpp(ip,iday)
- xdlai(iday) = xdlai(iday) + frac(ip) * xlai(ip,iday)
- xdfract(iday) = xdfract(iday)
+ iday2 = iday-1
+ if(iday .eq. 1) iday2 = nd
+ do ip = 1, npft
+ xdgpp(iday) = xdgpp(iday) + frac(ip) * xgpp(ip,iday)
+ xdnpp(iday) = xdnpp(iday) + frac(ip) * xnpp(ip,iday)
+ xdlai(iday) = xdlai(iday) + frac(ip) * xlai(ip,iday)
+ xdfract(iday) = xdfract(iday)
& + frac(ip) * xfract(ip,iday)
- do ipool = 1, npool
- xdbiom(iday) = xdbiom(iday) + frac(ip)*carbon(ip,ipool,iday)
- ybiomf(ip) = ybiomf(ip) + carbon(ip,ipool,iday)/float(nd)
- xbiom(ip,iday) = xbiom(ip,iday) + carbon(ip,ipool,iday)
- xemifire(ip,iday) = xemifire(ip,iday)
+ do ipool = 1, npool
+ xdbiom(iday) = xdbiom(iday) + frac(ip)*carbon(ip,ipool,iday)
+ ybiomf(ip) = ybiomf(ip) + carbon(ip,ipool,iday)/float(nd)
+ xbiom(ip,iday) = xbiom(ip,iday) + carbon(ip,ipool,iday)
+ xemifire(ip,iday) = xemifire(ip,iday)
& + emi_burn_veg(ip,ipool,iday)
- xcharvest(ip,iday) = xcharvest(ip,iday)
+ xcharvest(ip,iday) = xcharvest(ip,iday)
& + charvest(ip,ipool,iday)
- yemifiref(ip) = yemifiref(ip) + emi_burn_veg(ip,ipool,iday)
- enddo
-
- if(ip.le.npft0.or.ip.ge.npft0+ncrop+1) then
- ybiomag(ip) = ybiomag(ip)
+ yemifiref(ip) = yemifiref(ip) + emi_burn_veg(ip,ipool,iday)
+ enddo
+
+ if(ip.le.npft0.or.ip.ge.npft0+ncrop+1) then
+ ybiomag(ip) = ybiomag(ip)
& + (carbon(ip,1,iday)
- & + carbon(ip,2,iday)-root_biomass(ip,iday))/float(nd)
- ybiombg(ip) = ybiombg(ip)+root_biomass(ip,iday)/float(nd)
+ & + max(0.,carbon(ip,2,iday)-root_biomass(ip,iday)))
+ & /float(nd)
+ ybiombg(ip) = ybiombg(ip)+min(carbon(ip,2,iday),
+ & root_biomass(ip,iday))/float(nd)
c Ingrid
- ybiomtot(ip) = carbon(ip,1,iday)+carbon(ip,2,iday)
- else
- if (iphase(ip,iday2).eq.1.and.iphase(ip,iday).eq.2) then
- if(iday.eq.1) then
- root_biomass(ip,iday2) = root_ini(ip)
- carbon(ip,1,iday2) = ycar_ini(ip,1,ngt) !carbon_ini(ip,1)
- carbon(ip,2,iday2) = ycar_ini(ip,2,ngt) !carbon_ini(ip,2)
- endif
- ybiomag(ip) = carbon(ip,1,iday2)
- & + carbon(ip,2,iday2)-root_biomass(ip,iday2)
- ybiombg(ip) = root_biomass(ip,iday2)
+ ybiomtot(ip) = carbon(ip,1,iday)+carbon(ip,2,iday)
+ else
+ if (iphase(ip,iday2).eq.1.and.iphase(ip,iday).eq.2) then
+ if(iday.eq.1) then
+ carbon(ip,1,iday2) = ycar_ini(ip,1,ngt) !carbon_ini(ip,1)
+ carbon(ip,2,iday2) = ycar_ini(ip,2,ngt) !carbon_ini(ip,2)
+ rbmax = (rootsh(ip)/(1.+rootsh(ip)))
+ & *(carbon(ip,1,iday2)+carbon(ip,2,iday2))
+ root_ini(ip) = min(carbon(ip,2,iday2),rbmax)
+ root_biomass(ip,iday2) = root_ini(ip)
+ endif
+ ybiomag(ip) = carbon(ip,1,iday2)
+ & + max(0.,carbon(ip,2,iday2)-root_biomass(ip,iday2))
+ ybiombg(ip) = max(0.,root_biomass(ip,iday2))
c Ingrid
- ybiomtot(ip) = carbon(ip,1,iday2) + carbon(ip,2,iday2)
+ ybiomtot(ip) = carbon(ip,1,iday2) + carbon(ip,2,iday2)
endif
- endif
+ endif
c Ingrid
- yield(ip)=yield_fac(ip)*ybiomtot(ip)
+ yield(ip)=yield_fac(ip)*ybiomtot(ip)
- enddo
+ enddo
- if(xdgpp(iday).ne.0.) then
- xdfract(iday) = xdfract(iday) / xdgpp(iday)
- else
- if(abs(xdfract(iday)).gt.1.e-20) then
+ if(xdgpp(iday).ne.0.) then
+ xdfract(iday) = xdfract(iday) / xdgpp(iday)
+ else
+ if(abs(xdfract(iday)).gt.1.e-20) then
write(99961,*) 'error, fractionation:',iday,xdfract(iday),
& xdgpp(iday)
- endif
- xdfract(iday) = 0.
- endif
+ endif
+ xdfract(iday) = 0.
+ endif
enddo
@@ -327,7 +332,7 @@ c=======================================================================
tot_frac = 0.
do ip = 1, npft
- tot_frac = tot_frac + frac(ip)
+ tot_frac = tot_frac + frac(ip)
end do
c write(*,*)'npp_cal:after loop iday,ip, ybiom'
@@ -336,35 +341,35 @@ c fapar calculation
c=======================================================================
do iday = 1, nd
- xdfapar(iday) = 0.
- if(partoc(iday).gt.0.) then
- do ip = 1, npft
- faparl = 0.
- if(frac(ip) .gt. 1.e-7) then
- do ih = 1, nh2
- faparh = 0.
-
- faparh = faparh
+ xdfapar(iday) = 0.
+ if(partoc(iday).gt.0.) then
+ do ip = 1, npft
+ faparl = 0.
+ if(frac(ip) .gt. 1.e-7) then
+ do ih = 1, nh2
+ faparh = 0.
+
+ faparh = faparh
& + (sunhour(iday) * (fsun(ip,iday,ih)*apar(ip,iday,ih,1)
& + (1.-fsun(ip,iday,ih))*apar(ip,iday,ih,2))
& + (1.-sunhour(iday)) * apar(ip,iday,ih,3))
& *corr_hour(iday,ih)/radun
- faparl = faparl + faparh
- enddo
- faparl = faparl * xlai(ip,iday) * frac(ip)
- endif
- xdfapar(iday) = xdfapar(iday) + faparl
- enddo
- xdfapar(iday) = xdfapar(iday) / partoc(iday)
-
- if(tot_frac.gt.0.0) then
- xdfapar(iday)= xdfapar(iday) / tot_frac
- endif
- endif
- if(xdfapar(iday).lt.0. .or. xdfapar(iday).gt.1.) then
- write(99961,*) 'error L',iday,xdfapar(iday)
- endif
+ faparl = faparl + faparh
+ enddo
+ faparl = faparl * xlai(ip,iday) * frac(ip)
+ endif
+ xdfapar(iday) = xdfapar(iday) + faparl
+ enddo
+ xdfapar(iday) = xdfapar(iday) / partoc(iday)
+
+ if(tot_frac.gt.0.0) then
+ xdfapar(iday)= xdfapar(iday) / tot_frac
+ endif
+ endif
+ if(xdfapar(iday).lt.0. .or. xdfapar(iday).gt.1.) then
+ write(99961,*) 'error L',iday,xdfapar(iday)
+ endif
enddo
c write(*,*)'npp_cal:after loop iday fapar'
@@ -374,13 +379,13 @@ c calculation of g*Ca and g*Ci
c=======================================================================
do iday = 1, nd
- do ip = 1, npft
+ do ip = 1, npft
xdgca(iday) = xdgca(iday)
& + gca(ip,iday) * xlai(ip,iday) * frac(ip)
xdgci(iday) = xdgci(iday)
& + gci(ip,iday) * xlai(ip,iday) * frac(ip)
- enddo
+ enddo
enddo
c write(*,*)'npp_cal: end of subroutine'