phys/module_dust_emis.F

   1 !WRF:MEDIATION_LAYER:PHYSICS
   2 ! *** add new modules of schemes here
   3 !
   4 MODULE module_dust_emis
   5 CONTAINS
   6
   7 !====================================================================================
   8 ! simple dust emission scheme outside of WRF_CHEM
   9 !       revised from s/r mosaic_dust_gocartemis
  10 !         Mei Xu 2017Jan31 and additions by G. Thompson 2017Sep19
  11 !====================================================================================
  12   subroutine bulk_dust_emis (ktau,dt,num_soil_layers,u_phy,v_phy,          &
  13          rho_phy,alt,u10,v10,p8w,dz8w,smois,erod,                          &
  14          ivgtyp,isltyp,vegfra,albbck,xland,dx,g,                           &
  15          nifa2d,                                                           &
  16          ids,ide, jds,jde, kds,kde,                                        &
  17          ims,ime, jms,jme, kms,kme,                                        &
  18          its,ite, jts,jte, kts,kte                                         )
  19 !====================================================================================
  20 ! INPUT:
  21 ! ktau          - simulation time
  22 ! dt            - size of timestep
  23 ! u_phy,v_phy   - state variables of wind
  24 ! rho_phy,alt   - state variable of density and inverse of density
  25 ! dz8w,p8w      - model layer depth
  26 ! ivgtyp,isltyp - dominant vegetation, soil type
  27 ! vegfra        - vegetation fraction (0-100)
  28 ! albbck        - background sfc albedo used to increase erod where desert-like
  29 ! xland         - land mask (1: land, 2: water)
  30 ! dx, g         - model grid increment, gravity constant
  31 ! smois         - soil moisture
  32 ! erod          - fraction of erodible grid cell, for 1: Sand, 2: Silt, 3: Clay
  33 !               - read in from wrfinput
  34 ! OUTPUT:
  35 ! nifa2d        - dust generation in #/kg/s - to be used to update QNIFA (# kg-1)
  36 !------------------------------------------------------------------------------------
  37 ! Local variables:
  38 ! nmx           - number of dust size bins
  39 ! ilwi          - land/water flag
  40 ! bems          - binned dust emission in kg/timestep/cell
  41 ! totalemis     - dust emission from 0.1-10 um in radius, in ug/m2/s
  42 ! erodin        - fraction of erodible grid cel,
  43 ! w10m          - total wind speed at 10 m
  44 ! gwet          - volumetric soil moisture over porosity
  45 ! dxy           - model grid box area
  46 ! airmas        - grid box air mass
  47 ! airden        - air density
  48 ! DSRC          - source of each dust type           (kg/timestep/cell)
  49 !
  50 ! USED from module_data_gocart_dust:
  51 ! porosity      -
  52 ! ch_dust       - Constant to fudge the total emission of dust  (s2/m2)
  53 ! frac_s        - mass fraction of dust
  54 !====================================================================================
  55
  56   USE module_data_gocart_dust
  57
  58   IMPLICIT NONE
  59
  60    INTEGER,      INTENT(IN   ) :: ktau, num_soil_layers,           &
  61                                   ids,ide, jds,jde, kds,kde,               &
  62                                   ims,ime, jms,jme, kms,kme,               &
  63                                   its,ite, jts,jte, kts,kte
  64    INTEGER,DIMENSION( ims:ime , jms:jme )                  ,               &
  65           INTENT(IN   ) ::                                                 &
  66                                                      ivgtyp,               &
  67                                                      isltyp
  68    REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) ::                        &
  69                                                      vegfra,               &
  70                                                      albbck
  71    REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) ,      &
  72       INTENT(INOUT) ::                               smois
  73    REAL,  DIMENSION( ims:ime , jms:jme, 3 )                   ,               &
  74           INTENT(IN   ) ::    erod
  75    REAL,  DIMENSION( ims:ime , jms:jme )                   ,               &
  76           INTENT(IN   ) ::                                                 &
  77                                                      u10,                  &
  78                                                      v10,                  &
  79                                                      xland
  80    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  81           INTENT(IN   ) ::                                                 &
  82                                                         alt,               &
  83                                                      dz8w,p8w,             &
  84                                               u_phy,v_phy,rho_phy
  85
  86   REAL, INTENT(IN   ) :: dt,dx,g
  87   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: NIFA2D
  88 !
  89 ! local variables
  90 !
  91   integer :: i,j,k, nmx, month, n, m, maxflux
  92   real*8  w10m,gwet,den,diam,airden,airmas,erodin
  93   real*8  dxy
  94   real*8  converi
  95   REAL*8  u_ts0, u_ts, dsrc, srce
  96   real rhoa, g0, totalemis, pi, dustmas
  97
  98   converi=1.e9
  99   g0 = g*1.0E2
 100   pi = 3.14159
 101   dxy=dx*dx
 102   maxflux = 0
 103
 104   do j=jts,jte
 105   do i=its,ite
 106      nifa2d(i,j) = 0.0
 107   enddo
 108   enddo
 109
 110   nmx = 5
 111 ! consider 5 dust diameter/density as in module_data_gocart_dust
 112   do n = 1, nmx
 113
 114   den = den_dust(n)*1.0D-3            ! in g/cm3
 115   diam = 2.0*reff_dust(n)*1.0D2       ! in cm
 116   dustmas = (pi/6.) * den * (diam)**3 ! in g
 117   ch_dust(n,:)=1.0D-9  ! default is 1 ug s2 m-5 == 1.e-9 kg s2 m-5
 118   month = 3   ! it doesn't matter, ch_dust is not a month dependent now, a constant
 119   m = ipoint(n)  ! which of the 3 classes of fractional erod data is to be used; ipoint(3)=2
 120
 121   k=kts
 122   do j=jts,jte
 123   do i=its,ite
 124
 125     dsrc = 0.0
 126 !
 127 ! don't do dust over water!!!
 128 ! we still need to add in a check for snow cover on dusty ground. TO DO item
 129
 130     if(xland(i,j).lt.1.5) then
 131
 132       w10m=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
 133       airmas=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dxy/g   ! kg
 134
 135 ! we don't trust the u10,v10 values, if model layers are very thin near surface
 136       if(dz8w(i,kts,j).lt.12.)w10m=sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j))
 137
 138 ! consider using a gust factor of 33 percent higher than sustained wind (Cao and Fovell, 2018; also G. Bryan)
 139       w10m = w10m*1.33
 140
 141       erodin=erod(i,j,m)
 142
 143 ! increase erodability where the surface albedo is high to account better for real deserts
 144       if (erodin .gt. 1.E-8 .AND. albbck(i,j).gt.0.175 .and. vegfra(i,j).lt.12.5) then
 145          erodin = min(0.5d0, dble(erodin + 0.1*albbck(i,j)))
 146       endif
 147
 148 !  volumetric soil moisture over porosity
 149       gwet=smois(i,1,j)/porosity(isltyp(i,j))
 150       airden=rho_phy(i,kts,j)
 151       rhoa = airden*1.0D-3
 152
 153 ! Threshold velocity as a function of the dust density and the diameter from Bagnold (1941)
 154       u_ts0 = 0.13*1.0D-2*SQRT(den*g0*diam/rhoa)* &
 155               SQRT(1.0+0.006/den/g0/(diam)**2.5)/ &
 156               SQRT(1.928*(1331.0*(diam)**1.56+0.38)**0.092-1.0)
 157
 158 ! Case of surface dry enough to erode
 159       IF (gwet < 0.5) THEN  !  Pete's modified value
 160 !     IF (gwet < 0.2) THEN
 161          u_ts = max(0.0d0,dble(u_ts0*(1.2d0+2.0d-1*log10(max(1.0d-3, dble(gwet))))))
 162       ELSE
 163 ! Case of wet surface, no erosion
 164          u_ts = 100.0
 165       END IF
 166       srce = frac_s(n)*erodin*dxy  ! (m2)
 167 !     srce = 1.1*erodin*dxy  ! (m2)
 168       dsrc = max(0.0d0, dble(ch_dust(n,month)*srce*w10m*w10m *(w10m - u_ts)*dt)) ! (kg)
 169
 170 ! unit change from kg/timestep/cell to ug/m2/s
 171 !   totalemis=((dsrc)/dt)*converi/dxy
 172 !   totalemis=totalemis*1.01 !to account for the particles larger than 10 um based on assumed size distribution
 173
 174 ! convert dust flux to number concentration (#/kg/s)
 175 ! nifa2d = dsrc/dt / dustParticle_mas /cell_air_mass
 176
 177       nifa2d(i,j) = nifa2d(i,j) + dsrc/dt * 1000.0/dustmas/airmas
 178
 179     endif ! xland
 180     maxflux = MAX(maxflux,int(nifa2d(i,j)))
 181   enddo ! i
 182   enddo ! j
 183
 184   enddo ! n
 185 ! print*,'check maximum dust flux: ', maxflux
 186
 187 end subroutine bulk_dust_emis
 188
 189 END MODULE module_dust_emis