ecearth3/sources/tm5mp/proj/cb05/dry_deposition.F90

!#################################################################
!
! Purpose:
! --------
!
! This module  contains all subroutines needed for dry deposition calculations.
! The mean purpose is to perform on a high resolution of  1x1 degree:
!  
! 0. allocate the vd on the model resolution
!
! 1. read fields needed for further calculations in subroutines:
!    
! 2. calculate the tracer independent friction velocity, aerodynamic 
!    and sub-laminar resistance 
!
! 3. calculate fields needed for resistance calculations 
!
! 4. the tracer dependent vd
!
! 5. coarsen the vd
!
! 6. apply the coarsened deposition velocities to concentration field rm
!
! 7. deallocate vd
!
! Reference:
! ---------
! general documentationon ECMWF fields can be found on:
! http://www.ecmwf.int/research/ifsdocs/PHYSICS/
! ---------
! Authors:
! -------
! original code by Laurens Ganzeveld and Ad Jeuken (1997)
! revised and adapted for use in TM5 and use of ECMWF data 
! by F. Dentener and M. Krol (2003)
! 
!   24 Jan 2012 - P. Le Sager - modified for mpi lat-lon domain decomposition
!
!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

module dry_deposition

  use GO          , only : gol, goPr, goErr
  use dims        , only : nregions
  use chem_param  , only : ntrace
  use global_types, only : emis_data

  implicit none


  ! --- in/out ----------------------

  private

  public :: Dry_Deposition_Init, Dry_Deposition_Done
  public :: dd_surface_fields
  
  public :: vd
  public :: vda


  ! --- const -----------------------
  
  character(len=*), parameter  ::  mname = 'dry_deposition'


  ! --- var -------------------------

  ! vd      : final deposition velocities for all species on model resolution.
  ! vd_temp : to store the vd temporarily for one specie
  type(emis_data),dimension(nregions,ntrace) :: vd
  type(emis_data),dimension(nregions)        :: vd_temp, vd_temp_global

  ! vda     : final deposition velocities for all aerosols on model resolution.
  type(emis_data), pointer    ::  vda(:,:)


contains 


  !
  ! Allocate emission data structure "vd",
  ! final deposition velocities for all species on model resolution
  !

  subroutine Dry_Deposition_Init( status )

    use dims      , only : nregions, iglbsfc
    use chem_param, only : ndep
    use chem_param, only : nrdep
    use meteodata , only : lsmask_dat
    use meteodata , only : sr_ecm_dat, sr_ols_dat, ci_dat, sd_dat, swvl1_dat
    use meteodata , only : nsss_dat, ewss_dat, wspd_dat, slhf_dat, sshf_dat
    use meteodata , only : src_dat, d2m_dat, t2m_dat, ssr_dat
    use meteodata , only : nveg
    use meteodata , only : tv_dat, cvl_dat, cvh_dat
    use meteodata , only : Set
    use tm5_distgrid, only : dgrid, Get_DistGrid
    use ParTools,     only : isRoot

    ! --- in/out ------------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const -------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/Dry_Deposition_Init'
    
    ! --- local -------------------------------------

    integer  :: imr,jmr,region,n
    integer  :: iveg, i1, i2, j1, j2
    
    ! --- begin -------------------------------------
    
    ! deposition velocities ; 
    ! always allocated, filled with zeros if no tracers are to be deposited:
    do region = 1, nregions

       CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       do n=1,ntrace          
          allocate(vd(region,n)%surf(i1:i2, j1:j2))
          vd(region,n)%surf = 0.0          
       end do
       
    end do

    ! deposited tracers ?
    if ( ndep > 0 ) then

      ! select meteo to be used
      call Set( lsmask_dat(iglbsfc), status, used=.true. )
      call Set( sr_ecm_dat(iglbsfc), status, used=.true. )
      call Set( sr_ols_dat(iglbsfc), status, used=.true. )
      call Set(     ci_dat(iglbsfc), status, used=.true. )
      call Set(     sd_dat(iglbsfc), status, used=.true. )
      call Set(  swvl1_dat(iglbsfc), status, used=.true. )
      call Set(   ewss_dat(iglbsfc), status, used=.true. )
      call Set(   nsss_dat(iglbsfc), status, used=.true. )
      call Set(   wspd_dat(iglbsfc), status, used=.true. )
      call Set(   slhf_dat(iglbsfc), status, used=.true. )
      call Set(   sshf_dat(iglbsfc), status, used=.true. )
      call Set(    src_dat(iglbsfc), status, used=.true. )
      call Set(    d2m_dat(iglbsfc), status, used=.true. )
      call Set(    t2m_dat(iglbsfc), status, used=.true. )
      call Set(    ssr_dat(iglbsfc), status, used=.true. )
      call Set(    cvl_dat(iglbsfc), status, used=.true. )
      call Set(    cvh_dat(iglbsfc), status, used=.true. )
      do iveg = 1, nveg
        call Set( tv_dat(iglbsfc,iveg), status, used=.true. )
      end do

      ! temporary storage:
      do region = 1, nregions

         CALL GET_DISTGRID( DGRID(region), I_STRT=i1,  I_STOP=i2, &
                                           J_STRT=j1,  J_STOP=j2, &
                                           I_WRLD=imr, J_WRLD=jmr )

         allocate(vd_temp(region)%surf(i1:i2,j1:j2))
         
         if(isRoot) then
            allocate(vd_temp_global(region)%surf(imr,jmr))
         else
            allocate(vd_temp_global(region)%surf(1,1))
         end if
         
      end do
      
    end if   ! deposited tracers

    ! aerosols ?
    if ( nrdep > 0 ) then

       allocate( vda(nregions,nrdep) )
       
       do region = 1, nregions

          CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

          do n=1,nrdep

             allocate(vda(region,n)%surf(i1:i2,j1:j2))
             
             vda(region,n)%surf = 0.0
             
          end do
       end do
    end if
    
    ! ok
    status = 0

  end subroutine Dry_Deposition_Init
  
  
  ! ***


  !
  ! De-allocate emission data structure "vd",
  ! final deposition velocities for all species on model resolution
  !

  subroutine Dry_Deposition_Done( status )

    use dims      , only : nregions
    use chem_param, only : ndep
    use chem_param, only : nrdep

    ! --- in/out ------------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const -------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/Dry_Deposition_Done'
    
    ! --- local -------------------------------------

    integer  :: region, n
    
    ! --- begin -------------------------------------

    ! clear deposition velocities:
    do region = 1, nregions
      do n=1,ntrace
        deallocate(vd(region,n)%surf)
      end do
    end do

    ! clear temporary storage:
    if ( ndep > 0 ) then
      do region = 1, nregions
         deallocate(vd_temp(region)%surf)
         deallocate(vd_temp_global(region)%surf)
      end do
    end if   ! deposited tracers

    ! aerosols ?
    if ( nrdep > 0 ) then
      do region = 1, nregions
        do n = 1, nrdep
          deallocate( vda(region,n)%surf )
        end do
      end do
      deallocate( vda )
    end if
    
    ! ok
    status = 0

  end subroutine Dry_Deposition_Done


  !------------------------------------
  !
  ! Purpose:
  ! -------
  ! this subroutine reads and prepares all surface fields
  !
  ! External
  ! --------
  ! dd_get_3_hourly_surface_fields
  ! dd_calc_ustar_raero_rb
  ! dd_calc_rstom_rahcan
  ! dd_calc_inisurf
  ! dd_calc_rs
  ! dd_coarsen_vd
  !
  ! Reference
  ! ---------
  ! Ganzeveld and Lelieveld (1996) and references therein.
  !
  !------------------------------------

  subroutine dd_surface_fields( status )

    use binas          , only : vkarman
    use dims           , only : idate
    use dims           , only : okdebug
    use chem_param     , only : names
    use chem_param     , only : ndep
    use chem_param     , only : nrdep
#ifndef without_photolysis
    use photolysis_data, only : phot_dat
#endif
    use ParTools,        only : isRoot
    use dims           , only : iglbsfc, nlat180, nlon360
    use meteodata      , only : lsmask_dat
    use meteodata      , only : sr_ecm_dat, sr_ols_dat, ci_dat, sd_dat, swvl1_dat
    use meteodata      , only : ewss_dat, nsss_dat, wspd_dat, slhf_dat, sshf_dat
    use meteodata      , only : src_dat, d2m_dat, t2m_dat, ssr_dat
    use tm5_distgrid,    only : dgrid, Get_DistGrid, scatter, gather, dist_arr_stat
    
    ! --- in/out -----------------------------
    
    integer, intent(out)    ::  status

    ! --- const -------------------------------
    
    character(len=*), parameter         ::  rname = 'dd_surface_fields'

    ! --- local ------------------------------
    
    character(len=10)                   :: c_time
    real, dimension(:,:,:), allocatable :: soilph
    real, dimension(:,:),   allocatable :: lai, lai1
    real, dimension(:,:),   allocatable :: vgrat_low, vgrat_high
    real, dimension(:,:),   allocatable :: sstr, wind10m
    real, dimension(:,:),   allocatable :: ustar_loc,rb,raero
    real, dimension(:,:,:), allocatable :: vd11  
    real, dimension(:,:),   allocatable :: vd11a
    real, dimension(:,:),   allocatable :: rahcan
    real, dimension(:,:),   allocatable :: rstom
    real, dimension(:,:),   allocatable :: snow_cover
    real, dimension(:,:),   allocatable :: wet_skin
    real, dimension(:,:),   allocatable :: fws
    real, dimension(:,:),   allocatable :: rh2m
    real, dimension(:,:),   allocatable :: ddepvel_h2
#ifndef without_photolysis
    real, dimension(:,:),   allocatable :: ags
#endif
    real, dimension(:,:),   allocatable :: alpha, bubble, alphae
    real, dimension(:,:),   allocatable :: ustar_land, ustar_sea, sr_sea
    real, dimension(:,:),   allocatable :: global_field
    
    integer :: i, nsend, region,itrace
    integer :: i1, i2, j1, j2, imr, jmr

    ! no deposited tracers or aerosols ? then leave immediately:
    if ( ndep+nrdep == 0 ) then
      status=0; return
    end if

    ! From here until the remapping (coarsen), we work on the extra iglbsfc
    ! grid, on which the data to read are expected to be.
    CALL GET_DISTGRID( DGRID(iglbsfc), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, I_WRLD=imr, J_WRLD=jmr )

    if (isRoot) then
       allocate(global_field(imr,jmr))
    else
       allocate(global_field(1,1))
    end if
    
    allocate(soilph    (i1:i2, j1:j2,5))
    allocate(lai       (i1:i2, j1:j2)  )
    allocate(lai1      (i1:i2, j1:j2)  )
    allocate(vgrat_low (i1:i2, j1:j2)  )
    allocate(vgrat_high(i1:i2, j1:j2)  )
    allocate(sstr      (i1:i2, j1:j2)  )
    allocate(wind10m   (i1:i2, j1:j2)  )
    allocate(ustar_loc (i1:i2, j1:j2)  )
    allocate(raero     (i1:i2, j1:j2)  )
    allocate(rb        (i1:i2, j1:j2)  )
    allocate(rahcan    (i1:i2, j1:j2)  )
    allocate(rstom     (i1:i2, j1:j2)  )
    allocate(snow_cover(i1:i2, j1:j2)  )
    allocate(wet_skin  (i1:i2, j1:j2)  )
    allocate(fws       (i1:i2, j1:j2)  )
    allocate(rh2m      (i1:i2, j1:j2)  )
    allocate(ddepvel_h2(i1:i2, j1:j2)  )
#ifndef without_photolysis
    allocate(ags       (i1:i2, j1:j2)  )
#endif
    allocate(vd11      (i1:i2, j1:j2,ndep))

    if ( nrdep > 0 ) allocate(vd11a(i1:i2, j1:j2))

    allocate(ustar_land(i1:i2, j1:j2))
    allocate(ustar_sea (i1:i2, j1:j2))
    allocate(sr_sea    (i1:i2, j1:j2))
    allocate(alpha     (i1:i2, j1:j2))
    allocate(bubble    (i1:i2, j1:j2))
    allocate(alphae    (i1:i2, j1:j2))

    ! -- for output of time header on debug fields

    write(c_time,'(i4,3(i2))') idate(1:4)   !CMK corrected
    do i=1,10
       if (c_time(i:i)==' ') c_time(i:i)='0'
    end do
    
    !
    ! surface data sets may have the following characteristics
    ! 1. instanteneous. 
    !      The time written in the file is valid from t-1.5 to t+1.5
    ! 2. Accumulated. 
    !      The time written in the file is valid from t to t+3
    ! 3. Daily averaged. 
    !      The time on the file is always 00 hours, and valid until t+24
    !
    ! ***   It is essential to understand this timing when 
    !       reading and applying these sets.
    !
    ! fd2mk it has to be decided to 'save' fields like 
    ! vgrat_low and daily average surface fields
    ! for later use, or to read it every 3 hours time step
    !

    ! Read data
    call dd_get_soilph_lai           ! non-ECMWF provided data: soilph and lai
    IF_NOTOK_RETURN(status=1)
    
    call dd_get_surface_vegetation   ! vgrat_low, vgrat_high and lsm
    IF_NOTOK_RETURN(status=1)
    
    wind10m = wspd_dat(iglbsfc)%data(:,:,1)
    sstr = sqrt( ewss_dat(iglbsfc)%data(:,:,1)**2 + nsss_dat(iglbsfc)%data(:,:,1)**2 )

    ! AS: ensure being non-zero
    where( sstr .lt. 1.E-10 ) sstr = 1.E-10

    call dd_calc_ustar_raero_rb         ! raero, ustar_loc, rb
    IF_NOTOK_RETURN(status=1)
    
    call dd_calc_rstom_rahcan           ! rstom, rahcan

    call dd_calc_inisurf                ! snow_cover, wet_skin, fws, rh2m, ags

    call dd_calc_rs                     ! vd
    IF_NOTOK_RETURN(status=1)
    
    ! coarsen vd to model resolution
#ifndef without_photolysis
    call dd_coarsen_vd( vd11, ags, i1, j1, status )
#else
    call dd_coarsen_vd( vd11, i1, j1, status )
#endif
    IF_NOTOK_RETURN(status=1)

    deallocate(soilph)
    deallocate(lai)
    deallocate(lai1)
    deallocate(vgrat_low)
    deallocate(vgrat_high)

    deallocate(sstr)
    deallocate(wind10m)

    deallocate(ustar_loc)
    deallocate(raero)
    deallocate(rb)
    deallocate(rahcan)
    deallocate(rstom)
    deallocate(snow_cover)
    deallocate(wet_skin)
    deallocate(fws)
    deallocate(rh2m)
    deallocate(ddepvel_h2)
#ifndef without_photolysis
    deallocate(ags)
#endif
    deallocate(vd11)
    if ( nrdep > 0 ) deallocate( vd11a )

    deallocate(alpha)
    deallocate(bubble)
    deallocate(alphae)
    deallocate(ustar_land)
    deallocate(ustar_sea)
    deallocate(sr_sea)
    !FD23012004

    deallocate(global_field)


! PLS #ifdef MPI
! PLS     do region = 1, nregions
! PLS        nsend = im(region)*jm(region) 
! PLS        do itrace = 1, ntrace
! PLS           call mpi_bcast(vd(region,itrace)%surf, nsend, my_real, root_k, &
! PLS                com_lev, ierr)
! PLS        end do
! PLS        ! aerosols ?
! PLS        if ( nrdep > 0 ) then
! PLS          do itrace = 1, nrdep
! PLS            call mpi_bcast( vda(region,itrace)%surf, nsend, MY_REAL, root_k, com_lev, ierr )
! PLS          end do
! PLS        end if
! PLS #ifndef without_photolysis
! PLS        ! send albedo computed in dd_calc_inisurf to all other processors:
! PLS        call mpi_bcast(phot_dat(region)%albedo ,nsend, my_real, root_k, com_lev, ierr)
! PLS #endif
! PLS     end do
! PLS     if ( myid == root_k .and. okdebug ) then
! PLS        write (gol,'("dd_surface_fields: Deposition velocities broadcasted")') ; call goPr
! PLS     end if
! PLS #endif

    ! ok
    status = 0

  contains

    !---------------------------
    !
    ! Purpose
    ! -------
    ! Reads two independent datasets on soils and vegetation not provided by ECMWF
    ! 
    ! Reference:
    ! ---------
    ! Data provided by Laurens Ganzeveld
    !
    !--------------------------

    SUBROUTINE DD_GET_SOILPH_LAI

      use global_data, only : inputdir
      USE MDF, ONLY : MDF_Open, MDF_NETCDF, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close

      character(len=5) :: vname
      character(len=8) :: vname8
      character(len=9) :: vname9
      character(len=5),dimension(5)  :: name_soil = (/'spec1','spec2','spec3','spec4','spec5'/)
      integer :: i, mo, n, varid, fid
      real, dimension(:,:), allocatable :: pland
      
      ! work array
      if(isRoot)then
        allocate( pland(nlon360,nlat180) )
      else
        allocate( pland(1,1) )
      end if
      
      !
      !  -- Reading of input file with LAI data which is derived
      ! from the Olson ecosystem database (0.5 X 0.5 degrees), which
      ! discerns 46 ecosystems and their characteristics, e.g. LAI. 
      ! Fields are monthly averaged.
      ! 
      mo = idate(2) - 1

      if (isRoot)then
        CALL MDF_Open( TRIM(inputdir)//'/land/lsmlai.nc4', MDF_NETCDF, MDF_READ, fid, status )
        IF_NOTOK_RETURN(status=1)
      endif
      
      if (isRoot) then
        if ( mo < 10 ) then
          write(vname8,'("spec1__",i1)')mo
          CALL MDF_Inq_VarID( fid, vname8, varid, status ) 
          IF_NOTOK_RETURN(status=1)
        else
          write(vname9,'("spec1__",i2)')mo
          CALL MDF_Inq_VarID( fid, vname9, varid, status ) 
          IF_NOTOK_RETURN(status=1)
        endif
        CALL MDF_Get_Var( fid, varid, pland, status )
        IF_NOTOK_RETURN(status=1)
      endif
      call scatter( dgrid(iglbsfc), lai, pland, 0, status)
      IF_NOTOK_RETURN(status=1)

      if (isRoot)then
        CALL MDF_Close( fid, status )
        IF_NOTOK_RETURN(status=1)
      endif

      if ( okdebug ) then
        call dist_arr_stat(dgrid(iglbsfc), 'lai', lai, 0, status)
        IF_NOTOK_RETURN(status=1)
      end if
      
      !
      ! -- Reading of input file with soil pH data,which are derived from
      !    a 0.5 x 0.5 degree soil pH database, Batjes, 1995  
      !    over land soilph 1 to 5 should add up to 1.
      !    SOILPH(I,J,1) - soil pH <5.5
      !    SOILPH(I,J,2) - soil 5.5 <pH <7.3
      !    SOILPH(I,J,3) - soil 7.3< pH <8.5
      !    SOILPH(I,J,4) - soil 8.5 <pH
      !    SOILPH(I,J,5) - soil 4 < pH <8.5
      !
      if (isRoot)then
        CALL MDF_Open( TRIM(inputdir)//'/land/soilph.nc4', MDF_NETCDF, MDF_READ, fid, status )
        IF_NOTOK_RETURN(status=1)
      endif

      do n=1,5 
        write(vname,'("spec",i1)')n
        if (isRoot) then
          CALL MDF_Inq_VarID( fid, vname, varid, status ) 
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Get_Var( fid, varid, pland, status )
          IF_NOTOK_RETURN(status=1)
        endif
        call scatter( dgrid(iglbsfc), soilph(:,:,n), pland, 0, status)
        IF_NOTOK_RETURN(status=1)
      end do
      
       if (isRoot)then
         CALL MDF_Close( fid, status )
         IF_NOTOK_RETURN(status=1)
       endif

       ! work array
       deallocate(pland)
       
       if ( okdebug ) write(*,*) 'dd_get_soilph_lai: end'

    END SUBROUTINE DD_GET_SOILPH_LAI
    !
    !
    subroutine dd_get_surface_vegetation
      !-----------------------------------
      ! Purpose
      ! -------
      ! read ECMWF dataset for vegetation and landmask
      ! 
      !------------------------
      !
      ! routine to read fields that need to be initialised only once
      !
      use meteodata, only : nveg
      use meteodata, only : tv_dat, cvl_dat, cvh_dat
      use datetime,  only : date2tau, tau2date
      !
      ! following parameters are directly from Table 7.1 
      ! www.ecmwf.int/research/ifsdocs/Physics
      !
      ! cveg is fractional coverage given a certain type of vegetation present
      real,dimension(nveg),parameter       :: cveg=& 
           (/0.9, 0.85,  0.9,  0.9,  0.9, &
             0.99, 0.7,  0. ,  0.5,  0.9, &
             0.1,   9.,  0.6,  9.0,  9.0, &
             0.5,  0.5,  0.9,  0.9,  0.6 /)
      !high_low: e.g. high are trees, low schrub
      character(len=1),dimension(nveg),parameter :: high_low=&
           (/'L','L','H','H','H', &
             'H','L','O','L','L', &
             'L','O','L','O','O', &
             'L','L','H','H','O' /)
      !lai_veg assigned lai per type of vegetation. No seasonal cycle
      real,dimension(nveg),parameter       :: lai_veg=&
           (/3.0, 2.0, 5.0, 5.0, 5.0, &
             6.0, 2.0, 0.5, 1.0, 3.0, &
             0.5,-9.0, 4.0,-9.0,-9.0, &
             3.0, 1.5, 5.0, 2.5, 4.0 /)
      !evergreen is the vegetation seasonal or evergren
      integer,dimension(nveg),parameter :: evergreen =& 
           (/0,0,1,0,1, &
             0,0,0,0,0, &
             0,0,0,0,0, &
             1,0,0,0,0 /)
      character(len=30),dimension(nveg),parameter :: vegtype=(/&
           'Crops and mixed farming       ','Short grass                   ',&
           'Evergreen needle leaf trees   ','Deciduous needle leaf trees   ',&
           'Evergreen broadleaf trees     ','Deciduous broad leaf trees    ',&
           'Tall grass                    ','Desert                        ',&
           'Tundra                        ','Irrigated crops               ',&
           'Semidesert                    ','Ice caps and glaciers         ',&
           'Bogs and marshes              ','Inland water                  ',&
           'Ocean                         ','Evergreen shrubs              ',&
           'Deciduous shrubs              ','Mixed forest/woodland         ',&
           'Interrupted forest            ','Water and land mixtures       '/)
      !
      ! --
      !
      integer,dimension(6)                         :: idater_acc,idater

      character(len=2),dimension(nveg),parameter:: tvname=&
           (/'01','02','03','04','05', &
             '06','07','08','09','10', &
             '11','12','13','14','15', &
             '16','17','18','19','20' /)

      integer    :: i,j,nv,itaux

      ! ----------------- START

      !
      ! -- vgrat_low :  fractional coverage of low vegetation
      ! -- vgrat_high : fractional coverage of high vegetation
      ! -- lai1      : lai derived following ECMWF, 
      !                unfortunately there is no yearly cycle on this product
      ! -- fd Note:
      !    ewsmx : maximum field capacity could possibly be approximated 
      !            with the root depth
      !  
      !   ewsmx(:,:)=ewsmx(:,:)+cvl_sfc%data(:,:,1)*tv_sfc(nv)%data(:,:,1)/100.*
      !              root_depth(nv)/depth_tessel
      !   ewsmx(:,:)=ewsmx(:,:)+cvh_sfc%data(:,:,1)*tv_sfc(nv)%data(:,:,1)/100.*
      !              root_depth(nv)/depth_tessel
      !   Note that in the ECMWF model a constant value of 0.353 m3/m3 is used
      !   depth_tessel = 2.89  ! tessel soil model has four layers, 
      !   the lowest of four ends at 2.89 m 
      !   root depth is derived from ECMWF root distribution, 
      !   and is used to calculate wilting point
      !   real,dimension(nveg),parameter :: root_depth=(/& 
      !   0.4017,0.3455,0.4223,0.4391,0.4521,0.5479,0.4401,0.0700,0.1850,0.4017,&
      !   0.6737,0.0000,0.4912,0.0000,0.0000,0.5156,0.5156,0.5350,0.5350,0.0000/)
      ! -- fd not implemented

      vgrat_high = 0.0
      vgrat_low = 0.0
      lai1 = 0.0
      do nv = 1, 20 
         if ( high_low(nv) == 'L' ) then
            vgrat_low = vgrat_low + cvl_dat(iglbsfc)%data(:,:,1) * tv_dat(iglbsfc,nv)%data(:,:,1)/100.0 * cveg(nv)
            lai1 = lai1 + cvl_dat(iglbsfc)%data(:,:,1) * tv_dat(iglbsfc,nv)%data(:,:,1)/100.0 * cveg(nv) * lai_veg(nv) 
         end if
         if ( high_low(nv) == 'H' ) then
            vgrat_high = vgrat_high + cvh_dat(iglbsfc)%data(:,:,1) * tv_dat(iglbsfc,nv)%data(:,:,1)/100.0 * cveg(nv)
            lai1 = lai1 + cvh_dat(iglbsfc)%data(:,:,1) * tv_dat(iglbsfc,nv)%data(:,:,1)/100.0 * cveg(nv) * lai_veg(nv)
         end if
      end do !nv

      if ( okdebug ) then
         call dist_arr_stat(dgrid(iglbsfc), 'vgrat_low',  vgrat_low , 0, status)
         IF_NOTOK_RETURN(status=1)
         call dist_arr_stat(dgrid(iglbsfc), 'vgrat_high', vgrat_high, 0, status)
         IF_NOTOK_RETURN(status=1)
         call dist_arr_stat(dgrid(iglbsfc), 'lai1',       lai1      , 0, status)
         IF_NOTOK_RETURN(status=1)
         write(*,*) 'dd_get_surface_vegetation: end'
      end if

    end subroutine dd_get_surface_vegetation


    !--------------------------------
    !
    ! Purpose
    ! -------
    ! Calculate friction velocity (ustar), aerodynamic resistance (raero)
    ! and quasi laminar surface resistance (rb), 
    !
    ! Method
    ! ------
    ! uses the log normal wind profile information stored by ECMWF 
    ! in 10 meter wind
    ! to estimate a local ustar over land
    ! uses Charnock equation over sea to estimate ustar
    ! aerodynamic resistance is calculated from heat fluxes and ustar 
    ! using a constant reference height
    ! sub laminar resistance from ustar
    ! 
    ! External
    ! --------
    ! dumpfield
    !
    ! Reference
    ! ----------
    ! Ge Verver (personal communication, 2003)
    !------------------------

    subroutine dd_calc_ustar_raero_rb

      use binas, only : grav, vKarman

      ! Href: constant reference height for calculations of raero
      real, parameter :: Href=30.     
      ! some constants specific to calculation of ustar and raero
      real,parameter  :: alfa_charnock1 = 0.11
      real,parameter  :: rhoCp          = 1231.0
      real,parameter  :: rhoLv          = 3013000.0
      real,parameter  :: alfa_charnock2 = 0.018
      real,parameter  :: v_charnock     = 1.5e-5 !(m2/s)
      real,parameter  :: rz0            = 2.0
      real,dimension(:,:),allocatable :: sr_mix
      integer :: i, j

      real :: buoy, tstv, obuk, ra, y0, yra
      real :: xland, sr_land, sr_help


      allocate(sr_mix(i1:i2, j1:j2))
      do j=j1,j2
         do i=i1,i2
            
            xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.

            ! SEA:
            ! surface roughness  from Charnock equation
            ! friction velocity from surface stress
            !
            if(xland < 0.99) then 
               ustar_sea(i,j)=sqrt(sstr(i,j))
               sr_sea(i,j)=alfa_charnock1*v_charnock/ustar_sea(i,j) +  &
                           alfa_charnock2*sstr(i,j)/grav 
            else 
               ustar_sea(i,j)=0.0
               sr_sea(i,j)=0.0
            endif

            !
            ! LAND
            ! calculate the 'local' surface roughness for momentum that 
            ! is consistent with 10 m wind and the Olsson data base, 
            ! we assume that the windspeed at 75 m is independent of
            ! surface roughness
            !
            if ( xland > 0. ) then
               !>>> TvN
               ! Over land, the friction velocity ustar 
               ! is calculated from the 10 m wind speed, u10.
               ! In IFS u10 is a diagnostic variable, calculated to be
               ! compatible with "open-terrain" wind speed observations.
               ! Over rough or inhomogeneous terrain (z0M > 0.03 m),
               ! it is calculated using an aerodynamic roughness
               ! length typical for open terrain with grassland (0.03 m)
               ! (see IFS cy31r1 documentation, p. 46).

               ! In the expressions applied in TM5,
               ! the stability profile functions Psi_M have disappeared,
               ! and only the logarithmic terms are kept.
               ! Apart from this, the expression for ustar was incorrect:
               ! 10./sr_land should be 75./sr_land.
               ! This has now been corrected.

               ! Moreover, over islands and near coast lines, 
               ! zeros can occur in sr_ols_dat,
               ! which have to be removed.
               ! However, a lower bound of 1e-2 m = 1 cm
               ! seems too high for smooth surfaces,
               ! like deserts and ice caps.
               ! The minimum positive value in sr_ols_dat
               ! is 2.5e-4 m = 0.025 cm,
               ! which is obtained in parts of Antarctica.
               ! This is likely the result of regridding of
               ! a field with minimum value of 0.1 cm
               ! from 0.5x0.5 to 1x1 degrees.

               ! Please note that with the introduction of CY31R1,
               ! the orographic contribution to the aerodynamic roughness length
               ! in IFS has been replaced by an explicit specification of stress
               ! on model levels due to turbulent orographic form drag,
               ! and the climatological variable previously used
               ! has been replaced by a prognostic variable.
               ! In TM5, the climatological variable is still used,
               ! but it would be better to use the prognostic variable instead.

               ! Note also that the same calculation 
               ! is done in diffustion.F90.

               !sr_land=max(1e-2,sr_ols_dat(iglbsfc)%data(i,j,1))  ! occurs at Islands, etc.
               sr_land=max(1e-3,sr_ols_dat(iglbsfc)%data(i,j,1))
               sr_help=min(sr_ecm_dat(iglbsfc)%data(i,j,1),0.03)
               !fd ustar_land=vKarman*wind10m(i,j)/alog(10./sr_help)
               ! !ustar consistent with 'clipped' large scale roughness
               !ustar_land(i,j)=vKarman*wind10m(i,j)/alog(10./sr_help)*&
               !                alog(75./sr_help)/alog(10./sr_land)
               ustar_land(i,j)=vKarman*wind10m(i,j)/alog(10./sr_help)*&
                               alog(75./sr_help)/alog(75./sr_land)
               ! <<< TvN
            else
               ustar_land(i,j)=0.0
               sr_land=0.
            endif
            ustar_loc(i,j)=xland*ustar_land(i,j)+(1-xland)*ustar_sea(i,j)
            sr_mix(i,j) =xland*sr_land  +(1-xland)*sr_sea(i,j)

         end do !i
      end do !j

      !
      ! Calculation of quasi laminar resistance of the layer above the surface.
      ! E.g. Walcek, Atmos. Environ, 20, pp. 949-964, 1986, and earlier refs.
      !     
      do j=j1,j2
         do i=i1,i2
            !>>> TvN
            ! Instead of the current approach,
            ! it is proposed to use the surface roughness length 
            ! for heat directly as it is calculated prognostically 
            ! in IFS (GRIB number 245).
            !<<< TvN
            rb(i,j)=(rz0/(ustar_loc(i,j)*vkarman)) 
         end do
      end do

      !
      ! calculate the aerodynamic resistance
      !
      ! ra,the aerodynamic resistence in s/m over a layer with height z is
      ! the integral from z0 to z of 1/K(z) dz with K(z)=k.u*.z/phi(h)
      ! phi(h)=0.74(1-9z/l)**-0.5 for l<0 and 0.74+4.7(z/l) for l>0.
      ! by integration we get ra=0.74/(k.u*).(ln(z0/z)+ghi(z0)-ghi(z))
      ! ghi(z)=-6.4(z/l) for l>0 and 2.ln(y(z)+1) for l<0 with 
      ! y(z)=(1-9.z/l)**0.5.
 
      do j=j1,j2
         do i=i1,i2
            buoy = -sshf_dat(iglbsfc)%data(i,j,1) / rhoCp &
                   -0.61 * t2m_dat(iglbsfc)%data(i,j,1) * slhf_dat(iglbsfc)%data(i,j,1)/rhoLv
            tstv=-buoy/ustar_loc(i,j) 
            obuk=1e-6
            if( abs(tstv) .gt. tiny(tstv) ) then 
               obuk = ustar_loc(i,j)*ustar_loc(i,j)*t2m_dat(iglbsfc)%data(i,j,1)/(tstv*grav*vKarman)
            end if

            if( obuk > 0. ) then !   stable conditions
               ra = 0.74*( alog(Href/sr_mix(i,j)) + 6.4*(Href-sr_mix(i,j))/obuk )&
                    / (vKarman*ustar_loc(i,j)) 
            else              !   unstable
               y0 = sqrt(1.-9.*sr_mix(i,j)/obuk)+1.       
               yra = sqrt(1.-9.*Href/obuk)+1. 
               ra = 0.74*(alog(Href/sr_mix(i,j))+2.*(alog(y0)-alog(yra)))/ &
                    (vKarman*ustar_loc(i,j)) 
            end if
            raero(i,j)=max(10.,min(ra,1e10))  !
         end do !i
      end do !j

      if ( okdebug ) then
         call dist_arr_stat(dgrid(iglbsfc), 'ustar_loc', ustar_loc, 0, status )
         IF_NOTOK_RETURN(status=1)
         call dist_arr_stat(dgrid(iglbsfc), 'raero',     raero    , 0, status )
         IF_NOTOK_RETURN(status=1)
         call dist_arr_stat(dgrid(iglbsfc), 'rb',        rb       , 0, status )
         IF_NOTOK_RETURN(status=1)
         write(*,*) 'end calc_aero_ustar'
      end if

      deallocate(sr_mix)

    end subroutine dd_calc_ustar_raero_rb
    !
    ! 
    subroutine dd_calc_rstom_rahcan
      ! -----------------------------
      ! Purpose
      ! ---------
      ! Calculate water vapour stomatal resistance from the PAR
      ! (Photosythetically Active Radiation) which is 0.55 times
      ! the net short wave radiation (ssr) at the surface.
      ! Calculate rahcan from ustar and lai
      ! 
      ! External
      ! --------
      ! none
      !
      ! Reference
      ! ----------
      ! none
      !--------------------------------------------------
      !
      !  -- constants used for the calculation of the stomatal resistance 
      !     (see ECHAM3 manual)
      implicit none
      real,parameter  :: a=5000.
      real,parameter  :: b=10.
      real,parameter  :: c=100.
      real,parameter  :: vk=0.9
      real,parameter  :: vlt=1.    
      !
      real, parameter :: foresth=20. 
      integer         :: j, i
      real            :: vpar, d
      real            :: canht,laihelp

      !
      !  -- recalculated rstom for a LAI of 1 (see ECHAM3 manual)
      !
      do j=j1,j2
         do i=i1,i2
            !
            !    --  calculation of PAR from net short wave radiation
            !    --  radiation is corrected for daily cycle
            !
            rstom(i,j)=1e10    !high resistance during the night
            if (ssr_dat(iglbsfc)%data(i,j,1) > 0.) then
               vpar=max(1.,0.55*ssr_dat(iglbsfc)%data(i,j,1))
               d=(a+b*c)/(c*vpar) 
               rstom(i,j)=(vk*c)/((b/(d*vpar))*alog((d*exp(vk*vlt)+1.)/(d+1.))-&
                    alog((d+exp(-vk*vlt))/(d+1.)))
            end if !ssr > 0.
         end do !i
      end do !j
      !
      ! -- computation of in-canopy aerodynamic resistance from canopy height
      !    and the friction velocity and the Leaf Area Index 
      !    (see Erisman & Van Pul)
      !
      do j=j1,j2
         do i=i1,i2
            !
            !
            ! -- calculation of canopy height CANHT from effective fraction 
            !    of high vegetation assuming that forest has a canopy height 
            !    of 20 m (other vegetation 0 m)
            !    
            canht= vgrat_high(i,j)*foresth
            ! laihelp: the maximum values derived from ECMWF are more realistic    
            laihelp=min(lai1(i,j),lai(i,j))
            rahcan(i,j)=max(1.e-5,14.*laihelp*canht/ustar_loc(i,j))
         end do !j
      end do !i

      !cmk   if ( okdebug ) &
      !cmk    call dumpfield(0,'rahcan'//c_time//'.hdf',rahcan,'rahcan')
      !cmk   if ( okdebug ) call dumpfield(0,'rstom'//c_time//'.hdf',rstom,'rstom')
      if ( okdebug ) write(*,*) 'dd_calc_rstom_rahcan: end'
      !
    end subroutine dd_calc_rstom_rahcan
    !
    !
    subroutine dd_calc_inisurf
      ! ------------------------
      ! Purpose
      !---------
      ! Calculate some surface fields later needed in the calculation
      ! 
      ! External
      ! --------
      ! none
      !
      ! Reference
      ! ----------
      ! none
      !--------------------------------------------------
      use Binas, only: pi
      implicit none
      real,parameter    :: sncr     = 0.015       ! critical snow cover depth
      real,parameter    :: tstar    = 273.
      real,parameter    :: wlmax    = 2.e-4              
      real,parameter    :: ewsmx    = 0.323
      real,parameter    :: ewsmx_sat= 0.472
      real,parameter    :: wpwp     = 0.171
 !FD23012004
      real,parameter    ::  eff=0.5            ! parameters needed for sea/bubble formation
      real,parameter    ::  rdrop=0.005        ! cm
      real,parameter    ::  zdrop=10.0         ! cm
      real,parameter    ::  eps=0.6            ! parameter related to bubble formation
      real              ::  qdrop,s,phi,alpha1,vk1,vk2,alpharat   ! auxiliary parameters for bubble formation
      real              ::  sr_help            ! surface roughnes consistent with 10 m wind.
  !FD23012004

      !
      ! for albedo calculations
      real              :: snow, freesea, bare, desert
      real              :: e,es,wcr,wlmx,xland,ahelp,vgr,laihelp
      integer           :: i,j    ! auxiliary variables
      
      ! --- begin ---------------------------------------

      do j=j1,j2
         do i=i1,i2
            xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.0       ! land-sea fraction directly from ECMWF
            ! the maximum values derived from ECMWF are more realistic
            laihelp=min(lai1(i,j),lai(i,j))

            !
            ! -- calculation of monthly snow cover fraction sd using 
            !    constant critical snow depth
            !    alternatively this critical snow depth may be chosen as 
            !    a linear function of ln(sr_ecm)
            !    B. v.d. Hurk (2002) suggests 0.015 m (=sncr) for z0<0.25
            !    and 1 m for z0>5m and log-linear in between in between.
            !    factor 0.9 is introduced to account for the fact that 
            !    it is unlikely 
            !    that 'high' vegetation is completely snow covered
            ! 
            ! FD25.03.2003
            snow_cover(i,j)=min(xland,sd_dat(iglbsfc)%data(i,j,1)/sncr)

            !
            ! -- calculation of wet skin fraction from the skin reservoir 
            !    content src (prognostic variable in ECMWF)
            !    the vegetation fractions and their attributed LAI, 
            !    Wlmax which is the maximum amount of water that can be held 
            !    on one layer of leaves or bare ground, Wlmx is the 
            !    maximum skin reservoir content
            !   
            if ( xland > 0.0 ) then
               ! bare soil fraction small discrepancy when lakes are present
               bare=max(0.,1.-vgrat_high(i,j)-vgrat_low(i,j))
               wlmx=wlmax*(bare+laihelp)             
               wet_skin(i,j)=min(1.,src_dat(iglbsfc)%data(i,j,1)/wlmx)
            else
               wet_skin(i,j)=1.0   !for sea  CMK bug 01-07-2003
            end if
            !
            ! -- calculation of water stress factor FWS with Wsmx is the 
            !    field capacity.
            !
            !    Field capacity is defined as the maximum amount of water 
            !    that an column of soil can hold
            !    against gravity 24-48 hours after wetting of the soil
            !    Wcr is a critical value, Wpwp is the permanent wilting point 
            !    and Ws is the total amount of water available in the 
            !    root zone (ECHAM3 manual)
            !
            fws(i,j)=1e-5
            if (xland > 0) then
               wcr=0.5*ewsmx_sat
               ! used ewsmx_sat instead of ewsmx, is more consistent 
               ! with values of swvl1    
               if ( swvl1_dat(iglbsfc)%data(i,j,1) > wpwp &
                        .and. swvl1_dat(iglbsfc)%data(i,j,1) < wcr ) &
                    fws(i,j)=(swvl1_dat(iglbsfc)%data(i,j,1)-wpwp)/(wcr-wpwp)
               if ( swvl1_dat(iglbsfc)%data(i,j,1) > wcr ) fws(i,j)=1.
            end if
            !
            ! -- Computation of relative humidity (2 m) out of the air and dew 
            !    temperature at 2 m  which is used 
            !    
            es=exp(19.59*(t2m_dat(iglbsfc)%data(i,j,1)-tstar)/t2m_dat(iglbsfc)%data(i,j,1))
            e=exp(19.59*(1.-tstar/d2m_dat(iglbsfc)%data(i,j,1)))
            rh2m(i,j)=e/es
	    !
	    !-- Calculation of H2 deposition over arable field 
	    !   Sanderson et al., J. Atmos. Chem. 2000
	    !   Yonemura et al., JGR 2000
	    !   units: m/sec
	    !
	    ddepvel_h2(i,j) = min(10.e-4, max(1.e-7,(-41.39 * swvl1_dat(iglbsfc)%data(i,j,1) + 16.85) *1e-4)) 
	    
            !
            ! calculate albedo (needed for photolysis rates) from 
            ! different fractions
            ! 
            ! vgr: coverage by low and high vegetation
            vgr= vgrat_low(i,j)+vgrat_high(i,j)
            !
            ! --  if high vegetation is present and snow this may alter 
            !     the effective snow albedo: let high vegetation prevail
            snow = snow_cover(i,j)- max(0.0,snow_cover(i,j)+vgrat_high(i,j)-1.0) 

            bare = max(0.,1.0 - vgr - snow)
            ! soils with pH values larger than 7.3 
            desert = soilph(i,j,3) + soilph(i,j,4)
            ! when more desert is present than bare land...
            desert = desert - max(0.0,desert-bare)  
            bare = bare -desert 
            freesea = max(0.,1.-xland-ci_dat(iglbsfc)%data(i,j,1))
            !
#ifndef without_photolysis
            !AJS: uv-vis albedo computed from land use;
            !AJS: in future, ecmwf field should be used:
            !AJS:   name   :  UV visible albedo for direct radiation
            !AJS:   abbrev :  ALUVP  ;     unit   :  0-1
            !AJS:   code   :  15     ;     table  :  128
            ags(i,j) = freesea*0.05 + ci_dat(iglbsfc)%data(i,j,1)*0.70 + &
                 xland*(desert*0.10+bare*0.05+vgr*0.01+snow*0.70)   

            ! AS: make sure that it does not exceed 0.7 in cases where
            !     sea ice and lsmask are not fully consistent
            ags(i,j) = min(0.7,ags(i,j))
#endif

            if (freesea>0.01) then   !

               ! bubble bursting effect,see Hummelshoj, equation 10
               ! relationship by Wu (1988), note that Hummelshoj has not 
               ! considered the cunningham factor which yields a different 
               ! vb curve, with smaller values for small particles
               ! according to LG indeed 10 m windspeed (instead of 1 m windspeed in use in ECHAM5

               alpha(i,j)=MIN(1.0,MAX(1.e-10,1.7e-6*wind10m(i,j)**3.75))   ! set maximum!
               qdrop=5.*(100.*alpha(i,j))   !  100 is the flux of bubbles per cm^2/s (see Monohan(1988))

               bubble(i,j)=((100.*ustar_sea(i,j))**2.)/(100.*wind10m(i,j)) +  & 
                             eff*(2.*pi*rdrop**2.)*(2.*zdrop)*(qdrop/alpha(i,j))

               !--- Correction of particle radius for particle growth close to the 
               !    surface according to Fitzgerald, 1975, the relative humidity over 
               !    the ocean is restricted to 98% (0.98) due to the salinity

               s=MIN(0.98,rh2m(i,j))

              !fd23012004 beta(i,j)=EXP((0.00077*s)/(1.009-s)) THIS term was present in ECHAM code !
              !; max. value reached for this parameter is 1.04 and is ignored here.

              phi=1.058-((0.0155*(s-0.97))/(1.02-s**1.4))
              alpha1=1.2*EXP((0.066*s)/(phi-s))   
              vk1=10.2-23.7*s+14.5*s**2. 
              vk2=-6.7+15.5*s-9.2*s**2. 
              alpharat=1.-vk1*(1.-eps)-vk2*(1.-eps**2)
              alphae(i,j)=alpharat*alpha1

              !--- Over land no correction for large humidity close to the surface:

            else ! land surface
              alpha(i,j)=0.
              bubble(i,j)=0.
              alphae(i,j)=1.
            endif

         end do
      end do

      if ( okdebug ) write(*,*) 'after dd_cal_inisurf'

    end subroutine dd_calc_inisurf

    !--------------------
    !
    ! This subroutine calculates the surface resistance as
    ! a function of a series of resistances
    !
    !  purpose
    !  -------
    !  determine surface resistance "rsurf" for dry deposition 
    !  calculations
    !
    !  external
    !  --------
    !  
    !  reference
    !  ---------
    !  Ganzeveld and Lelieveld (1996) and references therein
    !
    !------------------------------------------------------------------
    !
    !  -- resistances and auxiliary variables
    !
    !------------------------------------------------------------------

    subroutine dd_calc_rs

      use binas     , only : vKarman, pi
      use chem_param, only : density_ref
      use chem_param, only : ndep, idep
      use chem_param, only : ddep_diffrb, ddep_rsoil, ddep_rwat, ddep_rws
      use chem_param, only : ddep_rsnow , ddep_rmes , ddep_rcut, ddep_diffcf
      use chem_param, only : diffcf_o3
      use chem_param, only : iso2, iso4, inh3, ihno3, ino, ino2, io3, ico
      use chem_param, only : nrdep, lur
      use toolbox   , only : coarsen_emission
      use ParTools,   only : isRoot

      integer,parameter  :: avg_field = 1

      real, parameter :: dynvisc=1.789e-4*2. ! g cm-1 s-1 CHECKED FD OK, there is temp. dependence Perry p. 3-248.
                                             ! unit is g cm-1 s+1 FD 
                                             ! checked with Seinfeld---> factor 2. came out (diameter ? radius)
      real, parameter :: cl=0.066*1e-4       ! mean free path [cm] (particle size also in cm)
      real, parameter :: bc= 1.38e-16        ! boltzman constant [g cm-2 s-1 K-1] (1.38e-23 J deg-1) =>binas
      real, parameter :: kappa=1.            ! shapefactor
      real, parameter :: visc=0.15           ! KINEMATIC molecular viscocity [cm2 s-1] 
                                             ! this is also function of temperature FD

      real                             :: xland1,xland2,xland3,xland4
      real                             :: xwat1,xwat2,xsum
      real                             :: rstomx,rsveg,ra1
      real                             :: vdland1,vdland2,vdland3,vdland4
      real                             :: rssoil,rsws,vdwat1,vdwat2
      real                             :: rssnow,rveg,rleaf,rcanopy
      real                             :: w10,zrsnhno3,ustarh
      real                             :: cvsh,cvwh,evgrath,tsoilph
      real                             :: xland,ts1,laihelp
      real,dimension(:,:), allocatable :: rwatso4
      real,dimension(:,:), allocatable :: rsoilnh3,rsoilso2,rsoilso4
      real,dimension(:,:), allocatable :: rsnowhno3,rsnowso2
      !
      integer :: jt,j,i,jdep, irdep
      character(len=8) :: adum

      real                     :: zr, vd_land, vd_sea, um, dc, cunning, relax, ust_land, st_land
      real                     :: sc,  vb_land, vi_land,  ust_sea, st_sea, re
      real                     :: vbsea, visea, vkdaccsea, vkd_sea, vkd_land
      real                     :: densaer

      real, allocatable   ::  curr_diffrb(:), curr_rsoil(:), curr_rwat(:), curr_rws(:)
      real, allocatable   ::  curr_rsnow (:), curr_rmes(:) , curr_rcut(:), curr_diffcf(:)
      ! Openmp parameters
      integer             :: js, je


      ! --- begin ---------------------------
      
      !
      ! *** deposition
      !
      
      if ( ndep > 0 ) then

        allocate(rwatso4  (i1:i2, j1:j2))
        allocate(rsoilnh3 (i1:i2, j1:j2))
        allocate(rsoilso2 (i1:i2, j1:j2))
        allocate(rsoilso4 (i1:i2, j1:j2))
        allocate(rsnowhno3(i1:i2, j1:j2))
        allocate(rsnowso2 (i1:i2, j1:j2))

        allocate( curr_diffrb(ndep) )
        allocate( curr_rsoil (ndep) )
        allocate( curr_rwat  (ndep) )
        allocate( curr_rws   (ndep) )
        allocate( curr_rsnow (ndep) )
        allocate( curr_rmes  (ndep) )
        allocate( curr_rcut  (ndep) )
        allocate( curr_diffcf(ndep) )


        !
        ! -- extension with the Wesely scheme, 1989 
        !    (Atmospheric Environ.,1989)
        !    in which the resistances of trace gases are estimated from 
        !    the reactivity relative to 
        !    ozone and the dissolvation relativeto SO2. The deposition process of 
        !    these two trace is used as a reference for the estimation.
        !
        ! calculate some tracer specific resistances
        !
        curr_diffrb = ddep_diffrb
        curr_rsoil  = ddep_rsoil  
        curr_rwat   = ddep_rwat  
        curr_rws    = ddep_rws   
        curr_rsnow  = ddep_rsnow 
        curr_rmes   = ddep_rmes  
        curr_rcut   = ddep_rcut  
        curr_diffcf = ddep_diffcf
        !
        do j=j1,j2
           do i=i1,i2
              ustarh=ustar_loc(i,j)

              ! -- calculation of the integrated resistance from the mass size 
              !    distribution for rural continental aerosol with an unimodal distr.
              !    with the mean mass radius at about 0.4 um. 
              !    (observations by Mehlmann (1986)
              !    as referenced in Warneck, 1988) and the friction velocity. 
              !    Polynominal fit!
              !    Over land the surface resistance of bare soil and snow 
              !    covered surfaces is calculated
              !
              ! -- following distributions are used:
              ! 1. deposition over land using  a rural continental size distribution
              !    i.e. mostly accumulation range aerosol
              ! 2. deposition over water using a marine size distribution 
              !    i.e. bimodal distribution with a 30% coarse fraction
              !    fjd this is most relevant when explicit chemistry in seasalt 
              !    is considered
              ! 3. deposition over water using rural continental size distribution
              !    the deposition of 'seasalt' sulfate may be calculated using
              !    the relation vd2=0.7*vd1+0.3*vdsssulfate
              !
              ! 1st distribution:
              rsoilso4(i,j) = max( 1., &
                   100./(0.08-0.57*ustarh+1.66*ustarh**2-0.36*ustarh**3) )
              ! 2nd distribution:
              ! rwatso4(i,j) =
              ! max(1.,100./(0.12-0.24*ustarh+5.25*ustarh**2 -1.84*ustarh**3))
              ! 3rd distribution: 
              rwatso4(i,j)  = max( 1., &
                   100./(0.07-0.1*ustarh+2.1*ustarh**2 -0.20*ustarh**3) )
              !   
              ! -- parameterization of snow/ice SO2 resistance as a function of the 
              !    snow/ice temperature, based on observations by Dasch and Cadle
              !
              ts1=max(200.,t2m_dat(iglbsfc)%data(i,j,1))
              !
              ! FD25032003 the measurements are not completely consistent; 
              ! put a lower
              ! limit of vd=0.10 cm/s for SO2 to avoid excessive accumulation
              ! if snow is melting high uptake
              !
              rsnowso2(i,j)=amin1(max(10.,10.**(-0.09*(ts1-273.)+2.4)),1.e+3)
              ! snow is melting changed (advise Wouter Greull) FD072003
              if ( ts1 > (273.+3.) .and. sshf_dat(iglbsfc)%data(i,j,1) > 10. ) rsnowso2(i,j)=50.

              !
              ! arbitrary attribution of soil resistance to NH3 deposition
              ! acid soils have lowest deposition
              !
              tsoilph=sum(soilph(i,j,1:5))
              rsoilnh3(i,j)=1e10
              if (tsoilph > 0.) &
                   rsoilnh3(i,j)=  max(25.,soilph(i,j,1)*25. +soilph(i,j,2)*25.&
                   +soilph(i,j,3)*200.+soilph(i,j,4)*200.+soilph(i,j,5)*100.)
              !
              ! -- Computation of rsoil for SO2 from soil pH three different 
              !    rh classes
              !    The rsoil values for each class are determined out of 
              !    regression and the average value of pH of the class. 
              !    Concerning rh, a wet, dry and very dry class 
              !    is discerned with the threshold value of 60% and 40%
              ! 
              rsoilso2(i,j)=1e10
              !
              ! -- for rh > 60 %
              !
              if ( tsoilph > 0. ) rsoilso2(i,j)=max(25.,(soilph(i,j,1)*115.+&
                   soilph(i,j,2)*65+soilph(i,j,3)*25.+ &
                   soilph(i,j,4)*25.+soilph(i,j,5)*70.)/tsoilph+&
                   amin1(max(0.,1000.*exp(269.-ts1)),1.e+5))

              ! -- for rh less than 60% and larger than 40 %

              if ( rh2m(i,j) < 0.6 .and. rh2m(i,j) > 0.4) &
                   !cmk    rsoilso2(i,j)=max(50.,(rsoilso2(i,j)*3.41-85.)+
                   !cmk amin1(max(0.,1000.*exp(269.-ts1)),1.e+5))
                   rsoilso2(i,j)=max(50.,rsoilso2(i,j)*3.41-85.)+ &
                   amin1(max(0.,1000.*exp(269.-ts1)),1.e+5)

              ! -- for rh less than 40 %

              if (rh2m(i,j).le.0.4) &
                   !cmk    rsoilso2(i,j)=max(50.,amin1(1000.,(rsoilso2(i,j)*3.41-&
                   !cmk    85.+((0.4-rh2m(i,j))/0.4)*1.e+5)+&
                   rsoilso2(i,j)=max(50.,amin1(1000.,(rsoilso2(i,j)*3.41-85.+ &
                   ((0.4-rh2m(i,j))/0.4)*1.e+3)+& 
                   ! 1e3, see paper Laurens, JGR
                   amin1(max(0.,1000.*exp(269.-ts1)),1.e+5)))

              ! -- Temperature correction term for HNO3 above snow surface and ice
              !    (Wesely, 1989), recomputation from K to 0C

              rsnowhno3(i,j)=amin1(max(10.,1000.*exp(269.-ts1)),1.e+5)   
              ! snow is melting changed (advise Wouter Greull) FD072003
              if ( ts1 > (273.+3.) .and. sshf_dat(iglbsfc)%data(i,j,1) > 10. ) rsnowhno3(i,j)=50. 
              !
              !
           end do !nlon
        end do !nlat


        !
        ! ***
        !

        vd11(:,:,:) = 1e-10
        !CMK rsurf(:,:,:)=1e+5 
        !
        !  -- Loop for tracers, all timesteps
        !
        do jdep=1,ndep
           jt = idep(jdep)
           !
           !  -- Only if a value for DIFFCF (diffusivity) is defined, the program 
           !   calculates a surface resistance
           !
           if ( curr_diffcf(jdep) > 1e-5 ) then
              !  
              do j=j1,j2
                 do i=i1,i2

                    xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.0
                    ! laihelp: the max values derived from ECMWF are more realistic
                    laihelp=min(lai(i,j),lai1(i,j))
                    !   -- Assigning of values calculated in previous routine
                    !   rsnow/rsoil of other components in data statement
                    if (jt == iso2) then
                       curr_rsnow(jdep)=rsnowso2(i,j)
                       curr_rsoil(jdep)=rsoilso2(i,j)
                    end if
                    if (jt == inh3) curr_rsoil(jdep)=rsoilnh3(i,j)
                    if (jt == iso4) then
                       !   set snow resistance of so4 equal to soil resistance
                       curr_rsnow(jdep)=rsoilso4(i,j) 
                       curr_rsoil(jdep)=rsoilso4(i,j)
                       curr_rwat(jdep)=rwatso4(i,j)
                    end if
                    !
                    if (jt == ihno3) curr_rsnow(jdep)=rsnowhno3(i,j)

                    !   -- Computation of stomatal resistance for component X, 
                    !      correction based on differences in molecular
                    !      diffusion of H2O and X and the water
                    !      stress is also considered, FWS

                    rstomx=curr_diffcf(jdep)*rstom(i,j)/fws(i,j)

                    !   -- Calculation of mesophyll resistance of NO 
                    !      and NO2 as function of
                    !      the stomatal resistance of ozone

                    if (jt == ino ) curr_rmes(jdep)= 5.*diffcf_o3*rstom(i,j)/fws(i,j)
                    if (jt == ino2) curr_rmes(jdep)=0.5*diffcf_o3*rstom(i,j)/fws(i,j)

                    rleaf=(1./((1./curr_rcut(jdep))+(1./(rstomx+curr_rmes(jdep)))))
                    ! linear upscaling from leaf scale to canopy scale 
                    ! applying the LAI derived from Olson 

                    rcanopy=rleaf/max(1e-5,laihelp)
                    rveg=(1./((1./(rahcan(i,j)+rb(i,j)*curr_diffrb(jdep)+curr_rsoil(jdep)))+ &
                         (1./rcanopy)))

		    !  -- Exception for CO: Use Sanderson et al. parameterization
		    !     CO dry dep. scaled from H2 dry dep., using soil wetness.
		    !     Note: This soil resistance is critical one, so one may
		    !     neglect other terms (cuticle, boundary layer, ...)
		    if (jt == ico) rveg = 1./(ddepvel_h2(i,j) * 0.6 ) 


                    !-- It can be assumed that HNO3 deposition only depends on ra
                    !   so that surface resistance = 0, however definition of 
                    !   minimal surface resistance in order to avoid 
                    !   extreme Vd values.

                    if (jt == ihno3) rveg=1.

                    !-- Computation of surface resistance Rs (s m-1) above land
                    !   Vd for surface with vegetation and already incorporated 
                    !   is the vegetation boudanry layer resistance

                    rsveg=rveg+rb(i,j)*curr_diffrb(jdep)

                    !-- Sulfate deposition calculation over vegetation according 
                    !   to the observations by Wesely (1985), 
                    !   see the Dry Deposition Inferential Model

                    if (jt == iso4) then 
                       !
                       ! fd removed the factor stheta 
                       ! (standard deviation of the wind direction)

                       if ( ssr_dat(iglbsfc)%data(i,j,1) > 250 ) then
                          rsveg=1./(0.01*ustar_loc(i,j))
                       else
                          rsveg=1./(0.002*ustar_loc(i,j))
                       end if
                       !
                       !   -- for particle deposition the rb is already 
                       !      implicitly included
                       rssoil=curr_rsoil(jdep)
                       !
                       !   assumed that the wet skin fraction consists of vegetation
                       !
                       rsws=rsveg 
                       rssnow=curr_rsnow(jdep)
                       !
                    else !jt neq iso4
                       ! 
                       !-- Rs for surface with bare soil
                       !
                       rssoil=curr_rsoil(jdep)+rb(i,j)*curr_diffrb(jdep)
                       !
                       !-- Wet skin fraction, 
                       !   It is assumed that the wet skin fraction 
                       !   mainly consists of vegetation thus laminar 
                       !   resistance for vegetation may be applied
                       !
                       rsws=curr_rws(jdep)+rb(i,j)*curr_diffrb(jdep)
                       !
                       !-- Snow coverage
                       !
                       rssnow=curr_rsnow(jdep)+rb(i,j)*curr_diffrb(jdep)
                       !
                    end if
                    !
                    !-- land surface deposition velocity
                    !

                    cvsh=snow_cover(i,j) - &
                         max(0.0,snow_cover(i,j)+vgrat_high(i,j)-1.0)
                    ! fd same assumption as photolysis

                    !cmk: need to correct the rsoil resistance here.
                    !     fraction will become bare soil, which in the 
                    !     vegetation deposition calculation
                    !     will result in a too low resistance, 
                    !     since the 'soil' part will be snow covered!

                    cvwh=min(xland,wet_skin(i,j) )
                    evgrath=min(xland,vgrat_low(i,j)+vgrat_high(i,j))

                    ra1=raero(i,j)
                    ! snow covered land fraction
                    xland1=cvsh
                    vdland1=1./(rssnow+ra1)
                    ! wet skin covered land fraction (not covered by snow)
                    xland2=max(0.,cvwh-xland1)
                    vdland2=1./(rsws+ra1)
                    ! vegetation covered land fraction 
                    ! (not covered by snow and wet skin)
                    xland3=max(0.,evgrath-xland1-xland2)
                    vdland3=1./(rsveg+ra1)
                    ! bare soil covered land fraction   
                    ! (landfraction not covered by snow, wet skin, or vegetation)
                    xland4=max(0.,xland-xland1-xland2-xland3)
                    vdland4=1./(rssoil+ra1)  

                    xwat1=min(1-xland,ci_dat(iglbsfc)%data(i,j,1))
                    vdwat1=1./(curr_rsnow(jdep)+rb(i,j)*curr_diffrb(jdep)+ra1)
                    xwat2 = max(0.,1.-xland-ci_dat(iglbsfc)%data(i,j,1))
                    vdwat2=1./(curr_rwat(jdep)+rb(i,j)*curr_diffrb(jdep)+ra1)
                    if (jt == iso4) then
                       vdwat1=1./(curr_rsnow(jdep)+ra1)
                       vdwat2=1./(curr_rwat(jdep)+ra1)
                    end if
                    xsum=xland1+xland2+xland3+xland4+xwat1+xwat2
                    !fd    if (xsum.ne.1.) then
                    !fd    print *,'sum',i,j,xsum,xland,xland1,xland2,'xl3',&
                    !fd    xland3,xland4,xwat1,xwat2,cvsh,cvwh,evgrath,ci_dat(iglbsfc)%data(i,j,1)
                    !fd    end if
                    !
                    !-- Computation of deposition velocity 
                    !
                    vd11(i,j,jdep)=xland1*vdland1+xland2*vdland2+xland3*vdland3+&
                         xland4*vdland4+xwat1*vdwat1+xwat2*vdwat2
                 end do !End of loop longitude
              end do  !End of loop latitude

              adum=names(jt)
              i=len_trim(adum)

              if ( okdebug ) then
                 call dist_arr_stat(dgrid(iglbsfc), 'vd'//adum(1:i), vd11(:,:,jdep), 0, status )
                 IF_NOTOK_RETURN(status=1)
              end if
              

           end if !  (curr_diffcf(jdep) > 1e-5) 
           !
        end do   !End of loop tracer

        deallocate(rwatso4)
        deallocate(rsoilnh3)
        deallocate(rsoilso2)
        deallocate(rsoilso4)
        deallocate(rsnowhno3)
        deallocate(rsnowso2)

        deallocate( curr_diffrb )
        deallocate( curr_rsoil  )
        deallocate( curr_rwat   )
        deallocate( curr_rws    )
        deallocate( curr_rsnow  )
        deallocate( curr_rmes   )
        deallocate( curr_rcut   )
        deallocate( curr_diffcf )

      end if  ! ndep > 0

      !
      ! *** aerosols
      !
      
      if ( nrdep > 0 ) then

        ! look up table for different aerosol radii:
        do irdep = 1, nrdep

!PLS !$OMP PARALLEL &
!PLS !$OMP  default (none) &
!PLS !$OMP  shared  (irdep, lur, wind10m, lsmask_dat, alphae, t2m_dat, &
!PLS !$OMP           ustar_land, raero, ustar_sea, sr_sea, alpha, bubble, &
!PLS !$OMP           vd11a) &
!PLS !$OMP  private (i, j, js, je, zr, vd_land, vd_sea, um, xland, cunning, &
!PLS !$OMP           dc, densaer, relax, sc, ust_land, st_land, vb_land, &
!PLS !$OMP           vkd_land, ust_sea, st_sea, re, vbsea, visea, &
!PLS !$OMP           vkdaccsea, vkd_sea, vi_land)
     !      call my_omp_domain (1, nlat180, js, je)
           do j=j1,j2
              do i=i1,i2

             zr = 2.0e-4 * lur(irdep)   ! diameter in cm !
             vd_land=0.
             vd_sea=0.

             um=wind10m(i,j)
             xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.0   ![]fraction

             !--- Cunningham factor: 

             cunning=1.+(cl/(alphae(i,j)*zr))*(2.514+0.800*EXP(-0.55*zr/cl))

             !--- Diffusivity:

             dc=(bc* t2m_dat(iglbsfc)%data(i,j,1)*cunning)/(3.*pi*dynvisc*(alphae(i,j)*zr)) ! [cm2/s] FD

             ! Relaxation:
             ! relax represents characteristic relaxation time scale [Seinfeld, p. 319]
             !      [     kg m-3 => g cm-3    ]
             densaer = density_ref   ! reference density (e.q. 1800 kg/m3)
             relax=cunning*densaer*1.E-3*((alphae(i,j)*zr)**2. ) &
                   /(18.*dynvisc*kappa)


             ! Sedimentation is calculated operator split in the subroutine sedimentation:
             !
             ! sedspeed=(((((alphae(i,j)*zr(i,j)))**2.)* &
             !          densaer(jl,klev,jmod,jrow)*1.E-3*grav*cunning)/(18.*dynvisc)) ! note grav should be in cm
             !    [     kg m-3 => g cm-3    ]

             ! Calculation of schmidt 

             sc      =visc/dc  ![cm2 s-1]/[cm2 s-1] dimensionless

             if (xland.gt.0.01) then
                !    note that in ECHAM there is a difference between vegetaton and snow/bare soil
                !    in TM5 there we assume this is incorporated in the 1x1 fields

                ust_land=ustar_land(i,j)

                !
                !    calculation of stokes numbers

                st_land  =max(0.1,(relax*(100.*ust_land)**2.)/visc)

                !--- Calculation of the dry deposition velocity
                !    See paper slinn and slinn, 1980, vd is related to d**2/3 
                !    over land, whereas over sea there is accounted for slip
                !    vb_ represents the contribution in vd of the brownian diffusion [cm s-1]
                !    and vi_ represents the impaction [cm s-1]
                ! 

                vb_land   =(1./vkarman)*((ust_land/um)**2)*100.*um*(sc**(-2./3.))       ![cm s-1]
                vi_land   =(1./vkarman)*((ust_land/um)**2)*100.*um*(10.**(-3./st_land)) ![cm s-1]
                vkd_land  =(vb_land+vi_land)*1e-2                                       ![m s-1]
                vd_land=1./(raero(i,j)+1./vkd_land)
                !      if (vkd_land.gt.1.) write(*,999) &
                !        'after vb',i,j,'jtype',jtype,'jmod',jmod,'vb',vb_land,'vi',vi_land,&
                !           'um10',um,'ust_land',ust_land,'dc',dc,'relax',relax,'schmidt',sc,'stokes',st_land
             endif ! xland.gt.0
             if (xland.lt.0.99) then
                !--- Over sea:
                !    Brownian diffusion for rough elements, see Hummelshoj
                !    re is the reynolds stress:

                ust_sea=ustar_sea(i,j)
                st_sea  =max(0.1,(relax*(100.*ust_sea)**2.)/visc)
                re  =(100.*ust_sea*100.*sr_sea(i,j))/visc     ! [cm/s]*[cm]/[cm2/s]
                vbsea     =(1./3.)*100.*ust_sea*((sc**(-0.5))*re**(-0.5))
                visea     =100.*ust_sea*10.**(-3./st_sea)
                vkdaccsea =vbsea+visea
                vkd_sea   =((1.-alpha(i,j))*vkdaccsea+alpha(i,j)*bubble(i,j))*1e-2 ! [m s-1]
                vd_sea=1./(raero(i,j)+1./vkd_sea) ! [m s-1]
       !    if (vkd_sea.gt.1.)  write(*,999) 'sea',i,j,'jtype',jtype,'jmod',jmod,'vb',vbsea,&
       !    'vi',visea,'vkd_sea',vkd_sea,'alpha*bubble',alpha(i,j)*bubble(i,j),'alpha',alpha(i,j)

            endif ! xland.lt.0.99

            vd11a(i,j)= min(0.1,(1.-xland)*vd_sea + xland*vd_land)  ! [m s-1] limit to 10 cm/s

          enddo ! j=1,nlat180
          enddo ! i=1,nlon360 
!$OMP END PARALLEL


          ! gather before coarsening
          call gather( dgrid(iglbsfc), vd11a, global_field, 0, status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) then
             call coarsen_emission('vd', imr, jmr, global_field, vd_temp_global, avg_field, status)
             IF_NOTOK_RETURN(status=1)
          end if
          
          ! scatter and store
          do region = 1, nregions

             call scatter( dgrid(region), vd_temp(region)%surf, vd_temp_global(region)%surf, 0, status)
             IF_NOTOK_RETURN(status=1)

             vda(region,irdep)%surf = vd_temp(region)%surf
             
          end do

        enddo    ! irdep

      end if   ! nrdep > 0
      
      !
      ! ***
      !

      if ( okdebug ) write(*,*) 'dd_calc_rs: end'

    end subroutine dd_calc_rs

  end subroutine dd_surface_fields


  ! ***


#ifndef without_photolysis
  subroutine dd_coarsen_vd( vd11, ags, I1, J1, status )
#else
  subroutine dd_coarsen_vd( vd11, I1, J1, status )
#endif

    use dims      , only : iglbsfc, okdebug
    use chem_param, only : ndep, idep, vd_ncopy, vd_copy_itarget, vd_copy_isource
    use chem_param, only : ihno4, ih2o2, ino2
    use chem_param, only : names
    use toolbox   , only : coarsen_emission
    use toolbox   , only : escape_tm
#ifndef without_photolysis
    use photolysis_data, only : phot_dat
#endif
    use tm5_distgrid,    only : dgrid, Get_DistGrid, gather, scatter
    use ParTools,        only : isRoot

    
    ! --- in/out -------------------------------
    
    real, intent(inout)   ::  vd11(I1:,J1:,:)
#ifndef without_photolysis
    real, intent(inout)   ::  ags(I1:,J1:)
#endif
    integer, intent(in)   ::  i1, j1
    integer, intent(out)  ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = 'dd_coarsen_vd'

    integer,parameter  :: avg_field = 1

    ! --- local -------------------------------

    integer             ::  i,j, region
    integer             ::  jt,jdep
    integer             ::  idno2,idh2o2,idhno4
    !character(len=8)    ::  adum
    !character(len=2)    ::  bdum
    real, allocatable   ::  vdhelp(:,:)
    integer :: i01, i02, j01, j02, imr, jmr
    
    ! --- begin ------------------------------


    CALL GET_DISTGRID( DGRID(iglbsfc), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02, I_WRLD=imr, J_WRLD=jmr )

    if (isRoot) then
       allocate(vdhelp(imr,jmr))
    else
       allocate(vdhelp(1,1))
    end if
    !
    ! --
    ! Scale the surface resistance of a number of trace gases and aerosols
    ! NOx deposition for all NOx components separately and scaling 
    ! of NOx afterwards
    !
    
    ! *** vd(HNO4) = ( vd(NO2) + vd(H2O2) ) / 2.0
    
    ! search deposited tracer indices for NO2, H2O2, and HNO4
    idno2  = -1
    idh2o2 = -1 
    idhno4 = -1
    ! try to reset ...
    do jdep = 1,ndep   
       select case(idep(jdep))
       case(ihno4)
          idhno4 = jdep
       case(ih2o2)
          idh2o2 = jdep
       case(ino2)
          idno2 = jdep
       case default
       end select
    end do
    ! HNO4 present ? then compute its deposition velocity:
    if ( idhno4 > 0 ) then
      ! check ...
      if ( any( (/idno2,idh2o2/) < 0 ) ) then
        write (gol,'("deposition velocity for HNO4 computed from NO2 and H2HO2, but at least one not found:")'); call goErr
        write (gol,'("  idno2, idh2o2  : ",2i4)') idno2, idh2o2; call goErr
        TRACEBACK; status=1; return
      end if
      ! vd(HNO4) = ( vd(NO2) + vd(H2O2) ) / 2.0
      vd11(:,:,idhno4) = ( vd11(:,:,idno2) + vd11(:,:,idh2o2) ) / 2.0
    end if
    
    ! ***
    
    do jdep=1,ndep
       jt = idep(jdep)

       call gather( dgrid(iglbsfc), vd11(:,:,jdep), vdhelp, 0, status)
       IF_NOTOK_RETURN(status=1)

       if (isRoot) then
          call coarsen_emission('vd', imr, jmr, vdhelp, vd_temp_global, avg_field, status)
          IF_NOTOK_RETURN(status=1)
       end if
       
       do region = 1,nregions

          call scatter( dgrid(region), vd_temp(region)%surf, vd_temp_global(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          
          vd(region,jt)%surf = vd_temp(region)%surf

          !if ( okdebug ) then 
          !   adum=names(jt)
          !   write(bdum,'(i2.2)') region
          !   i=len_trim(adum)
          !   !cmk  call dumpfield(0,'vd.hdf',&
          !   !cmk     vd(region,jt)%surf(:,:),'vd_'//adum(1:i)//bdum)
          !end if
          
       end do


    end do !jdep

    ! copy fields:
    do jdep = 1, vd_ncopy
      do region = 1, nregions
        vd(region,vd_copy_itarget(jdep))%surf = vd(region,vd_copy_isource(jdep))%surf
      end do
    end do

#ifndef without_photolysis
    
    call gather( dgrid(iglbsfc), ags, vdhelp, 0, status)
    IF_NOTOK_RETURN(status=1)

    if (isRoot) then
       ! coarsen ags for photolysis...(temp in vd_temp)
       call coarsen_emission('ags',imr, jmr, vdhelp, vd_temp_global, avg_field, status)
       IF_NOTOK_RETURN(status=1)
    end if
    
    do region = 1,nregions

       call scatter( dgrid(region), vd_temp(region)%surf, vd_temp_global(region)%surf, 0, status)
       IF_NOTOK_RETURN(status=1)
       
       ! store coarsend albedo in photolysis data:
       phot_dat(region)%albedo = vd_temp(region)%surf

       ! update albedo statistics:
       !JEW: phot_dat(region)%ags_av = phot_dat(region)%ags_av + phot_dat(region)%albedo
       !JEW: phot_dat(region)%nalb_av = phot_dat(region)%nalb_av + 1
    end do
    
#endif

    ! Done
    deallocate(vdhelp)    
    status = 0
    
  end subroutine dd_coarsen_vd


END MODULE DRY_DEPOSITION