Skip to content
Snippets Groups Projects
Commit b3fbe173 authored by Dethinne Thomas's avatar Dethinne Thomas
Browse files

I forgot to update one of the scripts

parent 9d028a4d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
couplage/CARAIB/ver01_Iv_couplage/carbon_npp_cal.f
+ 98
− 93
View file @ b3fbe173
......@@ -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'
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment