init_dst.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE init_dst (DST, U, &
                     hprogram,  &    ! Program calling unit
                  ksize_nature_p, & ! Number of nature points in a patch
                  kr_nature_p, &    ! Mask from patch --> nature vectors
                  kpatch, &         ! Maximum number of patches
                  pvegtype_patch  ) ! fraction (in a nature point) of a vegtype for a patch

!
USE modd_dst_n, ONLY : dst_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_dst_surf
USE modd_data_cover_par, ONLY : nvt_no, nvt_rock
!
USE modi_get_luout
USE modi_get_vegtype_2_patch_mask
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!PASSED VARIABLES
!
TYPE(dst_t), INTENT(INOUT) :: dst
TYPE(surf_atm_t), INTENT(INOUT) :: u
!
 CHARACTER(LEN=6), INTENT(IN)      :: hprogram              !Passing unit
!
INTEGER, DIMENSION(:), POINTER :: ksize_nature_p
INTEGER, DIMENSION(:,:), POINTER :: kr_nature_p
INTEGER, INTENT(IN)  :: kpatch
REAL, DIMENSION(:,:,:), POINTER :: pvegtype_patch
!
!LOCAL VARIABLES
 CHARACTER(LEN=4)    :: crgunit               ! type of log-normal geometric mean radius
INTEGER             :: jveg                  ! Counter for vegetation classes
INTEGER             :: jveg_in               ! Vegetation index
INTEGER             :: jpatch                ! Counter for patches
INTEGER             :: jmode                 ! Counter for dust modes
INTEGER             :: jmode_idx             ! Index for dust modes
INTEGER             :: iluout
INTEGER             :: isize_largest_dst
REAL(KIND=JPRB) :: zhook_handle

!get output listing unit
IF (lhook) CALL dr_hook('INIT_DST',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
!Allocate memory for the real values which will be used by the model
ALLOCATE(dst%XEMISRADIUS_DST(ndstmde))
ALLOCATE(dst%XEMISSIG_DST   (ndstmde))
ALLOCATE(dst%XMSS_FRC_SRC   (ndstmde))
!
!Get initial size distributions. This is cut and pasted
!from dead routine dstpsd.F90
!Check for different source parameterizations
IF(cemisparam_dst.EQ."She84")THEN
  crgunit   = 'MASS'
  xemisradius_ini_dst(:) = 0.5d6 * (/ 0.0111e-6, 2.524e-6, 42.10e-6 /)  ! [um]  Mass median radius She84 p. 75 Table 1
  xemissig_ini_dst(:) = (/ 1.89 , 2.0  , 2.13 /)             ! [frc] Geometric standard deviation She84 p. 75 Table 1
  xmss_frc_src_ini(:) = (/2.6e-6, 0.781, 0.219/)             ! [frc] Mass fraction She84 p. 75 Table 1
ELSEIF(cemisparam_dst.EQ."PaG77")THEN
  crgunit   = 'MASS'
  xemisradius_ini_dst(:) = 0.5d6 * (/0.27e-6  ,  5.6e-6  ,  57.6e-6 /) ! [um] Mass median radius PaG77 p. 2080 Table 1 
  xemissig_ini_dst(:) = (/ 1.88, 2.2  , 1.62 /)              ! [frc] Geometric standard deviation PaG77 p. 2080 Table 1
  xmss_frc_src_ini(:) = (/0.036, 0.957, 0.007/)              ! [frc] Mass fraction BSM96 p. 73 Table 2 (ad hoc)  
ELSEIF(cemisparam_dst.EQ."Dal87") THEN 
  ! D'Almeida, 1987 as default
  crgunit   = 'MASS'
  xemisradius_ini_dst(:) = 0.5d6 * (/ 0.832e-6 ,  4.82e-6 , 19.38e-6 /) ! [um] Mass median radius BSM96 p. 73 Table 2
  xemissig_ini_dst(:) = (/ 2.10, 1.90 , 1.60 /)              ! [frc] Geometric standard deviation BSM96 p. 73 Table 2
  xmss_frc_src_ini(:) = (/0.036, 0.957, 0.007/)              ! [frc] Mass fraction BSM96 p. 73 Table 2  
ELSEIF (cemisparam_dst.EQ."alf98".OR.cemisparam_dst.eq."EXPLI") THEN  !  Alfaro et al 1998 as default
  IF (cemisparam_dst.EQ."alf98") xflx_mss_fdg_fct = 6e-4
  IF (cemisparam_dst.EQ."EXPLI") xflx_mss_fdg_fct = 3.5e-4
  crgunit   = 'MASS'
  xemisradius_ini_dst(:) = 0.5*(/ 1.5, 6.7, 14.2 /)      ! [um] Mass median radius BSM96 p. 73 Table 2
  xemissig_ini_dst(:) = (/1.70, 1.60, 1.50/)          ! [frc] Geometric standard deviation BSM96 p. 73 Table 2
  xmss_frc_src_ini(:) = (/0.01, 0.19, 0.8 /)          ! [frc] Mass fraction BSM96 p. 73 Table 2  
ELSEIF (cemisparam_dst.EQ."AMMA ") THEN ! Default: New distribution from AMMA
  xflx_mss_fdg_fct = 105.e-4
  crgunit   = 'NUMB' 
  xemisradius_ini_dst(:) = 0.5*(/ 0.078, 0.641, 5.00 /)  ! [um] Number median radius 
  xemissig_ini_dst(:) = (/ 1.75,  1.76, 1.70/)        ! [frc] Geometric standard deviation 
  xmss_frc_src_ini(:) = (/0.008, 0.092, 0.99/)        ! [frc] Mass fraction  
ELSEIF (cemisparam_dst.EQ."CRUM ") THEN ! Default: New distribution from AMMA
  xflx_mss_fdg_fct = 10.e-4
  crgunit   = 'NUMB' 
  xemisradius_ini_dst(:) = 0.5*(/ 0.078, 0.641, 5.00 /)  ! [um] Number median radius 
  xemissig_ini_dst(:) = (/  1.75,   1.76, 1.70   /)   ! [frc] Geometric standard deviation   
  xmss_frc_src_ini(:) = (/0.0005, 0.0029, 0.9966/)    ! [frc] Mass fraction  
ELSE
  WRITE(iluout,*) " FATAL ERROR "
  WRITE(iluout,*) " YOU MUST DECIDE THE EMISSIUON PARAMETERIZATION, YOU USES "
  WRITE(iluout,*) " CEMISPARAM_DST = ",cemisparam_dst," AND IT IS NOT DEFINED "
  WRITE(iluout,*) " see init_dstn.f90 to see what dust parameterization is available. "
  CALL abor1_sfx("INIT_DST: CEMISPARAM_DST not defined ")
ENDIF
!
DO jmode=1,ndstmde
  jmode_idx=jorder_dst(jmode)
  !
  dst%XEMISSIG_DST   (jmode) = xemissig_ini_dst(jmode_idx)
  dst%XEMISRADIUS_DST(jmode) = xemisradius_ini_dst(jmode_idx)
  dst%XMSS_FRC_SRC   (jmode) = xmss_frc_src_ini(jmode_idx)
  !
  !Get emisradius, and at the same time convert to number median radius
  IF (crgunit=='MASS') &
    dst%XEMISRADIUS_DST(jmode) = dst%XEMISRADIUS_DST(jmode) * exp(-3.d0 * (log(dst%XEMISSIG_DST(jmode)))**2)  
  !
ENDDO
!
!Normalize the sum of the emissions to 1 so that all dust is
!put in one mode or the other
IF(sum(dst%XMSS_FRC_SRC(:)).LT.1.) dst%XMSS_FRC_SRC(:) = dst%XMSS_FRC_SRC(:) / sum(dst%XMSS_FRC_SRC(:))
!
!Allocate memory
!ALLOCATE(NVEGNO_DST)
!Set the number of classes that can emit dust (fxm: set this elsewhere)
nvegno_dst = 2
!
!Allocate memory for the vegtype-translator
ALLOCATE(dst%NVT_DST(nvegno_dst))
!
!Set the dust/vegtype translator vector
dst%NVT_DST(1)  = nvt_no
dst%NVT_DST(2)  = nvt_rock
!
!Allocate memory for roughness lengths of erodible surfaces
ALLOCATE(dst%Z0_EROD_DST(nvegno_dst))
!
!Set the roughness lengths corresponding to erodible surfaces
!Smooth roughness length is given to 1.d-5 (dstmbl.f90)
dst%Z0_EROD_DST(1) = 30.d-6    !m (30 um) 
dst%Z0_EROD_DST(2) = 200.d-6   !m (200 um) 
!
!Allocate memory for dust emitter surface vectors in patch vectors
IF (.NOT.ASSOCIATED(dst%NSIZE_PATCH_DST)) ALLOCATE(dst%NSIZE_PATCH_DST(nvegno_dst,kpatch))
!
DO jpatch = 1,kpatch
  DO jveg = 1,nvegno_dst
    !Count all the points in the patch where you have dust emitter vegetation
    dst%NSIZE_PATCH_DST(jveg,jpatch) = count(pvegtype_patch(:,dst%NVT_DST(jveg),jpatch) > 0.) 
  ENDDO
ENDDO
!
!Find the largest dust emitter vector in any patch
!ALLOCATE (NSIZE_LARGEST_DST)
isize_largest_dst = 0
DO jpatch=1,kpatch
  DO jveg = 1,nvegno_dst
    isize_largest_dst = max(isize_largest_dst,dst%NSIZE_PATCH_DST(jveg,jpatch))
  ENDDO
ENDDO
!
!Allocate memory for NR_PATCH_DST mask translate from patch vector to dust vector
ALLOCATE(dst%NR_PATCH_DST(isize_largest_dst,nvegno_dst,kpatch))
!
!Initialize the mask array
dst%NR_PATCH_DST(:,:,:)=0
!
!Get values from the dust emitter vegetation mask
DO jpatch=1,kpatch
  DO jveg=1,nvegno_dst
    jveg_in = dst%NVT_DST(jveg)          ! Get the real vegtype index
#ifdef RJ_OFIX
    CALL get_vegtype_2_patch_mask(iluout,    &
             dst%NSIZE_PATCH_DST(jveg,jpatch),             &!I Size of dust emitter vector
             ksize_nature_p(jpatch),                   &!I Size of patch vector
             u%NSIZE_NATURE,                             &!I Size of nature vector
!RJ: attempt to make this call generic
             kr_nature_p(:ksize_nature_p(jpatch),jpatch),&!I Mask from patch to nature
             pvegtype_patch,                           &!I Fraction of vegtype of nature point within jpatch 
             dst%NR_PATCH_DST(:dst%NSIZE_PATCH_DST(jveg,jpatch),jveg,jpatch),  &!O Part of mask array to fill with values
             kpatch,                                   &!I Number of possible patches
             jpatch,                                   &!I Index of patch in question
             jveg_in                                  &!I Index of vegtype in question
             )  
#else
    CALL get_vegtype_2_patch_mask(iluout,    &
             dst%NSIZE_PATCH_DST(jveg,jpatch),             &!I Size of dust emitter vector
             ksize_nature_p(jpatch),                   &!I Size of patch vector
             u%NSIZE_NATURE,                             &!I Size of nature vector
             kr_nature_p,                              &!I Mask from patch to nature
             pvegtype_patch,                           &!I Fraction of vegtype of nature point within jpatch 
             dst%NR_PATCH_DST(:dst%NSIZE_PATCH_DST(jveg,jpatch),jveg,jpatch),  &!O Part of mask array to fill with values
             kpatch,                                   &!I Number of possible patches
             jpatch,                                   &!I Index of patch in question
             jveg_in                                  &!I Index of vegtype in question
             )  
#endif
  ENDDO !Loop on patches
ENDDO    !Loop on veg-types
!
IF (lhook) CALL dr_hook('INIT_DST',1,zhook_handle)
!
END SUBROUTINE init_dst