wrfv2_fire/chem/module_gocart_dust_afwa.F

   1 MODULE GOCART_DUST_AFWA
   2 !
   3 ! AFWA dust routine created by Sandra Jones (AER and AFWA) and Glenn Creighton (AFWA).
   4 !
   5
   6   USE module_data_gocart_dust
   7
   8   IMPLICIT NONE
   9
  10   INTRINSIC max, min
  11
  12 CONTAINS
  13   subroutine gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy,  &
  14          v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,dustin,snowh,zs,   &
  15          ivgtyp,isltyp,vegfra,xland,xlat,xlong,gsw,dx,g,emis_dust,         &
  16          ust,znt,clay,sand,alpha,gamma,smtune,                             &
  17          ids,ide, jds,jde, kds,kde,                                        &
  18          ims,ime, jms,jme, kms,kme,                                        &
  19          its,ite, jts,jte, kts,kte                                         )
  20   USE module_configure
  21   USE module_state_description
  22
  23    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  24
  25    INTEGER,      INTENT(IN   ) :: julday, ktau,                            &
  26                                   ids,ide, jds,jde, kds,kde,               &
  27                                   ims,ime, jms,jme, kms,kme,               &
  28                                   its,ite, jts,jte, kts,kte
  29    INTEGER,DIMENSION( ims:ime , jms:jme )                  ,               &
  30           INTENT(IN   ) ::                                                 &
  31                                                      ivgtyp,               &
  32                                                      isltyp
  33    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),                &
  34          INTENT(IN ) ::                              moist
  35    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
  36          INTENT(INOUT ) ::                           chem
  37    REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL,           &
  38          INTENT(INOUT ) ::                                                 &
  39          emis_dust
  40    REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) ,     &
  41          INTENT(IN   ) ::                            smois
  42    REAL, DIMENSION( config_flags%num_soil_layers ) ,                       &
  43          INTENT(IN   ) ::                            zs
  44    REAL, DIMENSION( ims:ime , jms:jme, ndcls )             ,               &
  45          INTENT(IN   ) ::                            erod
  46    REAL, DIMENSION( ims:ime , jms:jme, 5 )                 ,               &
  47          INTENT(INOUT) ::                            dustin
  48    REAL, DIMENSION( ims:ime , jms:jme )                    ,               &
  49          INTENT(IN   ) ::                                                  &
  50                                                      u10,                  &
  51                                                      v10,                  &
  52                                                      gsw,                  &
  53                                                      vegfra,               &
  54                                                      xland,                &
  55                                                      xlat,                 &
  56                                                      xlong,                &
  57                                                      ust,                  &
  58                                                      znt,                  &
  59                                                      clay,                 &
  60                                                      sand,                 &
  61                                                      snowh
  62    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                         &
  63          INTENT(IN   ) ::                                                  &
  64                                                      alt,                  &
  65                                                      t_phy,                &
  66                                                      dz8w,p8w,             &
  67                                                      u_phy,v_phy,rho_phy
  68   REAL, INTENT(IN   ) :: dt,dx,g
  69
  70 ! Local variables
  71
  72   INTEGER :: nmx,smx,i,j,k,imx,jmx,lmx,lhave
  73   INTEGER,DIMENSION (1,1) :: ilwi
  74   REAL*8, DIMENSION (1,1) :: erodtot
  75   REAL*8, DIMENSION (1,1) :: gravsm
  76   REAL*8, DIMENSION (1,1) :: drylimit
  77   REAL*8, DIMENSION (5)   :: tc,bems
  78   REAL*8, DIMENSION (1,1) :: airden,airmas,ustar
  79   REAL*8, DIMENSION (1) :: dxy
  80   REAL*8, DIMENSION (3) :: massfrac
  81   REAL :: conver,converi,volsm
  82   REAL*8 :: zwant
  83   REAL, INTENT(IN   ) :: alpha, gamma, smtune
  84
  85   conver=1.e-9
  86   converi=1.e9
  87
  88 ! Number of dust bins
  89
  90   imx=1
  91   jmx=1
  92   lmx=1
  93   nmx=ndust
  94   smx=ngsalt
  95
  96   k=kts
  97   DO j=jts,jte
  98   DO i=its,ite
  99
 100 ! Don't do dust over water!!!
 101
 102     IF (xland(i,j) .lt. 1.5) THEN
 103       ilwi(1,1)=1
 104
 105 ! Total concentration at lowest model level. This is still hardcoded for 5 bins.
 106
 107       tc(1)=chem(i,kts,j,p_dust_1)*conver
 108       tc(2)=chem(i,kts,j,p_dust_2)*conver
 109       tc(3)=chem(i,kts,j,p_dust_3)*conver
 110       tc(4)=chem(i,kts,j,p_dust_4)*conver
 111       tc(5)=chem(i,kts,j,p_dust_5)*conver
 112
 113 ! Air mass and density at lowest model level.
 114
 115       airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g
 116       airden(1,1)=rho_phy(i,kts,j)
 117       ustar(1,1)=ust(i,j)
 118       dxy(1)=dx*dx
 119
 120 ! Total erodibility.
 121
 122       erodtot(1,1)=SUM(erod(i,j,:))
 123
 124 ! Mass fractions of clay, silt, and sand.
 125
 126       massfrac(1)=clay(i,j)
 127       massfrac(2)=1-(clay(i,j)+sand(i,j))
 128       massfrac(3)=sand(i,j)
 129
 130 ! Don't allow roughness lengths greater than 20 cm to be lofted.
 131 ! This kludge accounts for land use types like urban areas and
 132 ! forests which would otherwise show up as high dust emitters.
 133 ! This is a placeholder for a more widely accepted kludge
 134 ! factor in the literature, which reduces lofting for rough areas.
 135 ! Forthcoming...
 136
 137       IF (znt(i,j) .gt. 0.2) then
 138         ilwi(1,1)=0
 139       ENDIF
 140
 141 ! Do not allow areas with bedrock, lava, or land-ice to loft.
 142
 143       IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. &
 144           isltyp(i,j) .eq. 18) then
 145         ilwi(1,1)=0
 146       ENDIF
 147
 148 ! Another hack to ensure dust does not loft from areas with snow
 149 ! cover...because, well, that doesn't make sense.
 150
 151       IF (snowh(i,j) .gt. 0.01) then
 152         ilwi(1,1)=0
 153       ENDIF
 154
 155 ! Check LSM scheme.  If RUC, add drypoint back to smois to get back to volumetric soil
 156 ! moisture.  If not RUC, smois is volumetric soil moisture. Volumetric soil moisture
 157 ! can be tuned here with a multiple (dust_smtune) set in the namelist.  Note: There
 158 ! can be a significant difference between the volumetric soil moisture calculated
 159 ! here between RUC and NOAA, even from the same input soil moisture.  This is due to
 160 ! fact that RUC subtracts hygroscopic water (drypoint) and rounds any value less
 161 ! than 0.005 to 0.005.  When we add drypoint back to the RUC soil moisture, we are
 162 ! artificially moistening any areas with initial values less than 0.005.
 163
 164       sfc_select: SELECT CASE(config_flags%sf_surface_physics)
 165         CASE (RUCLSMSCHEME)
 166           volsm=max((smois(i,1,j)+drypoint(isltyp(i,j)))*smtune,0.)
 167         CASE DEFAULT
 168           volsm=max(smois(i,1,j)*smtune,0.)
 169       END SELECT sfc_select
 170
 171 ! Calculate gravimetric soil moisture and drylimit.
 172
 173       gravsm(1,1)=100*volsm/((1.-porosity(isltyp(i,j)))*(2.65*(1-clay(i,j))+2.50*clay(i,j)))
 174       drylimit(1,1)=14.0*clay(i,j)*clay(i,j)+17.0*clay(i,j)
 175
 176 ! Call dust emission routine.
 177
 178       call source_dust(imx, jmx, lmx, nmx, smx, dt, tc, ustar, massfrac, &
 179                        erodtot, ilwi, dxy, gravsm, airden, airmas, &
 180                        bems, g, drylimit, alpha, gamma)
 181
 182       IF(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
 183         dustin(i,j,1:5)=tc(1:5)*converi
 184       ELSE
 185         chem(i,kts,j,p_dust_1)=tc(1)*converi
 186         chem(i,kts,j,p_dust_2)=tc(2)*converi
 187         chem(i,kts,j,p_dust_3)=tc(3)*converi
 188         chem(i,kts,j,p_dust_4)=tc(4)*converi
 189         chem(i,kts,j,p_dust_5)=tc(5)*converi
 190       ENDIF
 191
 192 ! For output diagnostics
 193
 194       emis_dust(i,1,j,p_edust1)=bems(1)
 195       emis_dust(i,1,j,p_edust2)=bems(2)
 196       emis_dust(i,1,j,p_edust3)=bems(3)
 197       emis_dust(i,1,j,p_edust4)=bems(4)
 198       emis_dust(i,1,j,p_edust5)=bems(5)
 199
 200     ENDIF
 201   ENDDO
 202   ENDDO
 203
 204 end subroutine gocart_dust_afwa_driver
 205
 206
 207   SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,&
 208                          erod, ilwi, dxy, gravsm, airden, airmas, &
 209                          bems, g0, drylimit, alpha, gamma)
 210
 211 ! ****************************************************************************
 212 ! *  Evaluate the source of each dust particles size bin by soil emission
 213 ! *
 214 ! *  Input:
 215 ! *         EROD      Fraction of erodible grid cell                (-)
 216 ! *         ILWI      Land/water flag                               (-)
 217 ! *         GRAVSM    Gravimetric soil moisture                     (g/g)
 218 ! *         DRYLIMIT  Upper GRAVSM limit for air-dry soil           (g/g)
 219 ! *         ALPHA     Constant to fudge the total emission of dust  (1/m)
 220 ! *         GAMMA     Exponential tuning constant for erodibility   (-)
 221 ! *         DXY       Surface of each grid cell                     (m2)
 222 ! *         AIRMAS    Mass of air for each grid box                 (kg)
 223 ! *         AIRDEN    Density of air for each grid box              (kg/m3)
 224 ! *         USTAR     Friction velocity                             (m/s)
 225 ! *         MASSFRAC  Fraction of mass in each of 3 soil classes    (-)
 226 ! *         DT1       Time step                                     (s)
 227 ! *         NMX       Number of dust bins                           (-)
 228 ! *         SMX       Number of saltation bins                      (-)
 229 ! *         IMX       Number of I points                            (-)
 230 ! *         JMX       Number of J points                            (-)
 231 ! *         LMX       Number of L points                            (-)
 232 ! *
 233 ! *  Data (see module_data_gocart_dust):
 234 ! *         SPOINT    Pointer to 3 soil classes                     (-)
 235 ! *         DEN_DUST  Dust density                                  (kg/m3)
 236 ! *         DEN_SALT  Saltation particle density                    (kg/m3)
 237 ! *         REFF_SALT Reference saltation particle diameter         (m)
 238 ! *         REFF_DUST Reference dust particle diameter              (m)
 239 ! *         LO_DUST   Lower diameter limits for dust bins           (m)
 240 ! *         UP_DUST   Upper diameter limits for dust bins           (m)
 241 ! *         FRAC_SALT Soil class mass fraction for saltation bins   (-)
 242 ! *
 243 ! *  Parameters:
 244 ! *         BETAMAX   Maximum sandblasting mass efficiency          (-)
 245 ! *         CMB       Constant of proportionality                   (-)
 246 ! *         MMD_DUST  Mass median diameter of dust                  (m)
 247 ! *         GSD_DUST  Geometric standard deviation of dust          (-)
 248 ! *         LAMBDA    Side crack propogation length                 (m)
 249 ! *         CV        Normalization constant                        (-)
 250 ! *         G0        Gravitational acceleration                    (m/s2)
 251 ! *         G         Gravitational acceleration in cgs             (cm/s2)
 252 ! *
 253 ! *  Working:
 254 ! *         BETA      Sandblasting mass efficiency                  (-)
 255 ! *         U_TS0     "Dry" threshold friction velocity             (m/s)
 256 ! *         U_TS      Moisture-adjusted threshold friction velocity (m/s)
 257 ! *         RHOA      Density of air in cgs                         (g/cm3)
 258 ! *         DEN       Dust density in cgs                           (g/cm3)
 259 ! *         DIAM      Dust diameter in cgs                          (cm)
 260 ! *         DMASS     Saltation mass distribution                   (-)
 261 ! *         DSURFACE  Saltation surface area per unit mass          (m2/kg)
 262 ! *         DS_REL    Saltation surface area distribution           (-)
 263 ! *         SALT      Saltation flux                                (kg/m/s)
 264 ! *         DLNDP     Dust bin width                                (-)
 265 ! *         EMIT      Total vertical mass flux                      (kg/m2/s)
 266 ! *         EMIT_VOL  Total vertical volume flux                    (m/s)
 267 ! *         DSRC      Mass of emitted dust               (kg/timestep/cell)
 268 ! *
 269 ! *  Output:
 270 ! *         TC        Total concentration of dust        (kg/kg/timestep/cell)
 271 ! *         BEMS      Source of each dust type           (kg/timestep/cell)
 272 ! *         USTART    Threshold friction vel. (bin 7)               (m/s)
 273 ! *
 274 ! ****************************************************************************
 275
 276   INTEGER, INTENT(IN)   :: nmx,imx,jmx,lmx,smx
 277   INTEGER, INTENT(IN)   :: ilwi(imx,jmx)
 278   REAL*8, INTENT(IN)    :: erod(imx,jmx)
 279   REAL*8, INTENT(IN)    :: ustar(imx,jmx)
 280   REAL*8, INTENT(IN)    :: gravsm(imx,jmx)
 281   REAL*8, INTENT(IN)    :: drylimit(imx,jmx)
 282   REAL*8, INTENT(IN)    :: dxy(jmx)
 283   REAL*8, INTENT(IN)    :: airden(imx,jmx,lmx), airmas(imx,jmx,lmx)
 284   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 285   REAL*8, INTENT(OUT)   :: bems(imx,jmx,nmx)
 286   REAL, INTENT(IN)    :: g0,dt1
 287
 288   REAL*8    :: den(smx), diam(smx)
 289   REAL*8    :: dvol(nmx), dlndp(nmx)
 290 ! REAL*8    :: distr_dust(nmx)
 291   REAL*8    :: dsurface(smx), ds_rel(smx)
 292   REAL*8    :: massfrac(3)
 293   REAL*8    :: u_ts0, u_ts, dsrc, srce, dmass, dvol_tot
 294   REAL*8    :: emit, emit_vol
 295   REAL      :: rhoa, g
 296   REAL*8    :: salt, stotal
 297   INTEGER   :: i, j, m, n, s
 298
 299 ! Global tuning constant, alpha.  Sandblasting mass efficiency, beta.
 300 ! Beta maxes out for clay fractions above 0.2 => betamax.
 301
 302   REAL, INTENT(IN)  :: alpha
 303   REAL, PARAMETER :: betamax=5.25E-4
 304   REAL*8 :: beta
 305
 306 ! Experimental optional exponential tuning constant for erodibility.
 307 ! 0 < gamma < 1 -> more relative impact by low erodibility regions.
 308 ! 1 < gamma -> accentuates high erodibility regions.  Recommend this
 309 ! be set to 1 (default) unless looking for a way to increase spread
 310 ! within an ensemble framework.
 311
 312   REAL, INTENT(IN) :: gamma
 313
 314 ! Constant of proportionality from Marticorena et al, 1997 (unitless)
 315 ! Arguably more ~consistent~ fudge than alpha, which has many walnuts
 316 ! sprinkled throughout the literature.
 317
 318   REAL, PARAMETER :: cmb=1.0    ! Marticorena et al,1997
 319 ! REAL, PARAMETER :: cmb=2.61   ! White,1979
 320
 321 ! Parameters used in Kok distribution function. Advise not to play with
 322 ! these without the expressed written consent of someone who knows what
 323 ! they're doing. (See Kok, 2010 PNAS for details on this scheme).
 324
 325 ! REAL, PARAMETER :: mmd_dust=3.4D-6  ! median mass diameter (m)
 326 ! REAL, PARAMETER :: gsd_dust=3.0     ! geometric standard deviation
 327 ! REAL, PARAMETER :: lambda=12.0D-6   ! crack propogation length (m)
 328 ! REAL, PARAMETER :: cv=12.62D-6      ! normalization constant
 329
 330 ! Calculate saltation surface area distribution from sand, silt, and clay
 331 ! mass fractions and saltation bin fraction. This will later become a
 332 ! modifier to the total saltation flux.  The reasoning here is that the
 333 ! size and availability of saltators affects saltation efficiency. Based
 334 ! on Eqn. (32) in Marticorena & Bergametti, 1995 (hereon, MB95).
 335
 336   DO n=1,smx
 337     dmass=massfrac(spoint(n))*frac_salt(n)
 338     dsurface(n)=0.75*dmass/(den_salt(n)*reff_salt(n))
 339   ENDDO
 340
 341 ! The following equation yields relative surface area fraction.  It will only
 342 ! work if you are representing the "full range" of all three soil classes.
 343 ! For this reason alone, we have incorporated particle sizes that encompass
 344 ! the clay class, to account for its relative area over the basal
 345 ! surface, even though these smaller bins would be unlikely to play any large
 346 ! role in the actual saltation process.
 347
 348   stotal=SUM(dsurface(:))
 349   DO n=1,smx
 350     ds_rel(n)=dsurface(n)/stotal
 351   ENDDO
 352
 353 ! Calculate total dust emission due to saltation of sand sized particles.
 354 ! Begin by calculating DRY threshold friction velocity (u_ts0).  Next adjust
 355 ! u_ts0 for moisture to get threshold friction velocity (u_ts). Then
 356 ! calculate saltation flux (salt) where ustar has exceeded u_ts.  Finally,
 357 ! calculate total dust emission (tot_emit), taking into account erodibility.
 358
 359  g = g0*1.0E2                          ! (cm s^-2)
 360  emit=0.0
 361
 362  DO n = 1, smx                         ! Loop over saltation bins
 363    den(n) = den_salt(n)*1.0D-3         ! (g cm^-3)
 364    diam(n) = 2.0*reff_salt(n)*1.0D2    ! (cm)
 365    DO i = 1,imx
 366      DO j = 1,jmx
 367        rhoa = airden(i,j,1)*1.0D-3     ! (g cm^-3)
 368
 369      ! Threshold friction velocity as a function of the dust density and
 370      ! diameter from Bagnold (1941) (m s^-1).
 371
 372        u_ts0 = 0.13*1.0D-2*SQRT(den(n)*g*diam(n)/rhoa)* &
 373                SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ &
 374                SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0)
 375
 376      ! Friction velocity threshold correction function based on physical
 377      ! properties related to moisture tension. Soil moisture greater than
 378      ! dry limit serves to increase threshold friction velocity (making
 379      ! it more difficult to loft dust). When soil moisture has not reached
 380      ! dry limit, treat as dry (no correction to threshold friction
 381      ! velocity). GC
 382
 383        IF (gravsm(i,j) > drylimit(i,j)) THEN
 384          u_ts = MAX(0.0D+0,u_ts0*(sqrt(1.0+1.21*(gravsm(i,j)-drylimit(i,j))**0.68)))
 385        ELSE
 386          u_ts = u_ts0
 387        END IF
 388
 389      ! Saltation flux from Marticorena & Bergametti 1995 (MB95). ds_rel is
 390      ! the relative surface area distribution
 391
 392        IF (ustar(i,j) .gt. u_ts .and. erod(i,j) .gt. 0.0 .and. ilwi(i,j) == 1) THEN
 393          salt = cmb*ds_rel(n)*(airden(i,j,1)/g0)*(ustar(i,j)**3)* &
 394                 (1. + u_ts/ustar(i,j))*(1. - (u_ts**2)/(ustar(i,j)**2))  ! (kg m^-1 s^-1)
 395        ELSE
 396          salt = 0.D0
 397        ENDIF
 398
 399      ! Calculate total vertical mass flux (note beta has units of m^-1)
 400      ! Beta acts to tone down dust in areas with so few dust-sized particles that the
 401      ! lofting efficiency decreases.  Otherwise, super sandy zones would be huge dust
 402      ! producers, which is generally not the case.  Equation derived from wind-tunnel
 403      ! experiments (see MB95).
 404
 405        beta=10**(13.6*massfrac(1)-6.0)  ! (unitless)
 406        if (beta .gt. betamax) then
 407          beta=betamax
 408        endif
 409        emit=emit+salt*(erod(i,j)**gamma)*alpha*beta    ! (kg m^-2 s^-1)
 410      END DO
 411    END DO
 412  END DO                                 ! End do over saltation bins
 413
 414 ! Now that we have the total dust emission, distribute into dust bins using
 415 ! lognormal distribution (Dr. Jasper Kok, 2010), and
 416 ! calculate total mass emitted over the grid box over the timestep.
 417 !
 418 ! In calculating the Kok distribution, we assume upper and lower limits to each bin.
 419 ! For reff_dust=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/) (default),
 420 ! lower limits were ASSUMED at lo_dust=(/0.1D-6,1.0D-6,1.8D-6,3.0D-6,6.0D-6/)
 421 ! upper limits were ASSUMED at up_dust=(/1.0D-6,1.8D-6,3.0D-6,6.0D-6,10.0D-6/)
 422 ! These may be changed within module_data_gocart_dust.F, but make sure it is
 423 ! consistent with reff_dust values.  These values were taken from the original
 424 ! GOCART bin configuration. We use them here to calculate dust bin width, dlndp.
 425 ! dvol is the volume distribution. You know...if you were wondering. GC
 426
 427 ! dvol_tot=0.
 428 ! DO n=1,nmx  ! Loop over all dust bins
 429 !   dlndp(n)=LOG(up_dust(n)/lo_dust(n))
 430 !   dvol(n)=(2.0*reff_dust(n)/cv)*(1.+ERF(LOG(2.0*reff_dust(n)/mmd_dust)/(SQRT(2.)*LOG(gsd_dust))))*&
 431 !         EXP(-(2.0*reff_dust(n)/lambda)**3.0)*dlndp(n)
 432 !   dvol_tot=dvol_tot+dvol(n)
 433 !
 434 !  ! Convert mass flux to volume flux
 435 !   emit_vol=emit/den_dust(n) ! (m s^-1)
 436 ! END DO
 437 ! DO n=1,nmx  ! Loop over all dust bins
 438 !   distr_dust(n)=dvol(n)/dvol_tot
 439 ! END DO
 440
 441 ! Now distribute total vertical emission into dust bins and update concentration.
 442
 443  DO n=1,nmx  ! Loop over all dust bins
 444    DO i=1,imx
 445      DO j=1,jmx
 446       ! Calculate total mass emitted
 447         dsrc = emit*den_dust(n)*distr_dust(n)*dxy(j)*dt1  ! (kg)
 448         IF (dsrc < 0.0) dsrc = 0.0
 449
 450       ! Update dust mixing ratio at first model level.
 451         tc(i,j,1,n) = tc(i,j,1,n) + dsrc / airmas(i,j,1) ! (kg/kg)
 452         !bems(i,j,n) = dsrc  ! diagnostic
 453         bems(i,j,n) = 1000.*dsrc/(dxy(j)*dt1) ! diagnostic (g/m2/s)
 454      END DO
 455    END DO
 456  END DO
 457
 458 END SUBROUTINE source_dust
 459
 460
 461 END MODULE GOCART_DUST_AFWA