coupling_dstn.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 coupling_dst_n (DSTK, KK, PK, PEK, DK, &
       HPROGRAM,                 &!I [char] Type of ISBA version
       KI,                       &!I [nbr] number of points in patch
       KDST,                     &!I Number of dust emission variables
       PPS,                      &!I [Pa] surface pressure
       PQA,                      &!I [kg/kg] atmospheric specific humidity
       PRHOA,                    &!I [kg/m3] atmospheric density
       PPA,                      &!I [K] Atmospheric pressure
       PTA,                      &!I [K] Atmospheric temperature
       PU,                       &!I [m/s] zonal wind at atmospheric height 
       PUREF,                    &!I [m] reference height of wind
       PV,                       &!I [m/s] meridional wind at atmospheric height
       PZREF,                    &!I [m] reference height of wind
       PSFDST                    &!O [kg/m2/sec] flux of dust
       )  
  
!PURPOSE
!-------
!Recieve a full vector for a given patch containing dust emitter surface points
!Pack it to vectors for which are purely dust emitter surfaces
!Send back the vector PSFDST (kg/m2/s) coming from this patch.
!In COUPLING_ISBA this will again be weighted by fraction_of_patch_in_nature_point

!AUTHOR
!-------
!ALF GRINI <alf.grini@cnrm.meteo.fr>  2005
! Modification P. Tulet introduce friction velocity for mode repartition upon 
! Alfaro et al, 1998 (Geo. Res. Lett.)
!
!!      Modified    09/2012  : J. Escobar , SIZE(PTA) not allowed without-interface , replace by KI
!
!
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
USE modd_diag_n, ONLY : diag_t
!
USE modd_dst_n, ONLY : dst_t
!
USE modd_csts, ONLY : xrd                      ! [J/K/kg] The universal gas constant  
       

USE modd_dst_surf

USE modi_surface_cd
USE modi_cls_wind
USE modi_cls_t

USE modi_dustflux_get     ! Dust mobilization routines
USE modi_dustflux_get_mb  ! Dust mobilization routines (M. Mokhtari)  
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE

!INPUT
!
TYPE(dst_t), INTENT(INOUT) :: DSTK
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(diag_t), INTENT(INOUT) :: DK
!
 CHARACTER(LEN=*), INTENT(IN)       :: HPROGRAM       !I Name of program
INTEGER, INTENT(IN)                :: KI             !I Number of points in patch
INTEGER, INTENT(IN)                :: KDST           !I Number of dust emission variables
REAL, DIMENSION(KI), INTENT(IN)    :: PPS            !I [Pa] surface pressure
REAL, DIMENSION(KI), INTENT(IN)    :: PQA            !I [kg/kg] atmospheric specific humidity
REAL, DIMENSION(KI), INTENT(IN)    :: PRHOA          !I [kg/m3] atmospheric density
REAL, DIMENSION(KI), INTENT(IN)    :: PPA            !I [K] Atmospheric pressure
REAL, DIMENSION(KI), INTENT(IN)    :: PTA            !I [K] Atmospheric temperature
REAL, DIMENSION(KI), INTENT(IN)    :: PU             !I [m/s] zonal wind at atmospheric height 
REAL, DIMENSION(KI), INTENT(IN)    :: PUREF          !I [m] reference height of wind
REAL, DIMENSION(KI), INTENT(IN)    :: PV             !I [m/s] meridional wind at atmospheric height
REAL, DIMENSION(KI), INTENT(IN)    :: PZREF          !I [m] reference height of wind
REAL, DIMENSION(KI,KDST), INTENT(OUT) :: PSFDST      !O [kg/m2/sec] flux of dust for a patch
  
!LOCAL VARIABLES
REAL, DIMENSION(KI,NVEGNO_DST,NDSTMDE) :: ZSFDST_TILE  ![kg/m2] flux of dust for each vegetation types and each mode
INTEGER                            :: JVEG           ![idx] counter for vegetation types
REAL(KIND=JPRB) :: ZHOOK_HANDLE


!Initialize output which is total flux of dust (kg/m2/sec) from this patch
!This number will get weighted by fraction of patch later on (in coupling_isban)
IF (lhook) CALL dr_hook('COUPLING_DST_N',0,zhook_handle)
psfdst(:,:)=0.d0

!Initiate dust emissions for all points in patch
zsfdst_tile(:,:,:)=0.d0

DO jveg=1,nvegno_dst

  !Jump out of loop if no dust emitter points
  IF (dstk%NSIZE_PATCH_DST(jveg)==0) cycle

   CALL treat_surf(  &
      dstk%NSIZE_PATCH_DST(jveg),  &!I[idx] number of dust emitter points in patch
      dstk%NR_PATCH_DST (:,jveg)   &!I[idx] index translator from patch to dustemitter
            )  
   
ENDDO  !Loop on dust emitter vegetation
  
!Average dust flux from the two vegetation "tiles"
!Make sure you do this correctly so that area is OK.
!Need to have this as kg/m2/sec from PATCH because
!afterwards this is weighted by fraction of patch
 CALL avg_flux_dst(nvegno_dst)
  
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
IF (lhook) CALL dr_hook('COUPLING_DST_N',1,zhook_handle)
!
CONTAINS
    
SUBROUTINE treat_surf(KSIZE,KMASK)
!
USE modd_diag_n, ONLY : diag_t
!
IMPLICIT NONE

INTEGER, INTENT(IN)    :: KSIZE  !Size of dust emitter vector
INTEGER, DIMENSION(:), INTENT(IN)    :: KMASK !mask from patch vector to dust emitter vector

!ALLOCATABLE VARIABLES FOR ALL VEGETATION TYPES
REAL, DIMENSION(KSIZE) :: ZP_CLAY       ![frc] fraction of clay
REAL, DIMENSION(KSIZE) :: ZP_SAND       ![frc] fraction of sand
!
REAL, DIMENSION(KSIZE) :: ZP_PA         ![K] Atmospheric pressure
REAL, DIMENSION(KSIZE) :: ZP_TA         ![K] Atmospheric temperature
REAL, DIMENSION(KSIZE) :: ZP_QA         ![kg_{H2O}/kg_{air}] specific humidity
REAL, DIMENSION(KSIZE) :: ZP_RHOA_2M    ![kg/m3] atmospheric density close to surface
REAL, DIMENSION(KSIZE) :: ZP_SFMER      ![m/s] wind friction in meridional direction
REAL, DIMENSION(KSIZE) :: ZP_SFZON      ![m/s] wind friction in zonal direction
REAL, DIMENSION(KSIZE) :: ZP_U          ![m/s] zonal wind at atmospheric height
REAL, DIMENSION(KSIZE) :: ZP_V          ![m/s] meridional wind at atmospheric height
REAL, DIMENSION(KSIZE) :: ZP_VMOD       ![m/s] Wind at atmospheric height
!
REAL, DIMENSION(KSIZE) :: ZP_SFDST      ![kg/m2/s] dust flux (kg/m2/s)
REAL, DIMENSION(KSIZE,NDSTMDE) :: ZP_SFDST_MDE ![kg/m2/s] dust flux from modes
!
REAL, DIMENSION(KSIZE) :: ZP_TG         ![K] ground temperature
REAL, DIMENSION(KSIZE) :: ZP_WG         ![m3/m3] ground volumetric water content
REAL, DIMENSION(KSIZE) :: ZP_WSAT       ![m3/m3] saturation volumetric water content
REAL, DIMENSION(KSIZE) :: ZP_USTAR      ![m/s] wind friction speed
!
REAL, DIMENSION(KSIZE) :: ZP_Z0_EROD    ![m] roughness length for erodible surface
REAL, DIMENSION(KSIZE) :: ZP_DST_EROD      !Erodable Surface (COVER004=1 and COVER005=0.01)
REAL, DIMENSION(KSIZE,3) :: ZP_MSS_FRC_SRC ![%] dust mode fraction of emitted mass 
REAL, DIMENSION(KSIZE) :: ZP_INTER ![%] dust mode fraction of emitted mass 
REAL, DIMENSION(KSIZE) :: ZP_CD_DST     ! drag coefficient uses by DEAD
REAL, DIMENSION(KSIZE) :: ZH               ! 10m (wind altitude)
!
REAL, DIMENSION(KSIZE) :: ZP_MER10M     ![m/s] meridional wind at 10m
REAL, DIMENSION(KSIZE) :: ZP_WIND10M    ![m/s] wind at 10m
REAL, DIMENSION(KSIZE) :: ZP_ZON10M     ![m/s] zonal wind at 10m
REAL, DIMENSION(KSIZE) :: ZP_PS         ![Pa] surface pressure
REAL, DIMENSION(KSIZE) :: ZP_TS         ![K] surface temperature
REAL, DIMENSION(KSIZE) :: ZP_Q2M        ![kg_{H2O}/kg_{air}] speficic humidity
REAL, DIMENSION(KSIZE) :: ZP_HU2M       ![-] relative humidity
REAL, DIMENSION(KSIZE) :: ZP_RHOA       ![kg/m3] atmospheric density
REAL, DIMENSION(KSIZE) :: ZP_T2M        ![K] 2M temperature
REAL, DIMENSION(KSIZE) :: ZP_UREF       ![m] reference height for wind
REAL, DIMENSION(KSIZE) :: ZP_ZREF       ![m] reference height for wind
REAL, DIMENSION(KSIZE) :: ZP_CD         !I [] drag coefficient for momentum from isba
REAL, DIMENSION(KSIZE) :: ZP_CH         !I [] drag coefficient for heat
REAL, DIMENSION(KSIZE) :: ZP_CDN        ![-] drag coefficient neutral atm
REAL, DIMENSION(KSIZE) :: ZP_RI         !I [-] Richardson number
REAL, DIMENSION(KSIZE) :: ZP_Z0H        ![frc] Z0 heat with snow 
!
REAL, DIMENSION(5)   :: ZSEUIL
REAL, DIMENSION(6,2) :: ZPCEN
INTEGER                :: JJ, JS            !Counter for vector points
INTEGER                :: JMODE         
REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('COUPLING_DST_n:TREAT_SURF',0,zhook_handle)

!Pack patch-vectors to dust emitter vectors
!i.e. allocate memory for the packed (dust emitter) vectors
DO jj=1,ksize
  zp_clay(jj) = kk%XCLAY(kmask(jj),1)  
  zp_ps(jj) = pps(kmask(jj))
  zp_ts(jj) = dk%XTS  (kmask(jj))
  zp_qa(jj) = pqa(kmask(jj))
  zp_rhoa(jj) = prhoa(kmask(jj))
  zp_sand(jj) = kk%XSAND(kmask(jj),1)
  zp_pa(jj) = ppa(kmask(jj)) 
  zp_ta(jj) = pta(kmask(jj)) 
  zp_tg(jj) = pek%XTG (kmask(jj),1)
  zp_u(jj) = pu(kmask(jj))
  zp_uref(jj) = puref(kmask(jj)) 
  zp_v(jj) = pv(kmask(jj))
  zp_wg(jj) = pek%XWG (kmask(jj),1)
  zp_wsat(jj) = kk%XWSAT(kmask(jj),1)                      
  zp_zref(jj) = pzref(kmask(jj))    
  zp_cd(jj) = dk%XCD  (kmask(jj))
  zp_cdn(jj) = dk%XCDN (kmask(jj))
  zp_ch(jj) = dk%XCH  (kmask(jj))
  zp_ri(jj) = dk%XRI  (kmask(jj))
  zp_z0h(jj) = dk%XZ0H (kmask(jj))
ENDDO
!
!Manipulate some variables since we assume dust emission occurs over flat surface
zp_z0_erod(:) = dstk%Z0_EROD_DST(jveg)   !Set z0 to roughness of erodible surface
!
IF (jveg ==1) THEN
  zp_dst_erod(:) = 1.
ELSE
  zp_dst_erod(:) = 0.01
ENDIF
!
zp_cd_dst(:) = zp_cd(:)
!
IF (cvermod/='CMDVER') THEN
  !Re-calculate the drag over the dust emitter (erodable) surface. Since we
  !don't use the roughness length of the patch, but rather use specific roughness
  !lengths of dust emitter surfaces, the drag changes.  
  CALL surface_cd(zp_ri, zp_zref, zp_uref, zp_z0_erod, zp_z0h, zp_cd_dst, zp_cdn)
ENDIF
!
!Get total wind speed
zp_vmod(:) = sqrt(zp_u(:)**2 + zp_v(:)**2)  
!
zp_ustar(:) =  sqrt(zp_cd_dst(:))*zp_vmod
!Get zonal friction speed (m/s)
zp_sfzon(:) = - sqrt(zp_cd_dst(:))*zp_u(:)
!Get meridional surface stress
zp_sfmer(:) = - sqrt(zp_cd_dst(:))*zp_v(:)
!
!Get the 10m wind speed (needed for one operation in dust model)
!And get the density at 2m to feed into erosion model
zh(:)=10.
 CALL cls_wind(zp_u, zp_v, zp_uref, zp_cd_dst, zp_cdn, zp_ri, zh, &
                zp_zon10m, zp_mer10m                 )
zh(:)=2.
 CALL cls_t(zp_ta, zp_qa, zp_pa, zp_ps, zp_zref, zp_cd_dst, zp_ch, zp_ri,   &
              zp_ts, zp_z0h, zh, zp_t2m  )

!Get density at 2 meters rho=P/(RT)
zp_rhoa_2m(:) = zp_ps(:)/(zp_t2m(:)*xrd)
!
!Get wind speed at 10m heigh
zp_wind10m(:) = sqrt(zp_zon10m(:)**2 + zp_mer10m(:)**2)
!     
!the erodible surface
zp_wind10m(:) = max(zp_wind10m(:), 1e-2)
!
IF (cvermod=='CMDVER') THEN
  !
  CALL dustflux_get_mb(               &
       zp_ustar,                      & !I [m/s] wind friction speed over erodible surface
       zp_rhoa_2m,                    & !I [kg/m3] air density at surface
       zp_wg,                         & !I [m3/m3] volumetric water content of ground
       zp_z0_erod,                    & !I [m] roughness length of erodible surface
       zp_wsat,                       & !I [m3/m3] saturation water volumetric content
       zp_clay,                       & !I [frc] mass fraction of clay
       zp_sand,                       & !I [frc] mass fraction of sand
       zp_dst_erod,                   & !I [frc] erodabilitC) de la surface (cover004 = 1 et cocer005=0.5)
       zp_wind10m,                    & !I [m/s] wind at 10m 
       zp_sfdst,                      & !O [kg/m2/sec] flux of dust 
       ksize                          & !I [nbr] number of points for which we do calculation
       )
  !
ELSE
  !
  CALL dustflux_get(                        &
           zp_ustar,                        &!I [m/s] wind friction speed over erodible surface
           zp_rhoa_2m,                      &!I [kg/m3] air density at surface
           zp_wg,                           &!I [m3/m3] volumetric water content of ground
           zp_z0_erod,                      &!I [m] roughness length of erodible surface
           zp_wsat,                         &!I [m3/m3] saturation water volumetric content
           zp_clay,                         &!I [frc] mass fraction of clay
           zp_sand,                         &!I [frc] mass fraction of sand
           zp_wind10m,                      &!I [m/s] wind at 10m 
           zp_sfdst,                        &!O [kg/m2/sec] flux of dust 
           ksize                            &!I [nbr] number of points for which we do calculation
           )
ENDIF  
    
zp_mss_frc_src(:,:) = 0.

IF (cemisparam_dst == "EXPLI" .OR. cemisparam_dst == "AMMA ") THEN

  ! For the repartition of mass fraction we uses the real USTAR
  zp_ustar(:) =  sqrt(zp_cd_dst(:))*zp_vmod(:)

  zseuil      = (/0., 0.35, 0.42, 0.50 ,0.66/)

  IF (cemisparam_dst == "EXPLI") THEN
    ! Compute modes mass fraction of emitted dust upon Alfaro et al. 1998
    zpcen(:,1)  = (/0., 0., 0.01, 0.08, 0.15, 0.15/)
    zpcen(:,2)  = (/0., 0., 0.36, 0.43, 0.76, 0.76/)
    ! Case of u* < 35E-2 m/s
    ! only coarse mode is activated 
    ! mode 1 = 0
    ! mode 2 = 0
    ! mode 3 = 100      
    ! Case of  35E-2 m/s < u* < 42E-2 m/s
    !  0 % < mode 1 < 1 %
    !  0 % < mode 2 < 36 %
    !  63 % < mode 3 < 100 %
    ! Case of  42E-2 m/s < u* < 50E-2 m/s
    !  1 % < mode 1 < 8 %
    !  36 % < mode 2 < 43 %
    !  49 % < mode 3 < 63 %
    ! Case of  50E-2 m/s < u* < 66E-2 m/s
    !  8 % < mode 1 < 15 %
    !  43 % < mode 2 < 76 %
    !  9  % < mode 3 < 49 %
    ! Case of  u* > 66E-2 m/s 
    !   mode 1 = 15 %
    !   mode 2 = 76 %
    !   mode 3 = 9 %
  ELSEIF (cemisparam_dst == "AMMA ") THEN
    ! For the repartition of mass fraction we uses the real USTAR upon
    ! Alfaro/Gomes 2001, jgr, table 5, soil type Salty sand.
    ! Percentage of mass considered between the fin and accumulation modes
    zpcen(:,1) = (/0., 0., 0.0023, 0.0185, 0.0345, 0.0345/)
    zpcen(:,2) = (/0., 0., 0.0077, 0.0615, 0.1155, 0.1155/)
    ! Case of u* < 35E-2 m/s
    ! only coarse mode is activated 
    ! mode 1 = 0
    ! mode 2 = 0
    ! mode 3 = 100 
    ! Case of  35E-2 m/s < u* < 42E-2 m/s
    !  0 % < mode 1 < 0.23 %
    !  0 % < mode 2 < 0.77 %
    !  99 % < mode 3 < 100 % 
    ! Case of  42E-2 m/s < u* < 50E-2 m/s
    !  0.23 % < mode 1 < 1.85 %
    !  0.77 % < mode 2 < 6.15 %
    !  92 % < mode 3 < 99 %    
    ! Case of  50E-2 m/s < u* < 66E-2 m/s
    !  1.85 % < mode 1 < 3.45 %
    !  6.15 % < mode 2 < 11.55 %
    !  85  % < mode 3 < 92 %
    ! Case of  u* > 66E-2 m/s 
    !   mode 1 = 15 %
    !   mode 2 = 76 %
    !   mode 3 = 9 %
  ENDIF

  DO js = 1, SIZE(zseuil) - 1
    WHERE (zp_ustar(:) >= zseuil(js) .AND. zp_ustar(:) < zseuil(js+1))
      zp_inter(:)         = (zp_ustar(:) - zseuil(js)) / (zseuil(js+1) - zseuil(js))
      zp_mss_frc_src(:,1) = zpcen(js,1) + (zpcen(js+1,1) - zpcen(js,1)) * zp_inter(:)
      zp_mss_frc_src(:,2) = zpcen(js,2) + (zpcen(js+1,2) - zpcen(js,2)) * zp_inter(:) 
    END WHERE
  ENDDO

  WHERE (zp_ustar(:) >= zseuil(SIZE(zseuil)))
    zp_mss_frc_src(:,1) = zpcen(SIZE(zseuil)+1,1) 
    zp_mss_frc_src(:,2) = zpcen(SIZE(zseuil)+1,2)
  END WHERE

  zp_mss_frc_src(:,3) = 1. - zp_mss_frc_src(:,1) - zp_mss_frc_src(:,2)

ELSE 
  DO jmode = 1,ndstmde
    zp_mss_frc_src(:,jorder_dst(jmode)) = dstk%XMSS_FRC_SRC(jmode)
  ENDDO
END IF

DO jmode=1,ndstmde
  ! Explicit mass fraction from u*
  zp_sfdst_mde(:,jmode) = zp_sfdst(:)                        &  !Total mass flux of dust
                          * zp_mss_frc_src(:,jorder_dst(jmode)) !Fraction of dust going to each mode
ENDDO


!Transfer the fluxes to each tile
DO jmode=1,ndstmde
  DO jj=1,ksize
    zsfdst_tile(kmask(jj),jveg,jmode) = zp_sfdst_mde(jj,jmode)
  ENDDO
ENDDO

IF (lhook) CALL dr_hook('COUPLING_DST_n:TREAT_SURF',1,zhook_handle)

END SUBROUTINE treat_surf
   
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE avg_flux_dst(   &
    KTILE                  & ! Number of different dust emitter vegetations
           )  

! Purpose: Average all fluxes from independent dust emitter surfaces (KTILE)
! Remember: XP_VEGTYPE_PATCH is m^2_{emittersurface}/m^{nature}
! The goal is to obtain PSFDST in kg_{dust}/m^2_{patch}

!INPUT
INTEGER, INTENT(IN)     :: KTILE                    ! Number of different dust emitter vegetation types

!LOCAL
INTEGER                 :: II                       ! Counter for points in patch-vector
INTEGER                 :: JJ                       ! Counter for vegetation types
INTEGER                 :: JMODE                    ! Counter for modes
INTEGER                 :: JSV_IDX                  ! Index for scalar variable
INTEGER                 :: NMOMENT                  ! Number of moments
REAL                    :: VEGFRAC_IN_PATCH         ! fraction of total vegetation in this patch (m^2_{patch}/m^2_{nature})
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!Remember: XP_VEGTYPE_PATCH is m^2_{emittersurface}/m^{nature}
!The goal is to obtain PSFDST in kg_{dust}/m^2_{patch}
!
!Start loop on modes
IF (lhook) CALL dr_hook('AVG_FLUX_DST',0,zhook_handle)
!
nmoment = int(SIZE(psfdst,2) / ndstmde)
!
DO jmode=1,ndstmde
  !Make sure dust mass emission is put in index 2, 5, 8 etc if 3 moments per mode
  !Make sure dust mass emission is put in index 2, 4, 6 etc if 2 moments per mode
  IF (nmoment == 1) THEN
    jsv_idx = (jmode-1)*nmoment + 1  !Counter for scalar variable
  ELSE
    jsv_idx = (jmode-1)*nmoment + 2  !Counter for scalar variable
  END IF
        
  !Start loop on number of dust emitter surfaces
  DO jj=1,ktile
    !Loop on points inside patch
    DO ii=1,ki
               
      !Get sum of vegetation fraction in this patch
      !fxm: VERY BAD LOOP ORDER
      vegfrac_in_patch = sum(pk%XVEGTYPE_PATCH(ii,:))
               
      !Get production of flux by adding up the contribution 
      !from the different tiles (here, "tiles" are dust emitter surfaces)
      psfdst(ii,jsv_idx) = psfdst(ii,jsv_idx)                  & ![kg/m^2_{patch}/sec] dust flux per patch 
                           + (zsfdst_tile(ii,jj,jmode)         & ![kg/m^2_{emittersurface}/sec] Dust flux per surface area of dust emitter surface
                           * pk%XVEGTYPE_PATCH(ii,dstk%NVT_DST(jj))  & ![frc] m^2_{emittersurface}/m^2_{nature}
                           / vegfrac_in_patch )                  ![frc] m^2_{patch}/m^2_{nature}  
    ENDDO !loop on point in patch
  ENDDO    !loop on different dust emitter surfaces
         
ENDDO !Loop on modes

IF (lhook) CALL dr_hook('AVG_FLUX_DST',1,zhook_handle)
      
END SUBROUTINE avg_flux_dst
    
END SUBROUTINE coupling_dst_n