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

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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"
!
!------------------------------------------------------------------------------
!                    TM5                                                      !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE:     SEDIMENTATION
!
! !DESCRIPTION:  This module calculates sedimentation of aerosol particles. 
!                Also the deposition at the surface is calculated here since
!                 it uses similar routines. 
!
! It is assumed that the distribution is described by nmodes log-normal
! distributions as described in chem_param.F90. 
!
! Each mode has a number and mass and a sigma_lognormal. Number and mass are
! related and the mean aerosol radius can thus be calculated for each mode.
! 
! mass=pi*4./3.*radius**3.*number*exp(9./2.*ln2s) /1e18*density, with:
!
!     density = density of aerosol type
!     ln2s    = (alog(sigma_lognormal(mode)))**2
!     mass    = mode mass
!     number  = mode number
!
!    mode1 : accumulation
!    mode2 : coarse
!    mode3 : super coarse (ss) coarse
!
! For each mode a separate fall velocity is calculated according to the mass
! and number mean radii. Water take-up by seasalt particles is taken into
! account. This changes the density, radius, and sigma of the distribution.
!
! Also included is the deposition calculation. based on a lookup table
! calculated for a reference aerosol density (e.g. 1800 kg/m3) and a number of
! radii. This deposition curve is convoluted with the number/volume
! distribution of the aerosols.
!

! Again, for SS the water takeup is accounted for, and the effects on density,
! sigma and radius are calculated. The density has effect on the impaction
! term is the depotion calculation. This can be modeled by a shift in the
! radius. Thus the radii of the lookup table are adapted for density
! differences when impaction becomes important.
!
! routine Sedimentation_Init     : initialisation of sedimentation routine
! routine Sedimentation_Done     : free memory
! routine sedimentation_calcv    : calculate velocities
! routine deposition_calcv       : calculate velocities
! routine Sedimentation_Apply    : do the accutual sedimentation
!\\
!\\
! !INTERFACE: 
!
MODULE SEDIMENTATION
  !
  ! !USES:
  !
  use GO,            only: gol, goPr, goErr, goBug
  use GO,            only: GO_Timer_Def, GO_Timer_End, GO_Timer_Start
  use Binas,         only: grav, xmair, Avog
  use dims,          only: nregions
  use tm5_distgrid,  only: dgrid, Get_DistGrid, gather
  use partools,      only: isRoot
  use global_types,  only: d3_data, emis_data, aerosol_emis_data
  use chem_param,    only: ntracet, imsa, inh4, ino3_a
#ifdef with_m7
  use chem_param,    only: inus_n, iso4nus, isoanus, iais_n,  iso4ais, ibcais, ipomais, isoaais
  use chem_param,    only: iacs_n, iso4acs, ibcacs,  ipomacs, issacs,  iduacs, isoaacs
  use chem_param,    only: icos_n, iso4cos, ibccos,  ipomcos, isscos,  iducos, isoacos
  use chem_param,    only: iaii_n, ibcaii,  ipomaii, isoaaii, iaci_n,  iduaci, icoi_n,  iducoi 
#else
  use chem_param,    only: iso4
#endif
  use chem_param,    only: nmod
#ifdef with_budgets
  use budget_global, only: nbudg, nbud_vg, nzon_vg
#endif

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !  
  public Sedimentation_Init, Sedimentation_Done
  public Sedimentation_Apply
  public calculate_rh
  public aerosol_radius
  !
  ! !PUBLIC DATA MEMBERS:
  !  
  real,          dimension(nregions), public :: sum_drydep, sum_sedimentation   ! budgets for tracer #1 (significant on root only)
  type(d3_data), dimension(nregions), public :: rh, rha               ! rh (0-1), cloudfree vs. allsky

#ifdef with_m7
  integer, parameter, public                           :: nsed=33  ! number of tracers for which sedimentation is applied
  integer, parameter, dimension(nsed), public          :: ised = (/   inh4, ino3_a, imsa,     &
                                                              inus_n, iso4nus, isoanus, iais_n,  iso4ais, ibcais, ipomais, isoaais, &
                                                              iacs_n, iso4acs, ibcacs,  ipomacs, issacs,  iduacs, isoaacs,  &
                                                              icos_n, iso4cos, ibccos,  ipomcos, isscos,  iducos, isoacos,  &
                                                              iaii_n, ibcaii,  ipomaii, isoaaii, iaci_n,  iduaci, icoi_n,  &
                                                              iducoi /)
#else
  !  sedimentation of few gas-phase species
  integer, parameter, public                           :: nsed= 4  ! number of tracers for which sedimentation is applied
  integer, parameter, dimension(nsed), public          :: ised = (/ iso4, inh4, ino3_a, imsa /)
#endif

  ! TB This was made public, to allow getting indices in AerChemMIP output
  integer, dimension(ntracet),public                  :: sindex  ! pointer to match ntracet array with nsed
  !
  ! !PRIVATE DATA MEMBERS:
  !
  ! Mean values at the surface, per mode (i,j,mode)
  type(aerosol_emis_data), dimension(nregions) :: vn_deposition_mean    ! number mean deposition velocity (m/s)
  type(aerosol_emis_data), dimension(nregions) :: vn_sedimentation_mean ! number mean sedim. velocity at surface (m/s)
  type(aerosol_emis_data), dimension(nregions) :: vm_deposition_mean    ! mass mean deposition velocity (m/s)
  type(aerosol_emis_data), dimension(nregions) :: vm_sedimentation_mean ! mass mean sedim. velocity at surface (m/s)

  type(d3_data),   dimension(nregions,nmod) :: vn_sedimentation     ! number sedimentation velocities (Pa/s)
  type(emis_data), dimension(nregions,nmod) :: vn_deposition        ! number deposition velocity (Pa/s)
  type(d3_data),   dimension(nregions,nmod) :: vm_sedimentation     ! mass sedimentation velocities (Pa/s)
  type(emis_data), dimension(nregions,nmod) :: vm_deposition        ! mass deposition velocity (Pa/s)
  integer,         dimension(nregions)      :: n_deposition_mean    ! counter
  integer,         dimension(nregions)      :: n_sedimentation_mean ! counter
  logical, parameter                        :: vd_in_deposition = .false.  ! is sedimentation treated in deposition routine?
  integer, parameter                        :: ndp = 3 ! limit of the velocity to ndp times the layer thickness
                                                       ! iteration will account for a fall length through multiple layers
  real, parameter                           :: m_to_pa =  7.24e16*grav*xmair*1e3/Avog  !factor from m/s --> Pa/s
  !
  ! !PUBLIC TYPES:
  !
#ifdef with_budgets
  ! budget terms on grid basis:  
  type, public :: buddep_data
     real, dimension(:,:,:,:), pointer  :: sed    ! im,jm,nbud_vg,nsed
  end type buddep_data
  type(buddep_data),dimension(nregions),target,public  :: buddep_m7_dat
#endif

  integer :: itim_appl
  
  real                        :: sigma,      density
  character(len=*), parameter :: mname = 'sedimentation'
  !
  ! !REVISION HISTORY: 
  !    Feb 2004 by MK - 
  !    2 Apr 2010 - P. Le Sager - Optimize arrays de/allocation if m7 is not used.
  !   19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !  (1) Only relative humidity is used when not using m7.
  !  (2) sedimentation is defined at the bottom of the gridbox (see init)
  !  (3) Note the switch from emiss_data to aerosol_emiss_data type for the 
  !       *_mean variables (simplify and speedup mpi comm in sedimentation_done)
  !
  ! !TODO:
  !
  !EOP
  !------------------------------------------------------------------------------

CONTAINS

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SEDIMENTATION_INIT
  !
  ! !DESCRIPTION: allocate memory, map ised
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SEDIMENTATION_INIT ( status )
    !
    ! !USES:
    !
    use dims,       only : lm, iglbsfc
    use chem_param, only : ntracet
    use meteodata,  only : Set, spm_dat, temper_dat, humid_dat, cc_dat, phlb_dat
    use meteodata,  only : lsp_dat, cp_dat
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !    2 Apr 2010 - P. Le Sager - Added test on "with_m7"
    !   20 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Sedimentation_Init'

    integer               :: region, lmr, i
    integer               :: imode, n, i1, i2, j1, j2

    ! --- begin ---------------------------------

    ! meteo to be used
    do region = 1, nregions
       call Set( temper_dat(region), status, used=.true. )
       call Set(   phlb_dat(region), status, used=.true. )
       call Set(  humid_dat(region), status, used=.true. )
       call Set(     cc_dat(region), status, used=.true. )
       call Set(    lsp_dat(region), status, used=.true. )
       call Set(     cp_dat(region), status, used=.true. )
    end do

#ifdef with_m7

    ! Calculate mapping of ntracer to sedimentation arrays:
    sindex(:) = 0

    do i=1,nsed
       do n=1,ntracet
          if(ised(i) == n) sindex(n) = i
       enddo
    enddo

#endif

    ! Allocate
    do region = 1, nregions

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
       lmr=lm(region)
       
#ifdef with_m7      

#ifdef with_budgets
       sum_sedimentation(region) = 0.0
       sum_drydep(region)        = 0.0
       allocate( buddep_m7_dat(region)%sed(i1:i2, j1:j2, nbud_vg, nsed) )
       buddep_m7_dat(region)%sed = 0.0
#endif

       allocate(vn_deposition_mean   (region)%surf(i1:i2, j1:j2, nmod))      ! defined at bottom of surface layer
       allocate(vn_sedimentation_mean(region)%surf(i1:i2, j1:j2, nmod))      ! defined at bottom of surface layer
       allocate(vm_sedimentation_mean(region)%surf(i1:i2, j1:j2, nmod))      ! defined at bottom of surface layer
       allocate(vm_deposition_mean   (region)%surf(i1:i2, j1:j2, nmod))      ! defined at bottom of surface layer
       vn_deposition_mean   (region)%surf = 0.0
       vn_sedimentation_mean(region)%surf = 0.0
       vm_sedimentation_mean(region)%surf = 0.0
       vm_deposition_mean   (region)%surf = 0.0
       
       do imode = 1, nmod
          allocate(vn_sedimentation(region,imode)%d3  (i1:i2, j1:j2, lmr))   ! defined at bottom of layer
          allocate(vn_deposition   (region,imode)%surf(i1:i2, j1:j2     ))   ! defined at bottom of surface layer
          allocate(vm_sedimentation(region,imode)%d3  (i1:i2, j1:j2, lmr))   ! defined at bottom of layer
          allocate(vm_deposition   (region,imode)%surf(i1:i2, j1:j2     ))   ! defined at bottom of surface layer
          vn_sedimentation(region,imode)%d3   = 0.0
          vn_deposition   (region,imode)%surf = 0.0
          vm_sedimentation(region,imode)%d3   = 0.0
          vm_deposition   (region,imode)%surf = 0.0
       enddo

#endif

       allocate( rh(region)%d3(i1:i2, j1:j2, lmr)) 
       allocate(rha(region)%d3(i1:i2, j1:j2, lmr)) 
       n_deposition_mean(region)    = 0
       n_sedimentation_mean(region) = 0
    enddo

    call GO_Timer_Def( itim_appl, 'sedimentation appl', status )
    IF_NOTOK_RETURN(status=1)
    
    status = 0

  END SUBROUTINE SEDIMENTATION_INIT
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SEDIMENTATION_DONE
  !
  ! !DESCRIPTION: gather/aggregate/write budgets, and free memory
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SEDIMENTATION_DONE( status )
    !
    ! !USES:
    !
    use io_hdf,        only : DFACC_WRITE, io_write, DFNT_INT32
    use dims,          only : im, jm

#ifdef with_budgets
    use budget_global, only : budget_file_global, nbud_vg, budg_dat, nbudg, NHAB
    use file_hdf,      only : THdfFile, TSds
    use file_hdf,      only : Init, Done, WriteAttribute, WriteData, SetDim
    use Dims,          only : region_name
    use partools,      only : par_reduce_element
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)     ::  status
    !
    ! !REMARKS:
    !  
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Sedimentation_Done'

    Integer                                  :: region, i1, i2, j1, j2
#ifdef with_m7
#ifdef with_budgets
    type(THdfFile)                           :: io
    type(TSds)                               :: sds
    integer                                  :: j, i, n, nzone, nzone_v
    real, dimension(nbudg, nbud_vg, nsed)    :: budsedg
    real, dimension(nbudg, nbud_vg, nsed)    :: budsedg_all ! for buggy MPI
    real, dimension(:,:),    allocatable     :: budget2d    ! to gather budget terms on PEs
    real, dimension(:,:,:),  allocatable     :: budget3d
    real, dimension(:,:,:,:),allocatable     :: budget4d
#endif
#endif
    integer                                  :: imode

    ! --- begin ------------------------------

#ifdef with_m7
#ifdef with_budgets

    ! open file
    if( isRoot ) then
       call Init(io, budget_file_global, 'write', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io,'sum_drydep_m7',sum_drydep,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io,'sum_sedimentation_m7',sum_sedimentation,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io,'nsed',nsed,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io,'ised',ised,status)
       IF_NOTOK_RETURN(status=1)
    endif

    budsedg(:,:,:) = 0.0

    REG: do region=1,nregions   

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
              
       !-- bin aggregated sed budget
       do n=1,nsed
       do j=j1,j2
       do i=i1,i2
                nzone = budg_dat(region)%nzong(i,j)
                budsedg(nzone,:,n) = budsedg(nzone,:,n) + buddep_m7_dat(region)%sed(i,j,:,n)
       end do
       end do
       end do
       
       !-- write Non-Horizontally-Aggregated-Budgets
       if (NHAB) then

          ! gather sed budget and write to file
          if (isRoot )then
             allocate(budget4d(im(region), jm(region), nbud_vg, nsed))
          else
             allocate(budget4d(1,1,1,1))
          end if
              
          call GATHER( dgrid(region), buddep_m7_dat(region)%sed, budget4d, 0, status)
          IF_NOTOK_RETURN(status=1)

          if(isRoot )then
             call Init(Sds,io, 'buddep_dat_sed_'//region_name(region),(/im(region),jm(region),nbud_vg,nsed/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute(Sds,'ised',ised,status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             call SetDim(Sds, 2, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status)
             call SetDim(Sds, 3, 'nsed','sedimentation', (/(j, j=1,nsed)/), status)
             call WriteData(Sds, budget4d ,status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif

          deallocate(budget4d)
          
          
          ! gather deposition velocities and write to file
          if (isRoot ) then
             allocate(budget3d(im(region), jm(region), nmod))
          else
             allocate(budget3d(1,1,1))
          end if

          call GATHER( dgrid(region), vn_deposition_mean(region)%surf, budget3d, 0, status)
          IF_NOTOK_RETURN(status=1)

          if(isRoot ) then
             if (n_deposition_mean(region)==0)then
                budget3d = 0.
             else
                budget3d = budget3d / n_deposition_mean(region)
             end if
          
             call Init(Sds,io, 'vd_n_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget3d ,status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif


          ! gather and write deposition velocities (mass tracers)
          call GATHER( dgrid(region), vm_deposition_mean(region)%surf, budget3d, 0, status)
          IF_NOTOK_RETURN(status=1)
          
          if(isRoot ) then
             if (n_deposition_mean(region)==0)then
                budget3d = 0.
             else
                budget3d = budget3d / n_deposition_mean(region)
             end if
             
             call Init(Sds,io, 'vd_m_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget3d ,status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif


          ! gather and write sedimentation velocities
          call GATHER( dgrid(region), vn_sedimentation_mean(region)%surf, budget3d, 0, status)
          IF_NOTOK_RETURN(status=1)

          if(isRoot ) then

             if (n_sedimentation_mean(region)==0)then
                budget3d = 0.
             else
                budget3d = budget3d / n_sedimentation_mean(region)
             end if

             call Init(Sds,io, 'vs_n_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget3d ,status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif

          ! gather and write sedimentation velocities (mass tracers)
          call GATHER( dgrid(region), vm_sedimentation_mean(region)%surf, budget3d, 0, status)
          IF_NOTOK_RETURN(status=1)

          if(isRoot) then
             if (n_sedimentation_mean(region)==0)then
                budget3d = 0.
             else
                budget3d = budget3d / n_sedimentation_mean(region)
             end if

             call Init(Sds,io, 'vs_m_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
             IF_NOTOK_RETURN(status=1)
             call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, budget3d ,status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif

          deallocate(budget3d)
          
       endif ! NHAB
    enddo REG  ! regions

    !------- Gather and write aggregated budget

    ! Sum up contribution from each proc into root array
    call PAR_REDUCE_ELEMENT( budsedg, 'SUM', budsedg_all, status )
    IF_NOTOK_RETURN(status=1)

    if(isRoot )then
       call Init(Sds, io, 'budsed_m7', (/nbudg,nbud_vg,nsed/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call SetDim(Sds, 0, 'nbudg', 'horizontal region', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg', 'vertical layer', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'nsed', 'sedimentation', (/(j, j=1,nsed)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, budsedg_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
       call Done(io, status)
       IF_NOTOK_RETURN(status=1)
    endif

#endif /* BUDGET */
#endif /* M7 */


    !-------- Done
    do region = 1, nregions

#ifdef with_m7
#ifdef with_budgets
       deallocate( buddep_m7_dat(region)%sed)
#endif
       deallocate( vn_deposition_mean   (region)%surf )
       deallocate( vn_sedimentation_mean(region)%surf )
       deallocate( vm_deposition_mean   (region)%surf )
       deallocate( vm_sedimentation_mean(region)%surf )
              
       do imode = 1, nmod
          deallocate(vn_sedimentation(region,imode)%d3  )
          deallocate(vn_deposition   (region,imode)%surf)
          deallocate(vm_sedimentation(region,imode)%d3  )
          deallocate(vm_deposition   (region,imode)%surf)
       enddo
#endif

       deallocate(rh(region)%d3)
       deallocate(rha(region)%d3)
       
    enddo

    status = 0

  END SUBROUTINE SEDIMENTATION_DONE
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   DEPOSITION_CALCV
  !
  ! !DESCRIPTION: calculate velocities
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DEPOSITION_CALCV( region, status )
    !
    ! !USES:
    !
    use global_data,    only : mass_dat,         region_dat
    use meteodata,      only : spm_dat,          temper_dat, cc_dat
    use chem_param,     only : sigma_lognormal,  density_ref
    use chem_param,     only : mode_start,       mode_end,    names
    use dims,           only : at, bt, nsrce, tref, im, jm, lm, okdebug
    use chem_param,     only : nrdep, lur
    use dry_deposition, only : vd=>vda
    use binas,          only : pi
    use partools,       only : par_reduce
#ifdef with_m7
    use m7_data,    only : dens_mode,     rw_mode       ! density (kg/m3) r (m)
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)               :: region
    !                                   
    ! !OUTPUT PARAMETERS:               
    !                                   
    integer, intent(out)              :: status
    !
    ! !REVISION HISTORY: 
    !    2 Apr 2010 - P. Le Sager - 
    !   21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !  (1) Called only if m7 used.
    !   
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/DEPOSITION_CALCV'

    real, dimension(:,:,:),   pointer   :: p      ! surface pressure (Pa)                                  
    real, dimension(:,:,:),   pointer   :: t      ! temperature (K)
    real, dimension(:,:,:),   pointer   :: cc     ! cloud_cover (0-1)
    real, dimension(:,:,:,:), pointer   :: rm     ! tracer array (kg)
    real                                :: pb     ! pressure at bottom of box (Pa)
    real                                :: dt     ! general timestep (s)
    real                                :: dp     ! pressure difference layer
    integer                             :: i,j,l, mode, n, proc, itn, itracer
    real                                :: temp   ! temperature
    real                                :: to_pascal           
    real,    dimension(nmod)            :: lns              
    real                                :: radius           

    real,    dimension(nmod,2)          :: vd_mean_all, vd_max_all
    real,    dimension(nmod,2)          :: vd_mean,     vd_max           
    integer, dimension(nmod,2)          :: nvd,         nvd_all

    real,    dimension(nmod)            :: r_mean_all,  r_max_all
    real,    dimension(nmod)            :: r_mean,      r_max
    integer, dimension(nmod)            :: nr,          nr_all

    real,    dimension(nrdep)           :: d_aer         ! diameter vd loopup table (um)
    real,    dimension(nrdep)           :: nnumb,nvolume ! number and volume distribution
    real,    dimension(nrdep)           :: vdi           ! for the integration
    real,    dimension(nrdep)           :: vdi_def       ! for the integration
    real,    dimension(nrdep)           :: lure          ! effective loo
    real,    dimension(nrdep)           :: ddpi          ! integration bin-sizes
    real,    dimension(nrdep+1)         :: dlogdp,ddp    ! integration edges
    real,    dimension(nrdep)           :: logdp,logde   ! log(diameter)
    real                                :: dpg, ntot, vt       
    real                                :: lnsigma, rho_aer, sigma
    integer                             :: ir, ir1, nshift, nrd, nglob

    integer :: i1, i2, j1, j2
    !________________________start_________________________________ 

    nullify(rm)
    nullify(p)
    nullify(t)
    nullify(cc)
    rm => mass_dat  (region)%rm
    p  => spm_dat   (region)%data
    t  => temper_dat(region)%data
    cc => cc_dat    (region)%data

    dt = float(nsrce)/(2*tref(region))   ! timestep

    do mode =1,nmod
       lns(mode)      = log(sigma_lognormal(mode))
    enddo

    ! calculate the binsizes (um) around the radii of the pre-calculated vd's
    d_aer    = 2.0*lur      ! diameter (um)
    logdp    = log10(d_aer) ! log(diameter)

    n_deposition_mean(region) = n_deposition_mean(region) + 1

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )


    if (okdebug) then
       if (isRoot) then
          write(gol,*)  '________________________________________________________________________________' ; call goPr
          write(gol,*)  'Statistics for the deposition routine:' ; call goPr
          write(gol,*)  '________________________________________________________________________________' ; call goPr
          write(gol,*)  'level mode r_mean  r_max(um) vd_n(cm/s) vd_max(cm/s)  vd_m(cm/s) vd_mmax(cm/s)  ' ; call goPr
          write(gol,*)  '________________________________________________________________________________' ; call goPr
       end if
       vd_mean = 0.0
       r_mean  = 0.0
       vd_max  = 0.0
       r_max   = 0.0
       nvd     = 0
       nr      = 0
    endif

    !$OMP PARALLEL &
    !$OMP  default (none) &
    !$OMP  shared  (region, t, at, bt, p, dt, vd, &
    !$OMP           rm, rh, logdp, &
    !$OMP           vn_deposition_mean, vn_deposition, vm_deposition_mean, vm_deposition, &
    !$OMP           lur, sigma_lognormal, mode_start, mode_end) &
    !$OMP  shared  (okdebug, vd_mean, vd_max, nvd, r_max, r_mean, nr) &
#ifdef with_m7
    !$OMP shared   (rw_mode, dens_mode) &
#endif

    !$OMP  private (i, j, l, temp, pb, dp, to_pascal, vdi_def, &
    !$OMP           itn, radius, vdi, sigma, lnsigma, density, &
    !$OMP           lure, dlogdp, logde, nshift, nrd, &
    !$OMP           ddp, ddpi, d_aer, dpg, ntot, nnumb, nvolume, vt) 

    l = 1

    do j=j1,j2
       do i=i1,i2
          temp      = t(i,j,1)                                   ! at surface to temp box
          pb        = at(l) + bt(l)*p(i,j,1)                     ! pressure at bottom of box (Pa)
          dp        = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) ! layer thickness
          to_pascal = m_to_pa*dt*pb/temp                         ! convert from m/s ---> Pa/timestep

          do n=1,nrdep
             vdi_def(n) = vd(region,n)%surf(i,j)
          enddo

#ifdef with_m7
          M7MODES: do mode = 1, nmod

             itn = mode_start(mode)   ! position of number tracer

             ! compute radius:
             radius = rw_mode(region,mode)%d3(i,j,1)

             ! initial deposition velocities for increasing radia:
             vdi = vdi_def

             sigma   = sigma_lognormal(mode)
             lnsigma = log(sigma)
             density = dens_mode(region,mode)%d3(i,j,1)

             if(okdebug) then 
                if(radius > tiny(radius)) then
                   r_mean(mode) = r_mean(mode) + radius
                   nr(mode)     = nr(mode) + 1
                   r_max(mode)  = max(r_max(mode), radius)
                endif
             endif

             RADENS: if (radius > 1e-11 .and. density > 1e-2) then

                ! account for density different than density_ref of the look-up table (lur --> vdi):
                lure(:) = lur(:)
                logde(:) = logdp(:)
                do ir = 2, nrdep
                   if(vdi(ir) > vdi(ir-1)) exit  ! for bigger r's :  impaction dominates (density effects)
                   if ( ir == nrdep ) exit   ! trap upper boundary
                enddo
                do ir1 = ir, nrdep
                   lure(ir1) = lur(ir1)*sqrt(density_ref/density)
                   logde(ir1) = log10(2*lure(ir1))
                enddo

                ! compress look-up table such that radii are increasing monotonic:
                nshift = 0
                ir1 = ir
                do
                   if ( logde(ir1) > logde(ir-1) ) exit
                   nshift = nshift + 1
                   ir = ir -1
                   if(ir == 1) exit
                enddo
                nrd = nrdep - nshift
                if (nshift > 0) then 
                   do ir1 = ir, nrd
                      logde(ir1) = logde(ir1+nshift)
                      lure(ir1) = lure(ir1+nshift)
                      vdi(ir1) = vdi(ir1+nshift)
                   enddo
                endif

                ! do the integration of the shifted lookup table:
                dlogdp(1) = -3.0 
                do ir=2,nrd 
                   dlogdp(ir) = 0.5*(logde(ir-1)+logde(ir))   ! take middle of the log scale
                enddo
                dlogdp(nrd+1) = 3.0   ! 1000 um
                ddp(1:nrd+1)  = 10**dlogdp(1:nrd+1)
                ddpi(1:nrd) = ddp(2:nrd+1)-ddp(1:nrd)   ! integration intervals (um)
                d_aer(1:nrd) = 2.0*lure(1:nrd)

                ! perform convolution with log-normal distribution:
                dpg = 2*radius*1e6   ! diameter in um
                ntot = rm(i,j,1,itn)

                ! calculate the distribution (number and mass) over the deposition bins:
                if(ntot > 1.0 .and. radius > tiny(radius) ) then ! you need some aerosol!
                   do n=1,nrd
                      nnumb(n) = ntot/(sqrt(2.*pi)*d_aer(n)*lnsigma)*exp(-(log(d_aer(n))-log(dpg))**2/(2*lnsigma**2))
                      nvolume(n) = nnumb(n)*(pi/6.0)*d_aer(n)**3
                   enddo
                   vt = sum(nnumb(1:nrd)*ddpi(1:nrd)*vdi(1:nrd))/sum(nnumb(1:nrd)*ddpi(1:nrd))
                else   
                   vt = 0.0
                endif

                vn_deposition_mean(region)%surf(i,j,mode) = vn_deposition_mean(region)%surf(i,j, mode) + vt
                vn_deposition(region,mode)%surf(i,j) = min(to_pascal*vt,ndp*dp)   ! in Pa/timestep downwards

                if(okdebug) then 
                   if ( vt > tiny(vt) ) then
                      vd_mean(mode,1) = vd_mean(mode,1) + vt
                      vd_max(mode,1)  = max(vd_max(mode,1) , vt)
                      nvd(mode,1)     = nvd(mode,1) + 1
                   endif
                endif

                ! for mass:
                if(ntot > 1.0 .and. radius > tiny(radius) ) then ! you need some aerosol!
                   vt = sum(nvolume(1:nrd)*ddpi(1:nrd)*vdi(1:nrd))/sum(nvolume(1:nrd)*ddpi(1:nrd))
                else
                   vt = 0.0
                endif
                vm_deposition_mean(region)%surf(i,j, mode) = vm_deposition_mean(region)%surf(i,j, mode) + vt
                vm_deposition(region,mode)%surf(i,j) = min(to_pascal*vt,ndp*dp)   ! in Pa/timestep downwards

                if(okdebug) then 
                   if ( vt > tiny(vt) ) then
                      vd_mean(mode,2) = vd_mean(mode,2) + vt
                      vd_max(mode,2)  = max(vd_max(mode,2) , vt)
                      nvd(mode,2)     = nvd(mode,2) + 1
                   endif
                endif   ! 

             else
                vm_deposition(region,mode)%surf(i,j) = 0.0
                vn_deposition(region,mode)%surf(i,j) = 0.0
             endif RADENS

          end do M7MODES
          
#endif /* M7 */
       enddo  !i
    enddo     !j
    !$OMP END PARALLEL

    if(okdebug) then
       do mode = 1, nmod

          call Par_reduce( r_mean(mode)   , 'sum', r_mean_all(mode)   , status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( nr(mode)       , 'sum', nr_all(mode)       , status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( r_max(mode)    , 'max', r_max_all(mode)    , status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( vd_mean(mode,1), 'sum', vd_mean_all(mode,1), status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( nvd(mode,1)    , 'sum', nvd_all(mode,1)    , status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( vd_max(mode,1) , 'max', vd_max_all(mode,1) , status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( vd_mean(mode,2), 'sum', vd_mean_all(mode,2), status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( nvd(mode,2)    , 'sum', nvd_all(mode,2)    , status)
          IF_NOTOK_RETURN(status=1)
          call Par_reduce( vd_max(mode,2) , 'max', vd_max_all(mode,2) , status)
          IF_NOTOK_RETURN(status=1)

          if (isRoot) then
             if(nr_all(mode) > 0) then
                r_mean_all(mode) = r_mean_all(mode)*1e6/nr_all(mode)
             else
                r_mean_all(mode) = 0.
             end if

             if(nvd_all(mode,1) > 0) then
                vd_mean_all(mode,1) = 1e2*vd_mean_all(mode,1)/nvd_all(mode,1)
             else
                vd_mean_all(mode,1) = 0.
             end if

             if(nvd_all(mode,2) > 0) then
                vd_mean_all(mode,2) = 1e2*vd_mean_all(mode,2)/nvd_all(mode,2)
             else
                vd_mean_all(mode,2) = 0.
             end if

             print '(i6,i4,2f10.4,2f12.6,2f12.6)', 0, mode, r_mean_all(mode), 1e6*r_max_all(mode), &
                  vd_mean_all(mode,1), 1e2*vd_max_all(mode,1), vd_mean_all(mode,2), 1e2*vd_max_all(mode,2)
          end if
       enddo
       write(gol,*)  '________________________________________________________________________________' ; call goPr
    endif

    status=0

  END SUBROUTINE DEPOSITION_CALCV
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SEDIMENTATION_CALCV
  !
  ! !DESCRIPTION: calculate velocities
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SEDIMENTATION_CALCV(region, status)
    !
    ! !USES:
    !
    use global_data, only : mass_dat, region_dat
    use meteodata,   only : spm_dat,  temper_dat, cc_dat
    use chem_param,  only : sigma_lognormal, names, mode_start, mode_end
    use dims,        only : at, bt, nsrce, tref, im, jm, lm, okdebug
#ifdef with_m7
    use m7_data,     only : rw_mode, dens_mode
#endif
    use partools,    only : par_reduce
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)     :: region
    !                                   
    ! !OUTPUT PARAMETERS:               
    !                                   
    integer, intent(out)    :: status
    !
    ! !REVISION HISTORY: 
    !    2 Apr 2010 - P. Le Sager - 
    !   21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   (1) Called only if m7 used.
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Sedimentation_calcv'
    
    real, dimension(:,:,:),   pointer  :: p      ! surface pressure (Pa)
    real, dimension(:,:,:),   pointer  :: t      ! temperature (K)
    real, dimension(:,:,:),   pointer  :: cc     ! cloud cover (0-1)
    real, dimension(:,:,:),   pointer  :: q      ! humidity (kg/kg)
    real, dimension(:,:,:,:), pointer  :: rm     ! tracer array (kg)
    real                               :: pb     ! pressure at bottom of box (Pa)
    real                               :: dt     ! general timestep (s)
    real                               :: dp     ! pressure difference layer
    integer                            :: i,j,l, mode, proc, itn, n, nglob
    real                               :: temp   ! temperature
    real                               :: to_pascal
    real                               :: slinnfac 
    real                               :: zvis, rho_air, radius, Dpg, vt

    real,    dimension(nmod,2)         :: vd_mean_all, vd_max_all
    real,    dimension(nmod,2)         :: vd_mean,     vd_max           
    integer, dimension(nmod,2)         :: nvd,         nvd_all

    real,    dimension(nmod)           :: r_mean_all,  r_max_all
    real,    dimension(nmod)           :: r_mean,      r_max
    integer, dimension(nmod)           :: nr,          nr_all

    real                               :: lnsigma, sigma
    integer                            :: i1, j1, i2, j2

    !________________________start_________________________________

    nullify(rm)
    nullify(p)
    nullify(cc)
    nullify(t)
    p  => spm_dat   (region)%data
    t  => temper_dat(region)%data
    cc => cc_dat    (region)%data
    rm => mass_dat  (region)%rm
    
    dt = float(nsrce)/(2*tref(region))   ! timestep

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    ! transfer sedimentation vel. from m/s to Pa/s (note define this as positive = falling)
    ! dp = -rho*g*dz 
    ! rho is calculated from temperature and pressure....
    ! 7.24e16*p/T (#/cm3) (pV=nRT--> nR/V = p/T )
    ! #/cm3 ---> kg/m3  f = xmair*1e3/Avog
    ! Thus: rho (kg/m3) = 7.24e16*p/T  * xmair*1e3/Avog
    ! and (dPa) = -rho*g*(vdep)*dt = -7.24e16*grav*xmair*1e3/Avog * (p/T) *dt

    n_sedimentation_mean(region) = n_sedimentation_mean(region) + 1

    if(okdebug.and.isRoot) then
       write(gol,*)  '________________________________________________________________________________' ; call goPr
       write(gol,*)  'Statistics for the sedimentation routine:' ; call goPr
       write(gol,*)  '________________________________________________________________________________' ; call goPr
       write(gol,*)  'level mode r_mean  r_max(um) vd_n(cm/s) vd_max(cm/s)  vd_m(cm/s) vd_mmax(cm/s)  ' ; call goPr
       write(gol,*)  '________________________________________________________________________________' ; call goPr
    endif

#ifdef with_m7
    !$OMP PARALLEL &
    !$OMP  default (none) &
    !$OMP  shared  (region, t, at, bt, p, dt,  &
    !$OMP           rm, lm, & 
    !$OMP           vn_sedimentation_mean, vn_sedimentation, vm_sedimentation_mean, vm_sedimentation, &
    !$OMP           sigma_lognormal, mode_start, mode_end) &
    !$OMP  shared  (okdebug, vd_mean, vd_max, nvd, r_max, r_mean, nr) &
    !$OMP shared   (rw_mode, dens_mode) &
    !$OMP  private (i, j, l, temp, pb, dp, zvis, rho_air, &
    !$OMP           to_pascal, mode, itn, radius, sigma, lnsigma, &
    !$OMP           density, slinnfac, dpg, vt)

    LEVS: do l=1, lm(region)

       if(okdebug) then
          vd_mean = 0.0
          r_mean = 0.0
          vd_max = 0.0
          r_max = 0.0
          nvd  = 0
          nr  = 0
       endif

       do j=j1, j2
       do i=i1, i2

             if(l == 1) then 
                temp = t(i,j,1)                    ! at surface to temp box
             else
                temp = 0.5*(t(i,j,l)+t(i,j,l-1))   ! interpolate to bottom of box
             endif

             pb        = at(l) + bt(l)*p(i,j,1)                     ! pressure at bottom of box (Pa)
             dp        = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) ! layer thickness
             zvis      = 0.0000014963*temp*sqrt(temp)/(temp+120.)   ! viscocity  (kg/ms)
             rho_air   = 7.24e16*pb/temp  * xmair*1e3/Avog          ! kg/m3
             to_pascal = m_to_pa*dt*pb/temp                         ! convert from m/s ---> Pa/timestep

             M7MODES: do mode = 1, nmod

                radius = rw_mode(region,mode)%d3(i,j,l)

                sigma = sigma_lognormal(mode)
                lnsigma = log(sigma)
                density = dens_mode(region,mode)%d3(i,j,l)

                slinnfac = exp(2.0*lnsigma*lnsigma)

                if(okdebug) then
                   if(radius > tiny(radius)) then
                      r_mean(mode) = r_mean(mode) + radius
                      nr(mode) = nr(mode) + 1
                      r_max(mode) = max(r_max(mode), radius)
                   endif
                endif

                ! for number: take first mode (Seinfeld, 1986, pg 286)
                dpg = radius*2.0*exp(lnsigma*lnsigma)   ! diameter in m  
                vt =  fall1(pb,Dpg,zvis,temp,rho_air,density)

                if(okdebug) then
                   if ( vt > tiny(vt) ) then
                      vd_mean(mode,1) = vd_mean(mode,1) + vt
                      vd_max(mode,1) = max(vd_max(mode,1) , vt)
                      nvd(mode,1)     = nvd(mode,1) + 1
                   endif
                endif

                vn_sedimentation(region,mode)%d3(i,j,l) = min(to_pascal*vt,ndp*dp)   ! in Pa/timestep downwards

                if(l == 1) then
                   vn_sedimentation_mean(region)%surf(i,j,mode) = &
                        vn_sedimentation_mean(region)%surf(i,j,mode) + vt
                endif

                ! for mass: the third moment
                dpg = radius*2.0*exp(3*lnsigma*lnsigma)  ! diameter in m
                vt =  fall1(pb,Dpg,zvis,temp,rho_air,density)

                if(okdebug) then
                   if ( vt > tiny(vt) ) then
                      vd_mean(mode,2) = vd_mean(mode,2) + vt*slinnfac
                      vd_max(mode,2) = max(vd_max(mode,2) , vt*slinnfac)
                      nvd(mode,2)     = nvd(mode,2) + 1
                   endif
                endif

                vm_sedimentation(region,mode)%d3(i,j,l) = min(to_pascal*vt*slinnfac,ndp*dp)  ! multiply with slinnfac

                if(l == 1) then
                   vm_sedimentation_mean(region)%surf(i,j,mode) = &
                        vm_sedimentation_mean(region)%surf(i,j,mode) + vt*slinnfac
                endif

             enddo M7MODES
          enddo ! i
       enddo ! j

       if(okdebug) then
          do mode = 1, nmod

             call Par_reduce( r_mean(mode)   , 'sum', r_mean_all(mode)   , status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( nr(mode)       , 'sum', nr_all(mode)       , status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( r_max(mode)    , 'max', r_max_all(mode)    , status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( vd_mean(mode,1), 'sum', vd_mean_all(mode,1), status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( nvd(mode,1)    , 'sum', nvd_all(mode,1)    , status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( vd_max(mode,1) , 'max', vd_max_all(mode,1) , status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( vd_mean(mode,2), 'sum', vd_mean_all(mode,2), status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( nvd(mode,2)    , 'sum', nvd_all(mode,2)    , status)
             IF_NOTOK_RETURN(status=1)
             call Par_reduce( vd_max(mode,2) , 'max', vd_max_all(mode,2) , status)
             IF_NOTOK_RETURN(status=1)

             if (isRoot) then
                if(nr_all(mode) > 0) then
                   r_mean_all(mode) = r_mean_all(mode)*1e6/nr_all(mode)
                else
                   r_mean_all(mode) = 0.
                end if

                if(nvd_all(mode,1) > 0) then
                   vd_mean_all(mode,1) = 1e2*vd_mean_all(mode,1)/nvd_all(mode,1)
                else
                   vd_mean_all(mode,1) = 0.
                end if

                if(nvd_all(mode,2) > 0) then
                   vd_mean_all(mode,2) = 1e2*vd_mean_all(mode,2)/nvd_all(mode,2)
                else
                   vd_mean_all(mode,2) = 0.
                end if

                print '(i6,i4,2f10.4,2f12.6,2f12.6)', l, mode, r_mean_all(mode), 1e6*r_max_all(mode), &
                     vd_mean_all(mode,1), 1e2*vd_max_all(mode,1), vd_mean_all(mode,2), 1e2*vd_max_all(mode,2)
             end if
          enddo
          write(gol,*)  '________________________________________________________________________________' ; call goPr
       endif
       
    enddo LEVS
    !$OMP END PARALLEL
    
#endif /* M7 */

    status = 0
    
  END SUBROUTINE SEDIMENTATION_CALCV
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SEDIMENTATION_APPLY
  !
  ! !DESCRIPTION: Sedimentation calculation based on pre-calculated
  !               v_sedimentation.  
  !               Based on dynamw routine, including the slopes.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SEDIMENTATION_APPLY( region, status )
    !
    ! !USES:
    !
    use global_data, only : mass_dat, region_dat
    use meteodata  , only : spm_dat
    use chem_param,  only : ra
    use dims,        only : lm, at, bt, limits, sec_month, nsrce, tref, okdebug
    use chem_param,  only : mode_nm, mode_tracers, mode_tracers_sed,names, mode_nm_sed
    use partools,    only : par_reduce, par_reduce_element
#ifdef with_m7
    use emission_data, only : bb_lm
#ifndef without_emission
    use emission_data, only : emis_mass, emis_number
#endif
#endif
#ifdef with_budgets
    use emission_data, only : budemi_dat, sum_emission
    use ebischeme,     only : buddep_dat
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)               :: region
    !                                   
    ! !OUTPUT PARAMETERS:               
    !                                   
    integer, intent(out)              :: status
    !
    ! !REVISION HISTORY: 
    !   19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !  (1) effective only if m7 used...
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Sedimentation_Apply'

    real,dimension(:,:,:,:),pointer   :: rm,rxm,rym,rzm
    real,dimension(:,:,:),pointer     :: vsold
    real,dimension(:,:  ),pointer     :: vdold
    real,dimension(:,:,:),pointer     :: p   ! surface pressure
    real,dimension(:,:,:),allocatable :: f, pf, fx, fy, vs
    real,dimension(:,:,:),allocatable :: emit
    real,dimension(:,:),allocatable   :: vd 
    integer                           :: i, j, l, n, lmr, mode
    real, parameter                   :: one = 1.0
    real                              :: gamma, gam, l_gam
    real                              :: dp, dtime, month2dt
    real                              :: rmold, rmnew, rmplus, fdep, fsed
    integer,parameter                 :: iter_limit=200
    integer                           :: n_advz, iter
    integer                           :: nzone_v
    integer                           :: ls, le, ii, inmode
    real                              :: l_sum(3), g_sum(3)
    integer :: i1, j1, i2, j2, iflag, jflag, lflag

    !________________________start_________________________________
    ! Big if-condition over all routine

    ! start timing:
    call GO_Timer_Start( itim_appl, status )
    IF_NOTOK_RETURN(status=1)
    
#ifdef with_m7

    ! Bounds
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )    
    lmr = lm(region)
    
    allocate( f  ( i1:i2, j1:j2, lmr  ))    ! flux of tracer (kg)
    allocate(pf  ( i1:i2, j1:j2, lmr  ))    ! rxz  moment flux
    allocate(fx  ( i1:i2, j1:j2, lmr  ))    ! rmx flux
    allocate(fy  ( i1:i2, j1:j2, lmr  ))    ! rmy flux
    allocate(vs  ( i1:i2, j1:j2, lmr  ))
    allocate(vd  ( i1:i2, j1:j2       ))
    allocate(emit( i1:i2, j1:j2, bb_lm))

    call sedimentation_calcv(region, status) ! calculate sedimentation....
    IF_NOTOK_RETURN(status=1)

    call deposition_calcv(region, status)    ! calculate deposition.....
    IF_NOTOK_RETURN(status=1)
    
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(p)
    p   =>  spm_dat(region)%data
    rm  => mass_dat(region)%rm
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm

    if(okdebug) then
       write(gol,*) '    call sedimentation'; call goPr
    end if
    
    ! If requested, limit vertical slopes such that no non-negative tracer masses should occur
    ls = 1  ; le = lmr
    if (limits) then
       do n = 1, ubound (rzm, 4)
       do l = ls, le
       do j = j1, j2
       do i = i1, i2

          rzm(i,j,l,n) = max(min(rzm(i,j,l,n),rm(i,j,l,n)), -rm(i,j,l,n))

       end do
       end do
       end do
       end do
    endif

    ! determine timestep
    dtime=float(nsrce)/(2*tref(region))   

    ! conversion to emission per timestep
    month2dt=dtime/sec_month

    ! ================= 
    ! Loop over tracers
    ! ================= 
    do mode =1,nmod
       !do inmode=0,mode_nm(mode)
       do inmode=0,mode_nm_sed(mode)

          n = mode_tracers_sed(inmode,mode)
          !n = mode_tracers(inmode,mode)

          !------------- reset
          n_advz=1
          nullify(vsold)
          nullify(vdold)
          if (inmode == 0) then    ! number or mass tracer
             vsold => vn_sedimentation(region,mode)%d3 
             vdold => vn_deposition(region,mode)%surf
          else
             vsold => vm_sedimentation(region,mode)%d3 
             vdold => vm_deposition(region,mode)%surf
          endif

          !------------- find # iter needed for no CFL violation
          cfl: do

             vs(:,j1:j2,:) = vsold(:,j1:j2,:)/float(n_advz)
             vd(:,j1:j2)   = vdold(:,j1:j2)  /float(n_advz)

             advz: do iter=1,n_advz

                ! reset gamma
                l_gam = 0.

                do l=  1, lmr
                do j= j1, j2
                do i= i1, i2
                   ! vs (Pa) always downwards
                   ! calculate the flux at the bottom of box:
                   ! note that we define this as negative (vs>0, f<0)
                   ! pressure difference in layer dp (Pa)
                   dp =  (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1))   ! positive..
                   if(l /= 1) then
                      gamma=-vs(i,j,l)/dp  !always negative (according to negative cm dynamw)
                   else
                      gamma=-(vs(i,j,l)+vd(i,j))/dp  !always negative (according to negative cm dynamw)
                   endif
                   
                   ! replace 'one' with 0.99999 (NB/PLS)
                   if (abs(gamma)>=0.9999999) then
                      l_gam=max(l_gam,abs(gamma))
       !               exit advz   ! commented for consistency with TM5
                   end if
                   
                enddo
                enddo
                enddo

             enddo advz

             call Par_REDUCE( l_gam, 'MAX', gam, status, all=.true.)
             

             ! if CFL violation increase iteration counter and re-start cfl loop
             ! but check if the number of iterations becomes too large
             
             if(abs(gam)>=0.9999999) then  ! PLS: replace "one" with 0.9999999 to be consistent with above

!Commented since "exit advz" is commented, and hence not meaningfull                
!                 if(okdebug)then
!
!                    iflag=min(i,i2)
!                    jflag=min(j,j2)
!                    lflag=min(l,lmr)
!                   
!                    ! i-j-l information correct if from the processor where violation happened 
!                    print *,' --- CFL violation: gamma=',gam,vs(iflag,jflag,lflag),'    in (region,i,j,l,n)= ',region,iflag,jflag,lflag,n
!                    print*,' #iterations:', n_advz
!                 endif

                n_advz = n_advz + 1
                if ( n_advz > iter_limit ) then
                   status=1
                   IF_NOTOK_RETURN(status=1)
                end if
                cycle cfl
             else        ! situation OK, no CFL: use this n_advz
                exit cfl
             endif

          enddo cfl

          !------------- Apply number of iterations to calculate new tracer distributions

#ifdef with_budgets
          l_sum = 0.0
#endif

#ifndef without_emission
          if (n==ino3_a.or.n==inh4.or.n==imsa) then
             emit(:,j1:j2,:) = 0.0
          else if(inmode == 0) then
             emit(:,j1:j2,:) = 0.0
             do ii=1,mode_nm(mode)   ! add up all number emissions in the mode...
                emit(:,j1:j2,:) = emit(:,j1:j2,:) + emis_number(region,mode)%d4(:,j1:j2,:,ii)*month2dt/float(n_advz)
             enddo
          else ! this is a 'mass' emission with index nmode
             emit(:,j1:j2,:) = emis_mass(region,mode)%d4(:,j1:j2,:,inmode)*month2dt/float(n_advz)
          endif
#else
          emit = 0.0
#endif

          ! do the iteration
          do iter=1,n_advz

             ! compute f, pf, fx and fy
             do l=  1, lmr
             do j= j1, j2
             do i= i1, i2
                ! vs (Pa) always downwards
                ! calculate the flux at the bottom of box:
                ! note that we define this as negative (vs>0, f<0)
                ! pressure difference in layer dp (Pa)
                dp =  (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1))   ! positive..
                if(l /= 1) then
                   gamma=-vs(i,j,l)/dp  !always negative (according to negative cm dynamw)
                else
                   gamma=-(vs(i,j,l)+vd(i,j))/dp  !always negative (according to negative cm dynamw)
                endif
                f(i,j,l)  = gamma*(rm(i,j,l,n)-(one+gamma)*rzm(i,j,l,n))      !kg (neg.)
                pf(i,j,l) = -vs(i,j,l)*(gamma*gamma*rzm(i,j,l,n)-3.*f(i,j,l)) !neg.
                fx(i,j,l) = gamma*rxm(i,j,l,n)   !kg  (neg.)
                fy(i,j,l) = gamma*rym(i,j,l,n)   !kg  (neg.)

                ! - Cannot happen unless there is a bug (pls)
                !if (abs(gamma)>=one) then
                !   status=1
                !   IF_NOTOK_RETURN(status=1)
                !end if
                
             enddo
             enddo
             enddo

             ! calculate new tracer mass, and tracer mass slopes
             ! update rm, rzm, rxm and rym in interior layers of the column
             l = lmr
             do j = j1, j2
             do i = i1, i2

                dp =  (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1))
                rm (i,j,l,n)=rm(i,j,l,n) + f(i,j,l)
                rzm(i,j,l,n)=rzm(i,j,l,n)+  &
                             (        pf(i,j,l) &
                             -(-vs(i,j,l)) *rzm(i,j,l,n)            &
                             +3.*( ( -vs(i,j,l)) *rm (i,j,l,n)            &
                             -( f(i,j,l))*dp ))/dp
                   if(limits) then 
                      rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n))
                   endif
                   rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j,l))
                   rym(i,j,l,n)=rym(i,j,l,n)+(fy(i,j,l))
#ifdef with_budgets
                   nzone_v=nzon_vg(l) 
                   !if (n==40)write(1332,*) f(i,j,l)*1e3/ra(n)
                   buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) =   &
                   buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) -  f(i,j,l)*1e3/ra(n)   ! in mole
                   ! Downward flux is negative. And it is defined at the bottom of the box. A positive flux
                   ! is positive and would be an income for the respective grid cell. A minus sign takes care
                   ! We want to define the sedimentation as a cost.
#endif
                enddo
             enddo
             
             do l = lmr-1, 2, -1
             do j = j1, j2
             do i = i1, i2
                dp =  (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1))
                rmold = rm(i,j,l,n)
                if (l .le. bb_lm) then
                   rmplus = rmold + f(i,j,l)-f(i,j,l+1) + emit(i,j,l)
                else   
                   rmplus = rmold + f(i,j,l)-f(i,j,l+1) 
                endif
                rm(i,j,l,n) = rmplus
                rzm(i,j,l,n)=rzm(i,j,l,n)+  &
                           (        pf(i,j,l)-pf(i,j,l+1) &
                           -(vs(i,j,l+1)-vs(i,j,l)) *rzm(i,j,l,n)            &
                           +3.*( (-vs(i,j,l+1)-vs(i,j,l)) *rm (i,j,l,n)            &
                           -(f(i,j,l+1)+ f(i,j,l)) *dp ))/dp
                if(limits) then 
                   rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n))
                endif
                rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j,l)-fx(i,j,l+1))
                rym(i,j,l,n)=rym(i,j,l,n)+(fy(i,j,l)-fy(i,j,l+1))
#ifdef with_budgets
                nzone_v=nzon_vg(l)    
                buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) =   &
                buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) -  (f(i,j,l)-f(i,j,l+1))*1e3/ra(n)   ! in mole
                ! The sedimentation is calculated as 'income' again. With the minus sign, those are costs.
#ifndef without_emission
                if ( l <= bb_lm ) then
                   budemi_dat(region)%emi(i,j,nzone_v,n) =   &
                        budemi_dat(region)%emi(i,j,nzone_v,n) + emit(i,j,l)*1e3/ra(n)   ! in mole
                end if
#endif
#endif
             enddo
             enddo
             enddo
             
             l = 1
             do j = j1, j2
             do i = i1, i2
                dp =  (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1))
                if(vd_in_deposition) then
                   write(gol,*)' apply_sedimentation: vs in deposition disabled'; call goBug
                   status=1
                   IF_NOTOK_RETURN(status=1)
                end if
                
                !if(vd_in_deposition) then
                !   rm (i,j,l,n)=rm(i,j,l,n) - f(i,j,l+1)
                !   rzm(i,j,l,n)=rzm(i,j,l,n)+  &
                !       (       -pf(i,j,l+1) &
                !              -(vs(i,j,l+1)) *rzm(i,j,l,n)            &
                !        +3.*( (-vs(i,j,l+1)) *rm (i,j,l,n)            &
                !               -(f(i,j,l+1))*dp ))/dp
                !   if(limits) then 
                !      rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n))
                !   endif
                !   rxm(i,j,l,n)=rxm(i,j,l,n)+(-fx(i,j,l+1))
                !   rym(i,j,l,n)=rym(i,j,l,n)+(-fy(i,j,l+1))
                !else
                !                   rm (i,j,l,n)=rm(i,j,l,n) + f(i,j,l)-f(i,j,l+1)
                !                   rzm(i,j,l,n)=rzm(i,j,l,n)+  &
                !                       (        pf(i,j,l)-pf(i,j,l+1) &
                !                              -(vs(i,j,l+1)-vs(i,j,l)) *rzm(i,j,l,n)            &
                !                        +3.*( (-vs(i,j,l+1)-vs(i,j,l)) *rm (i,j,l,n)            &
                !                               -(f(i,j,l+1)+ f(i,j,l))*dp ))/dp
                !                   if(limits) then 
                !                      rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n))
                !                   endif

                ! apply to rxm, rym the fluxes:
                rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j,l)-fx(i,j,l+1))
                rym(i,j,l,n)=rym(i,j,l,n)+(fy(i,j,l)-fy(i,j,l+1))

                ! for rm: apply emissions, sedimentation flux from above
                ! deposition and sedimentation at surface: Backward Eularian
                rmold = rm(i,j,l,n)
                rmplus = (rmold - f(i,j,l+1) + emit(i,j,l))   ! note f is negative !
                rmnew = rmplus/(1. + (vd(i,j) + vs(i,j,1))/dp)
                rm(i,j,l,n) = rmnew
                if(rmold > 1e-8) rzm(i,j,l,n) = rzm(i,j,l,n)*rmnew/rmold

                ! budget:
#ifdef with_budgets
                if((vs(i,j,1) + vd(i,j)) > 1e-12) then
                   fsed =  -(rmplus-rmnew)*vs(i,j,1)/(vs(i,j,1) + vd(i,j))   ! is negative
                   fdep =  -(rmplus-rmnew)*vd(i,j)  /(vs(i,j,1) + vd(i,j))
                else
                   fsed = 0.0
                   fdep = 0.0
                endif
                
                if(n == 1) then 
                   l_sum(1) = l_sum(1) + f(i,j,l+1) - fsed ! goes into sum_sedimentation
                   l_sum(2) = l_sum(2) - fdep              ! goes into sum_drydep
#ifndef without_emission
                   l_sum(3) = l_sum(3) + emit(i,j,1)       ! goes into sum_emission
#endif
                endif
                nzone_v=nzon_vg(l)    
                buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) =   &
                buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) - (fsed - f(i,j,l+1) )*1e3/ra(n)   ! in mole
                ! The comment line says that fsed is negative. As we regard sedimentation as a cost, we define
                ! a positive sedimentation as a loss of material. A negative fsed is loss of material, so we
                ! need to have -fsed. The (upward) flux on the top of the grid cell is a cost as well. 
                buddep_dat(region)%dry(i,j,n) =   &
                     buddep_dat(region)%dry(i,j,n) - fdep*1e3/ra(n)   ! in mole
#ifndef without_emission
                budemi_dat(region)%emi(i,j,1,n) =   &
                     budemi_dat(region)%emi(i,j,1,n) + emit(i,j,1)*1e3/ra(n)   ! in mole
#endif
#endif /* BUDGETS */
                
             enddo !i
             enddo !j

          enddo  ! iter !

#ifdef with_budgets
          call PAR_REDUCE_ELEMENT( l_sum, 'SUM', g_sum, status )

          if(isRoot ) then 
             sum_sedimentation(region) = sum_sedimentation(region) + g_sum(1)
             sum_drydep       (region) = sum_drydep       (region) + g_sum(2)
#ifndef without_emission
             sum_emission     (region) = sum_emission     (region) + g_sum(3)
#endif
          end if
#endif

       enddo ! loop over tracers in mode
    enddo ! loop over modes

    deallocate( f)
    deallocate(pf)
    deallocate(fx)
    deallocate(fy)
    deallocate(vs)
    deallocate(vd)
    deallocate(emit)

    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(vsold)
    nullify(vdold)

#endif /* M7 */

    ! end timing:
    call GO_Timer_End( itim_appl, status )
    IF_NOTOK_RETURN(status=1)

    status = 0

  END SUBROUTINE SEDIMENTATION_APPLY
  !EOC

#ifdef with_m7
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    FALL1
  !
  ! !DESCRIPTION: function to calculate the fall velocity from the particle 
  !               diameter, as a function of density, temperature and pressure.
  !\\
  !\\
  ! !INTERFACE:
  !
  REAL FUNCTION FALL1( press, zmd, zvis, t, zdensair, densaer_p) result(vt)
    !
    ! !INPUT PARAMETERS:
    !
    real, intent(in)      :: press      ! pressure in Pa (or bar?)
    real, intent(in)      :: zmd        ! median radius of aerosol (m)
    real, intent(in)      :: zvis       ! viscosity  (kg/(sm))
    real, intent(in)      :: t          ! temperature (K)
    real, intent(in)      :: zdensair   ! density air (kg/m3)
    real, intent(in)      :: densaer_p  ! density aerosol particles (kg/m3)
    !
    ! !REVISION HISTORY: 
    !    2 Apr 2010 - P. Le Sager - 
    !
    ! !REMARKS:
    ! Called in Sedimentation_Apply, only if m7 used.
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    real :: zlair
    
    ! ----------- start
    if (zmd > tiny(zmd)) then 
       vt=2./9.*(densaer_p-zdensair) * grav/zvis*(zmd/2.)**2.

       zlair=0.066*(1.01325e5/press)*t/293.15*1e-6

       !--- With cunnigham slip- flow correction (S&P, Equation 8.34):

       vt = vt * (1.+ 1.257*zlair/zmd*2. +  0.4*zlair/zmd*2. *   exp(-1.1/(zlair/zmd*2.)) )
    else
       vt = 0.0
    endif

  END FUNCTION FALL1
  !EOC
  
#endif /* M7 */  

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    AEROSOL_RADIUS
  !
  ! !DESCRIPTION: function that calculates the median aerosol radius (m),
  !               given the total mass and number of a log-normal distribution.   
  !\\
  !\\
  ! !INTERFACE:
  !
  REAL FUNCTION AEROSOL_RADIUS(mtot, ntot, sigma, rho_aer) result(radius)
    !
    ! !USES:
    !
    use Binas, only: Pi
    !
    ! !INPUT PARAMETERS:
    !
    real, intent(in)   :: mtot    ! total mass (kg)
    real, intent(in)   :: ntot    ! total number (#)
    real, intent(in)   :: sigma   ! the sigma of the log-normal aerosol size distribution
    real, intent(in)   :: rho_aer ! the density of the aerosol (kg/m3)
    !
    ! !REVISION HISTORY: 
    !    2 Apr 2010 - P. Le Sager - 
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    real               :: lns

    if(mtot >  tiny(mtot) .and. ntot > tiny(ntot)) then
       lns = alog(sigma)
       radius = (mtot/(rho_aer*4.0*pi/3.0 *ntot *exp((9./2.)*lns**2)))**(1./3.)
    else
       radius = 0.0
    endif

  END FUNCTION AEROSOL_RADIUS
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALCULATE_RH
  !
  ! !DESCRIPTION: calculate the rh, with 100% rh assumption in cloudy part.
  !               In the cloud free part, the rh thus is smaller, and the water
  !               uptake by aerosols will be smaller.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CALCULATE_RH
    !
    ! !USES:
    !
    use chem_param, only : xmh2o
    use meteodata,  only : phlb_dat
    use MeteoData,  only : temper_dat, humid_dat, cc_dat
    use dims,       only : nregions,   im, jm ,lm
    !
    ! !REVISION HISTORY: 
    !   19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:), pointer    :: phlb   ! pressure at border (Pa)
    real, dimension(:,:,:), pointer    :: t      ! temperature (K)
    real, dimension(:,:,:), pointer    :: q      ! humidity (kg h2o / kg air)
    real, dimension(:,:,:), pointer    :: cc     ! cloud cover (0-1)
    real                               :: tr, wv, airn, h2on, rrh, s, ccs
    integer                            :: region, i,j,l, i1,i2,j1,j2

    do region = 1, nregions

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       nullify(phlb)
       nullify(t)
       nullify(q)
       nullify(cc)

       phlb => phlb_dat(region)%data
       t    => temper_dat(region)%data
       q    => humid_dat(region)%data
       cc   => cc_dat(region)%data


       do l=1, lm(region)
          do j=j1, j2
             do i=i1, i2
                tr = 1. - 373.15/t(i,j,l)
                wv=exp((((-.1299*tr-.6445)*tr-1.976)*tr+13.3185)*tr)
                airn = 7.24e16*0.5*(phlb(i,j,l) + phlb(i,j,l+1))/t(i,j,l)   ! somethings seem redundent here!
                h2on = q(i,j,l)*airn*xmair/xmh2o                            ! leave it for readability!
                rrh = h2on*t(i,j,l)/(1013.25*wv*7.24e16)
                s = 0.01*max(0.0, min(rrh, 99.9 ) )       ! 0-0.999 scale!
                rha(region)%d3(i,j,l) = s
                ! scale relative humidity to cloudfree part
                ! assuming 100% rh in the cloudy part, but never smaller than 0.75!
                if (s > 0.90) then 
                   ccs = cc(i,j,l)
                   if((1. - ccs) > tiny(ccs)) s = max(0.75, (s-ccs)/(1. - ccs)) 
                endif
                rh(region)%d3(i,j,l) = s
             enddo
          enddo
       enddo

       nullify(phlb)
       nullify(t)
       nullify(q)
       nullify(cc)

    enddo

  END SUBROUTINE CALCULATE_RH
  !EOC

END MODULE SEDIMENTATION