stomata.f90

CLM集合卡曼滤波数据同化算法

!     assimilation rate. ! pay attention on this iteration 

      co2m = pco2m/psrf                               ! mol mol-1
      co2a = pco2a/psrf

      range    = pco2m * ( 1. - 1.6/gradm ) - gammas

      do ic = 1, 6
      pco2y(ic) = 0.
      eyy(ic) = 0.
      enddo

      ITERATION_LOOP: do ic = 1, 6                                                         

      call sortin(eyy, pco2y, range, gammas, ic)      ! legacy from SiB2
      pco2i =  pco2y(ic)

!-----------------------------------------------------------------------
!                      NET ASSIMILATION
!     the leaf assimilation (or gross photosynthesis) rate is described 
!     as the minimum of three limiting rates:
!     omc: the efficiency of the photosynthetic enzyme system (Rubisco-limited);
!     ome: the amount of PAR captured by leaf chlorophyll;
!     oms: the capacity of the leaf to export or utilize the products of photosynthesis.
!     to aviod the abrupt transitions, two quadratic equations are used:
!             atheta*omp^2 - omp*(omc+ome) + omc*ome = 0
!         btheta*assim^2 - assim*(omp+oms) + omp*oms = 0
!-----------------------------------------------------------------------
 
      atheta = 0.877
      btheta = 0.95

      omc = vm   * ( pco2i-gammas ) / ( pco2i + rrkk ) * c3 + vm * c4 
      ome = epar * ( pco2i-gammas ) / ( pco2i+2.*gammas ) * c3 + epar * c4 
      oms   = omss * c3 + omss*pco2i * c4 

      sqrtin= max( 0., ( (ome+omc)**2 - 4.*atheta*ome*omc ) ) 
      omp   = ( ( ome+omc ) - sqrt( sqrtin ) ) / ( 2.*atheta )
      sqrtin= max( 0., ( (omp+oms)**2 - 4.*btheta*omp*oms ) )
      assim = ( ( oms+omp ) - sqrt( sqrtin ) ) / ( 2.*btheta )

      assimn= ( assim - respc)                        ! mol m-2 s-1
                                                                               
!-----------------------------------------------------------------------
!                      STOMATAL CONDUCTANCE
!
!  (1)   pathway for co2 flux
!                                                  co2m
!                                                   o
!                                                   |   
!                                                   |    
!                                                   <  |
!                                        1.37/gsh2o >  |  Ac-Rd-Rsoil
!                                                   <  v
!                                                   |
!                                     <--- Ac-Rd    |
!     o------/\/\/\/\/\------o------/\/\/\/\/\------o
!    co2i     1.6/gsh2o     co2s    1.37/gbh2o     co2a
!                                                   | ^                        
!                                                   | | Rsoil                        
!                                                   | |                        
!
!  (2)   pathway for water vapor flux 
!
!                                                  em
!                                                   o
!                                                   |   
!                                                   |    
!                                                   <  ^
!                                           1/gsh2o >  | Ea
!                                                   <  |
!                                                   |
!                                     ---> Ec       !
!     o------/\/\/\/\/\------o------/\/\/\/\/\------o
!     ei       1/gsh2o      es       1/gbh2o       ea
!                                                   | ^                        
!                                                   | | Eg
!                                                   | |                        
!
!  (3)   the relationship between net assimilation and tomatal conductance :
!        gsh2o = m * An * [es/ei] / [pco2s/p] + b
!        es = [gsh2o *ei + gbh2o * ea] / [gsh2o + gbh2o]
!        ===>
!        a*gsh2o^2 + b*gsh2o + c = 0
!
!-----------------------------------------------------------------------

!-->  co2a = co2m - 1.37/max(0.446,gah2o) * (assimn - 0.)     ! mol mol-1
      co2s = co2a - 1.37*assimn/gbh2o                  ! mol mol-1

      co2st = min( co2s, co2a )
      co2st = max( co2st,1.e-5 )
                               
      assmt = max( 1.e-12, assimn )
      hcdma = ei*co2st / ( gradm*assmt ) 

      aquad = hcdma                     
      bquad = gbh2o*hcdma - ei - bintc*hcdma
      cquad = -gbh2o*( ea + hcdma*bintc )  
                                             
      sqrtin= max( 0., ( bquad**2 - 4.*aquad*cquad ) )
      gsh2o = ( -bquad + sqrt ( sqrtin ) ) / (2.*aquad) 

      es  = ( gsh2o-bintc ) * hcdma                   ! pa 
      es  = min( es, ei )
      es  = max( es, 1.e-2)

      gsh2o = es/hcdma + bintc                        ! mol m-2 s-1                                      

      pco2in = ( co2s - 1.6 * assimn / gsh2o )*psrf   ! pa
      eyy(ic) = pco2i - pco2in                        ! pa

!-----------------------------------------------------------------------        
                                                                               
      if( abs(eyy(ic)) .lt. 0.1 ) exit
                                                                               
      enddo ITERATION_LOOP

! convert gsh2o (mol m-2 s-1) to resistance rst ( s m-1)
      rst   = min( 1.e6, 1./(gsh2o*tlef/tprcor) )     ! s m-1                  


      END SUBROUTINE stomata


      SUBROUTINE sortin( eyy, pco2y, range, gammas, ic )                        
                                                                               
!=======================================================================        
!     arranges successive pco2/error pairs in order of increasing pco2.         
!     estimates next guess for pco2 using combination of linear and             
!     quadratic fits.                                                           
!
!     original author: P. J. Sellers (SiB2)
!-----------------------------------------------------------------------        
      implicit none
 
      integer, intent(in) :: &
      ic           !

      real, INTENT(in) :: &
      range,      &!
      gammas       !

      real, INTENT(inout), dimension(6) :: &
      eyy,        &!
      pco2y        !                                            

!----- Local -----------------------------------------------------------        
      integer i, j, n, i1, i2, i3, is, isp, ix
      real a, b, pmin, emin
      real pco2b, pco2yl, pco2yq
      real ac1, ac2, bc1, bc2, cc1, cc2
      real bterm, aterm, cterm
!-----------------------------------------------------------------------        
                                                                               
      if( ic .ge. 4 ) go to 500                                                 
      pco2y(1) = gammas + 0.5*range                                             
      pco2y(2) = gammas + range*( 0.5 - 0.3*sign(1.0,eyy(1)) )                  
      pco2y(3) = pco2y(1) - (pco2y(1)-pco2y(2))/(eyy(1)-eyy(2)+1.e-10)*eyy(1)   
                                                                               
      pmin = min( pco2y(1), pco2y(2) )                                        
      emin = min(   eyy(1),   eyy(2) )                                        
      if ( emin .gt. 0. .and. pco2y(3) .gt. pmin ) pco2y(3) = gammas            
      go to 200                                                                 
500   continue                                                                  
                                                                               
      n = ic - 1                                                                
      do 1000 j = 2, n                                                          
      a = eyy(j)                                                                
      b = pco2y(j)                                                              
      do 2000 i = j-1,1,-1                                                      
      if(eyy(i) .le. a ) go to 100                                              
      eyy(i+1) = eyy(i)                                                         
      pco2y(i+1) = pco2y(i)                                                     
2000  continue                                                                  
      i = 0                                                                     
100   eyy(i+1) = a                                                              
      pco2y(i+1) = b                                                            
1000  continue                                                                  
                                                                               
      pco2b = 0.                                                                
      is    = 1                                                                 
      do 3000 ix = 1, n                                                         
      if( eyy(ix) .lt. 0. ) pco2b = pco2y(ix)                                   
      if( eyy(ix) .lt. 0. ) is = ix                                             
3000  continue                                                                  
      i1 = is-1                                                                 
      i1 = max(1, i1)                                                          
      i1 = min(n-2, i1)                                                        
      i2 = i1 + 1                                                               
      i3 = i1 + 2                                                               
      isp   = is + 1                                                            
      isp = min( isp, n )                                                      
      is = isp - 1                                                              
                                                                               
      pco2yl=pco2y(is) - (pco2y(is)-pco2y(isp))/(eyy(is)-eyy(isp))*eyy(is)
                                                                               
!----------------------------------------------------------------------         
!   method using a quadratic fit                                                
!----------------------------------------------------------------------         
                                                                               
      ac1 = eyy(i1)*eyy(i1) - eyy(i2)*eyy(i2)                                   
      ac2 = eyy(i2)*eyy(i2) - eyy(i3)*eyy(i3)                                   
      bc1 = eyy(i1) - eyy(i2)                                                   
      bc2 = eyy(i2) - eyy(i3)                                                   
      cc1 = pco2y(i1) - pco2y(i2)                                               
      cc2 = pco2y(i2) - pco2y(i3)                                               
      bterm = (cc1*ac2-cc2*ac1)/(bc1*ac2-ac1*bc2)                               
      aterm = (cc1-bc1*bterm)/ac1                                               
      cterm = pco2y(i2) - aterm*eyy(i2)*eyy(i2) - bterm*eyy(i2)                 
      pco2yq= cterm                                                             
      pco2yq= max( pco2yq, pco2b )                                            
      pco2y(ic) = ( pco2yl+pco2yq)/2.                                           
                                                                               
200   continue                                                                  
                                                                               
      pco2y(ic) = max ( pco2y(ic), 0.01 )                                     
                                                                               
      END SUBROUTINE sortin