! sim.f
      double precision S,Q,Rn,U,Us,Ta,ea,Ca,du
      double precision T,D,Cc,H,lE,A,Ra,Rc
      imax=1
      do 100 i=1,imax,1
!        read(*,*,end=101) S,Rn,Ta,Ea
        S=1000.
        Rn=800.
        Ta=20.
        Ea=10.
        Q=S*2.5
        U=1.0
        Us=0.4
        Ca=380.
        call main(Q,Rn,U,Us,Ta,Ea,Ca,T,D,Cc,H,lE,A,Ra,Rc)
	write(*,1000) Q,Rn,U,Us,Ta,Ea,Ca,T,D,Cc,H,lE,A,Ra,Rc
1000    format(15f20.4)
100   continue
101   continue
      end
!******************************************************************************
! 入力項目: Q, Rn, U, Us, Ta, ea, Ca
!           光合成有効放射、純放射、風速、摩擦速度
!           気温、水蒸気圧、大気CO2濃度
! 出力項目： T, D, Cc, H, lE, A, Ra, Rc
!           葉温、飽差、葉内CO2濃度、顕熱フラックス、潜熱フラックス
!           光合成速度、空気力学的抵抗、群落抵抗
      subroutine main(Q,Rn,U,Us,Ta,ea,Ca,T,D,Cc,H,lE,A,Ra,Rc)
      double precision Q,Rn,U,Us,Ta,ea,Ca
      double precision Ra,Rc,H,lE,A,T,D,Cc
      double precision ga,gc,gt,E,A0,T0
      call cal_Ra(U,Us,Ra)
      T=31.
      do 100 i=1,100,1
        T0=T
        call cal_D(T,ea,D)
        call cal_Rc(Q,T,D,Rc)
        call cal_T(Rn,Ta,Ra,Rc,D,T,H,lE)
        if(T.eq.T0) goto 200
100   continue
200   continue
      call cal_E(lE,T,E)
      call cal_g(Ra,Rc,T,ga,gc,gt)
      A=1.
      do 300 j=1,100,1
        A0=A
        call cal_Cc(gt,E,Ca,A,Cc)
        call cal_A(Q,T,Cc,A)
        if(A.eq.A0) goto 400
300   continue
400   continue
!      A=-44.*A/1000.
      if(abs(T).gt.9999999.) T=-9999999.
      if(abs(D).gt.9999999.) D=-9999999.
      if(abs(Cc).gt.9999999.) Cc=-9999999.
      if(abs(H).gt.9999999.) H=-9999999.
      if(abs(lE).gt.9999999.) lE=-9999999.
      if(abs(A).gt.9999999.) A=-9999999.
      if(abs(Ra).gt.9999999.) Ra=-9999999.
      if(abs(Rc).gt.9999999.) Rc=-9999999. 
      end
!******************************************************************************
! cal_Cc: gt, E, Ca, A -> Cc
! gt: 全コンダクタンス、E: 水蒸気フラックス、Ca: CO2濃度、A: 光合成
! Cc: 細胞間隙CO2濃度
      subroutine cal_Cc(gt,E,Ca,A,Cc)
      double precision gt,E,Ca,A,Cc
      if((gt.eq.0.).and.(E.eq.0.)) then
        Cc=Ca
      else
        Cc=((gt-E*1000./18./2.)*Ca-A)/(gt+E*1000./18./2.)
      end if
      end
!******************************************************************************
! cal_g: Ra, Rc, T -> ga, gc, gt 
! Ra: 空気力学的抵抗、Rc: 群落抵抗、T: 葉温、ga: 空気力学的コンダクタンス、
! gc: 群落コンダクタンス、gt: 全コンダクタンス
      subroutine cal_g(Ra,Rc,T,ga,gc,gt)
      double precision Ra,Rc,T,ga,gc,gt,P,R,T0,Rt
! 定数
      P=101300.
      R=8.31
      T0=273.15
!
      if(Ra.gt.9999.) then
        ga=0.
      else
        ga=P/Ra/R/(T+T0)
      end if
      if(Rc.gt.9999.) then
        gc=0.
      else 
        gc=P/Rc/R/(T+T0)
      end if
      Rt=1.6**(2/3)*Ra+1.6*Rc
      if(Rt.gt.9999.) then
	gt=0.
      else
	gt=P/Rt/R/(T+T0)
      end if
      end
!******************************************************************************
! cal_T: Rn, Ta, Ra, Rc, D -> T, H, lE
! Rn: 純放射、Ta: 気温、Ra: 空気力学的抵抗、Rc: 群落抵抗、D: 飽差、T: 葉温
      subroutine cal_T(Rn,Ta,Ra,Rc,D,T,H,lE)
      double precision Rn,Ta,Ra,Rc,D,T
      double precision Rou,Gamma,H,lE
! 定数
      Cp=1.01
      Rou=1289.76-4.11*T
      Gamma=0.000636*T+0.646
!
      if(Rc+Ra.lt.9999.) then
        lE=Rou*Cp*D/Gamma/(Ra+Rc)
      else
        lE=0.
      end if
      H=Rn-lE
      T=Ta+H*Ra/Rou/Cp
      if(T.lt.-20.) T=-20.
      if(T.gt.50.)  T=50.
      end
!******************************************************************************
! cal_E: lE, T -> E
! lE: 潜熱、T: 葉温、E: 水蒸気フラックス
      subroutine cal_E(lE,T,E)
      double precision lE,T,E,l
      l=(2501.-2.37*T)*1000.
      E=lE/l
      end
!******************************************************************************
! cal_Q: S -> Q
! S: 日射、Q: 光合成有効放射
      subroutine cal_Q(S,Q)
      double precision S,Q
      Q=S*2.5
      end 
!******************************************************************************
! cal_D: T, e -> D
! T: 葉温、e: 水蒸気圧、D: 飽差
      subroutine cal_D(T,ea,D)
      double precision T,ea,D
      D=6.1078*exp(17.268*T/(T+237.3))-ea
      end       
!*****************************************************************************
! cal_Ra: U, Us -> Ra
! U: 風速、Us: 摩擦速度、Ra: 空気力学的抵抗
      subroutine cal_Ra(U,Us,Ra)
      double precision U,Us,Ra
      if(Us.le.0.) then
        Ra=9999.9
      else 
        Ra=U/Us/Us
      end if
      end
!******************************************************************************
! cal_Rc: Q, T, D -> Rc
! Q: 光合成有効放射、T: 葉温、D: 飽差、Rc: 群落抵抗
      subroutine cal_Rc(Q,T,D,Rc)
      double precision T,Q,Rc,D
      double precision F_Q,F_T,F_D
      double precision Rcmin,a,To,Th,Tl,b,gc
!      double precision Pa,T0,R
! パラメータ ------------------------------------------------------------------
      Rcmin=146.9
      a=0.00000676
      To=12.0
      Th=50.0
      Tl=-5.0
      b=0.00831
! -----------------------------------------------------------------------------
! 定数
!      Pa=101300.
!      T0=273.15
!      R=8.31 
! 関数 F(Q), F(T), F(D)
      if(Q.ge.0.) then
        F_Q=Q/(Q+1/a/Rcmin)
      else
        F_Q=0.
      end if
      if((Tl.lt.-999.).and.(Th.gt.999.)) then
	F_T=1.
      else
        if((T.gt.Tl).and.(T.lt.Th)) then
          F_T=(T-Tl)/(To-Tl)*(((Th-T)/(Th-To))
     &                  **((Th-To)/(To-Tl)))
        else
          F_T=0.
	end if
      end if
      F_D=1-b*D
      if(F_D.lt.0.) F_D=0.
! 群落コンダクタンスgc と 群落抵抗Rc
      gc=F_Q*F_T*F_D/Rcmin
      if(gc.lt.0.00001) then
        Rc=9999.9
      else 
        Rc=1/gc
      end if
      end
!******************************************************************************
! cal_A: Q, T, Cc -> A
! Q: 光合成有効放射、T: 葉温, A: 光合成
      subroutine cal_A(Q,T,Cc,A)
      double precision Q,T,Cc,A
      double precision Vc,Wj,Wc,Wp,Rd
      double precision ef,Rou,Rd25
      double precision P,R,T0,Kc,Ko
      double precision pC,pG,J,Vcmax
! パラメータ ------------------------------------------------------------------
      ef=0.073
      Rou=0.22
      Rd25=0.12
      c=1.
! -----------------------------------------------------------------------------
! 定数
      P=101300.
      R=8.31
      T0=273.15
      Kc=46.*exp((T+T0-298.)*59356./(298.*R*(T+T0)))
      Ko=33000.*exp((T+T0-298.)*35948./(298.*R*(T+T0)))
      Vcmax=218.2*Rou*
     &  exp((T+T0-298)*58520/(298*R*(T+T0)))/
     &  (1+exp((703*(T+T0)-220000)/(R*(T+T0))))
      pC=P*Cc/1000000.
      pG=2205*Kc/Ko
      if(pC.lt.pG) pC=pG
      J=ef*Q/2.
      Wj=J/(4.5+10.5*pG/pC)
      Wc=Vcmax*pC/(pC+Kc*(1.+21000./Ko))
      Wp=Vcmax/2.
      if((T+T0-298.).eq.0.) then
        Rd=Rd25
      else
        Rd=Rd25*exp((T+T0-298.)*66405.*c/(298.*R*(T+T0)))/
     &  (1+exp(1.3*(T-55.)))
      end if
      call min_(Wj,Wc,Wp,Vc)
      A=Vc*(1.-pG/pC)*44./1000.-Rd
      end
!******************************************************************************
! min: x1, x2, x3 -> x
      subroutine min_(x1,x2,x3,x)
      double precision x1,x2,x3,x
      if(x1.lt.x2) then 
        x=x1
      else
        x=x2
      end if
      if(x3.lt.x) x=x3
      end
!******************************************************************************     


