********************************************************************
* Equation sequence to calculate 
* species-specific deposition velocities
* according to the method of Wesely (1989)
*                                                                 *
* SPECIES-SPECIFIC  DEPOSITION PARAMETERS 
*  dH2Odx      => D_H2O/D_x          
*  kH (kHenry) => constante de Henry  
*  reacf       => facteur de reactivite 
*                                      
* OTHER INPUTS :                                                         *
*   temp     : temperature (K)                                     *
*   wind     : wind speed (m/s)            
*
* OUTPUT :                                                        *
*   Vd(i)        : deposition veolcity
*
* EXPLANATION OF PARAMETERS FOR SURFACE(i) (tabulated below)
*   rik(i)   : minimum stomatal resistance to H2Ovap 
*   rlu0k(i) : leaf cuticular resistance
*   rack(i)  : transfer resistance (depends on canopy height/density
*   rgsSk(i) : resistance de piegage au sol pour SO2              *
*   rgsOk(i) : resistance de piegage au sol pour O3               *
*   rclSk(i) : resistance des composants de la partie inferieure  *
*              de la canopee (ecorce, brindilles...) pour SO2     *
*   rclOk(i) : resistance des composants de la partie inferieure  *
*              de la canopee (ecorce, brindilles...) pour O3      *
*   z0k(i)   : roughness depth
*
*******************************************************************
* SEASONAL TABLES of RESISTANCES by SURFACE TYPE
*******************************************************************
season1 cult   decid   conif   urban
rac       10    1000    2000     100
rclS    9999    9000    3000    9999
rclO    1000     400    1000    9999
rgsO     150     200     200     300
rgsS     150     500     500     400
ri      9999    9999     250    9999
rlu0    9999    9000    4000    9999
z0      0.05       1       1       1
*******************************************************************
season2 cult   decid   conif   urban
rac      200    2000    2000     100
rclS    2000    2000    2000    9999
rclO    1000    1000    1000    9999
rgsO     150     200     200     300
rgsS     150     500     500     400
ri        60      70     130    9999
rlu0    2000    2000    2000    9999
z0      0.05       1       1       1
*******************************************************************

* INTERNAL PARAMETERS *
      pi=3.141592
* Von Karman constant
      cvk=0.4
* reference temperature (273 K)
      T0=273.
* diffusion coefficient for H2O at 273 K (m2/s)
      Dwat=0.23E-4
* air viscosity (m2/s)
      xmu=1.461E-5
* slope at ground (angle)
      gslope=0.
* reference height (height at which the wind is given) (m)
      zs=10.

* CALCULATIONS *
* solar irradiance (W/m2) in terms of solar zenith angle (degrees)
      W = MAX ( 1370.*(0.25+0.50*0.6)*cos(pi*zen/180.) , 0. )

* (LOOP OVER THE DIFFERENT SURFACES) *

* CORRECT FOR TEMPERATURE DEPENDENCE
        rt = 1000.*EXP(-(temp-273.)-4.)
        Rlu0 = Rlu0 + rt
        RclS = RclS + rt
        RclO = RclO + rt
        RgsS = RgsS + rt
        RgsO = RgsO + rt

* set maximum and minimum values
        IF (Ri.GE.9998.)   Ri=1E5        
        IF (Rlu0.GE.9998.) Rlu0=1E5        
        IF (Rac.GE.9998.)  Rac=1E5        
        IF (RgsS.GE.9998.) RgsS=1E5        
        IF (RgsO.GE.9998.) RgsO=1E5        
        IF (RclS.GE.9998.) RclS=1E5        
        IF (RclO.GE.9998.) RclO=1E5 
        IF (rac.LT.1.)     Rac=1.
        IF (rgsS.LT.1.)    RgsS=1.

* default value        
      tempfac=100.
* if 0'C < temperatures < 40 'C and Ri < 9999.:
      tempfac = 400./( (temp-273.)*(40.-(temp-273.)) )

* default value        
      Rs = Ri
* if Ri < 9999.:
      Rs = Ri*(1.+(200./(W+0.1))**2.)*tempfac

      Rdc = 100.*( 1+1000./(W+10.) )/(1.+1000.*gslope)

* aerodynamic resistance
* -------------------------------
* winter (season 1)
          xmin = 200.
          xmax = 50.
* summer (season 2)
          xmin = 50.
          xmax = 250.

        IF (zen.GE.90.) THEN
          HV = - xmin
        else
          HV = xmax*COS(pi*zen/180.)-xmin
        ENDIF

* friction velocity (without correction at first)
      ustar = wind * cvk / log( zs / z0 )

* iteration to evaluate correction function
      iter = 0
      DO WHILE(iter < 10 .AND. err > 0.01)
        iter = iter+1

* existing value of ustar, for comparison
        xx = ustar

* correction for wind, light, temperature 
        Z_L = -0.032 * HV / ((ustar**3.) * temp)

        IF (Z_L < 0.) THEN
          Z_L = MAX (Z_L, -0.99)
          ZZ = log(-Z_L)
          phiH = exp(0.032+0.448*ZZ-0.132*ZZ*ZZ)
        ELSE
          Z_L = MIN (Z_L, 0.99)
          phiH = -5.*Z_L
        ENDIF

* new (corrected) value of ustar
        ustar = wind*cvk/(log(zs/z0)-phiH)

* difference from the previous estimate of ustar
        err = abs((ustar-xx)/xx)

      END DO  

* reevaluate correction function if necessary
   IF (Z_L.LT.0) phiH = exp(0.598 + 0.39 * ZZ - 0.09 * ZZ**2.)

* evaluate Ra

   ustark = (ustar * cvk)
   Ra = (log(zs/z0)-phiH)/ ustark

* laminar resistance 

   DH2O = Dwat * (temp/T0)**1.75
   Rb = (2./ustark)*((xmu/DH2O)**(2./3.))*dH2Odx**(2./3.)) 

* --------------------------
* Rc, surface resistance (CHOOSE ONE OF FOLLOWING 3 OPTIONS)
* surface resistance for O3
   Rc(O3) = 1./(1./(Rs*depdatspe(j,1)) + 1./rlu0 + 1./(Rdc+RclO) + 1./(Rac+RgsO))

* surface resistance for SO2
   Rc(SO2) = 1./ (1./(Rs*depdatspe(j,1)) + 1./rlu0 + 1./(Rdc+RclS) + 1./(Rac+RgsS))

* surface resistance for all other species
   Rsm = Rs * dH2Odx + 1./( kH/3000. + 100.*reacf )
   Rlu = rlu0/( kH*1E-5 + reacf )
   Rcl = 1./( (kH/1E5/RclS) + (reacf/RclO) )
   Rgs = 1./( (kH/1E5/RgsS) + (reacf/RgsO) )
   Rc  = 1./( (1./Rsm)+(1./Rlu)+(1./(Rdc+Rcl))+(1./(Rac+Rgs)) )
* --------------------------
                
* surface resistance minimum and maximum values
   Rc = MAX (MIN (Rc,9999.),10.)

* => deposition velocity (cm/s)
  Vdtyp = 100./ (Ra + Rb + Rc) 

* sum species-specific deposition velocities over all surfaces present
  Vd = Vd + (xsurf(i)*Vdtyp(i),i=1:4)

* END *