chem/module_gocart_dust_afwa.F

   1 MODULE GOCART_DUST_AFWA
   2 !
   3 ! AFWA dust routine
   4 ! Created by Sandra Jones (AER and AFWA) and Glenn Creighton (AFWA).
   5 !
   6 ! A. Ukhov, 11 March 2021, Now "emis_dust" is accumulated dust
   7 ! emission (kg/m2). Before was instantenious flux (g m^-2 s^-1).
   8 ! Code cleanup, remove unused variables.
   9
  10   USE module_data_gocart_dust
  11
  12   IMPLICIT NONE
  13
  14   INTRINSIC max, min
  15
  16 CONTAINS
  17   SUBROUTINE gocart_dust_afwa_driver(dt,config_flags,alt,  &
  18          chem,rho_phy,smois,u10,v10,dz8w,erod,erod_dri,dustin,snowh,   &
  19          isltyp,vegfra,lai_vegmask,xland,g,emis_dust,      &
  20          ust,znt,clay_wrf,sand_wrf,clay_nga,sand_nga,afwa_dustloft,                 &
  21          tot_dust,tot_edust,vis_dust,alpha,gamma,smtune,ustune,            &
  22          ids,ide, jds,jde, kds,kde,                                        &
  23          ims,ime, jms,jme, kms,kme,                                        &
  24          its,ite, jts,jte, kts,kte                                         )
  25   USE module_configure
  26   USE module_state_description
  27   USE module_data_sorgam, ONLY: factnuma,factnumc,soilfac
  28
  29    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  30
  31    INTEGER,      INTENT(IN   ) :: ids,ide, jds,jde, kds,kde,               &
  32                                   ims,ime, jms,jme, kms,kme,               &
  33                                   its,ite, jts,jte, kts,kte
  34    INTEGER,DIMENSION( ims:ime , jms:jme )                  ,               &
  35           INTENT(IN   ) ::    isltyp
  36
  37    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
  38          INTENT(INOUT ) ::                           chem
  39    REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL,           &
  40          INTENT(INOUT ) ::    emis_dust
  41
  42    REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) ,     &
  43          INTENT(IN   ) ::                            smois
  44    REAL, DIMENSION( ims:ime , jms:jme, ndcls )             ,               &
  45          INTENT(IN   ) ::                            erod,erod_dri
  46    REAL, DIMENSION( ims:ime , jms:jme, 5 )                 ,               &
  47          INTENT(INOUT) ::                            dustin
  48    REAL, DIMENSION( ims:ime , jms:jme )                    ,               &
  49          INTENT(IN   ) ::                            u10,                  &
  50                                                      v10,                  &
  51                                                      vegfra,               &
  52                                                      lai_vegmask,          &
  53                                                      xland,                &
  54                                                      ust,                  &
  55                                                      znt,                  &
  56                                                      clay_wrf,             &
  57                                                      sand_wrf,             &
  58                                                      clay_nga,             &
  59                                                      sand_nga,             &
  60                                                      snowh
  61    REAL, DIMENSION( ims:ime , kms:kme , jms:jme )   ,                      &
  62          INTENT(IN   ) ::                            alt,                  &
  63                                                      dz8w,                 &
  64                                                      rho_phy
  65   REAL, DIMENSION( ims:ime , jms:jme )               ,                     &
  66          INTENT(  OUT) ::                            afwa_dustloft,        &
  67                                                      tot_edust
  68   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                          &
  69          INTENT(  OUT) ::                            tot_dust,             &
  70                                                      vis_dust
  71   REAL, INTENT(IN   ) :: dt,g
  72
  73   ! Local variables
  74
  75   INTEGER :: nmx,smx,i,j,k,imx,jmx,lmx
  76   INTEGER,DIMENSION (1,1) :: ilwi
  77   REAL*8, DIMENSION (1,1) :: erodtot
  78   REAL*8, DIMENSION (1,1) :: vegmask
  79   REAL*8, DIMENSION (1,1) :: gravsm
  80   REAL, DIMENSION( ims:ime , jms:jme ) :: clay,sand
  81   REAL*8, DIMENSION (1,1) :: drylimit
  82   REAL*8, DIMENSION (5)   :: tc,bems
  83   REAL*8, DIMENSION (1,1) :: airden,ustar
  84   REAL*8, DIMENSION (3) :: massfrac
  85   REAL*8                :: volsm
  86   REAL :: conver,converi
  87   REAL :: psi,ustart,w10
  88   REAL*8 :: zwant
  89   REAL, INTENT(IN   ) :: alpha, gamma, smtune, ustune
  90   INTEGER :: smois_opt
  91   real    :: dz_lowest
  92
  93   conver=1.e-9
  94   converi=1.e9
  95
  96   ! Number of dust bins
  97
  98   imx=1
  99   jmx=1
 100   lmx=1
 101   nmx=ndust
 102   smx=ngsalt
 103
 104   k=kts
 105   DO j=jts,jte
 106   DO i=its,ite
 107
 108     ! Masked value for afwa_dustloft
 109
 110     afwa_dustloft(i,j)=-99.
 111
 112     ! Don't do dust over water!!!
 113
 114     IF (xland(i,j) .lt. 1.5) THEN
 115       ilwi(1,1)=1
 116
 117       ! Total concentration at lowest model level. This is still hardcoded
 118       ! for 5 bins.
 119
 120       IF(config_flags%chem_opt == CB05_SORG_AQ_KPP .or. &
 121          config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP ) then
 122         tc(1)=0.0
 123         tc(2)=0.0
 124         tc(3)=0.0
 125         tc(4)=0.0
 126         tc(5)=0.0
 127       ELSE
 128         tc(1)=chem(i,kts,j,p_dust_1)*conver
 129         tc(2)=chem(i,kts,j,p_dust_2)*conver
 130         tc(3)=chem(i,kts,j,p_dust_3)*conver
 131         tc(4)=chem(i,kts,j,p_dust_4)*conver
 132         tc(5)=chem(i,kts,j,p_dust_5)*conver
 133       END IF
 134
 135       ! Air mass and density at lowest model level.
 136
 137       airden(1,1)=rho_phy(i,kts,j)
 138       ustar(1,1)=ust(i,j)
 139       dz_lowest = dz8w(i,1,j)
 140
 141       ! Friction velocity tuning constant (Note: recommend 0.7 for PXLSM,
 142       ! else use 1.0.  This was created due to make the scheme compatible
 143       ! with the much stronger friction velocities coming out of PXLSM).
 144
 145       ustar(1,1)=ustar(1,1) * ustune
 146
 147       ! Total erodibility.  Determine what DSR we're using.
 148       ! Note, the DRI erodibility dataset is an optional 1km resolution dataset
 149       ! but currently is only available over southwest Asia. If running
 150       ! a domain outside of SWA, this will fill in the missing data from
 151       ! DRI with the Ginoux dataset, but there will be inconsistencies. GAC
 152       ! Erodibility is broken up into 3 bins, sum them up here.
 153
 154       IF (config_flags%dust_dsr .eq. 1) then  ! DRI DSR
 155           IF (erod_dri(i,j,1).ge.0) THEN  ! Where DRI is defined
 156              erodtot(1,1) = SUM(erod_dri(i,j,:))
 157           ELSE ! Outside where DRI not defined, use Ginoux
 158              erodtot(1,1)=SUM(erod(i,j,:))
 159           ENDIF
 160       ELSE ! Ginoux DSR
 161           erodtot(1,1)=SUM(erod(i,j,:))
 162       ENDIF
 163
 164       ! Set the vegmask variable to the desired vegation mask at this gridpoint
 165
 166       IF (config_flags%dust_veg .eq. 1) then
 167
 168          ! 12-month vegetation fraction
 169          ! If user chose this 12-month greenfrac vegetation mask option,
 170          ! cut off everything above 5%
 171
 172          IF (vegfra(i,j) .ge. 5.) then
 173              vegmask(1,1) = 0.0
 174          ELSE
 175              vegmask(1,1) = 1.0
 176          ENDIF
 177       ELSE IF (config_flags%dust_veg .eq. 2) then
 178
 179          ! 8-day MODIS LAI vegmask
 180          ! 1 = no veg, produce dust; 0 = vegation, do not produce dust
 181
 182          vegmask(1,1) = lai_vegmask(i,j)
 183       ELSE
 184
 185          ! Default choice = static ginoux vegmask
 186
 187          IF (erod(i,j,1) .eq. 0) THEN
 188             vegmask(1,1) = 0.0
 189          ELSE
 190             vegmask(1,1) = 1.0
 191          ENDIF
 192       ENDIF
 193
 194       ! Remove vegetated areas (vegmask=0) from total erodibility.
 195
 196       erodtot(1,1) = erodtot(1,1) * vegmask(1,1)
 197
 198       ! Option to use an optional high resolution soil type database from NGA.
 199       ! Option 0 = WRF (default); Option 1 = NGA
 200       ! Note NGA dataset is currently only available over southwest Asia
 201       ! Until a global dataset becomes available, option 0 is recommended
 202       ! for consistency.  If option 1 is chosen for a domain outside of
 203       ! SWA, this logic will fill in the areas missing from NGA with the
 204       ! defaults from WRF. It will work, but it will be inconsistent. GAC
 205
 206       IF (config_flags%dust_soils .eq. 1) then
 207          IF (clay_nga(i,j) .ge.0) then
 208              clay(i,j) = clay_nga(i,j)
 209              sand(i,j) = sand_nga(i,j)
 210          ELSE
 211              clay(i,j) =clay_wrf(i,j)
 212              sand(i,j) =sand_wrf(i,j)
 213          ENDIF
 214       ELSE
 215          clay(i,j) =clay_wrf(i,j)
 216          sand(i,j) =sand_wrf(i,j)
 217       ENDIF
 218
 219       ! Mass fractions of clay, silt, and sand.
 220
 221       massfrac(1)=clay(i,j)
 222       massfrac(2)=1-(clay(i,j)+sand(i,j))
 223       massfrac(3)=sand(i,j)
 224
 225       ! Don't allow roughness lengths greater than 20 cm to be lofted.
 226       ! This kludge accounts for land use types like urban areas and
 227       ! forests which would otherwise show up as high dust emitters.
 228       ! This is a placeholder for a more widely accepted kludge
 229       ! factor in the literature, which reduces lofting for rough areas.
 230       ! Forthcoming...
 231
 232       IF (znt(i,j) .gt. 0.2) then
 233         ilwi(1,1)=0
 234       ENDIF
 235
 236       ! Do not allow areas with bedrock, lava, or land-ice to loft.
 237
 238       IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. &
 239           isltyp(i,j) .eq. 18) then
 240         ilwi(1,1)=0
 241       ENDIF
 242
 243       ! Another hack to ensure dust does not loft from areas with snow
 244       ! cover...because, well, that doesn't make sense.
 245
 246       IF (snowh(i,j) .gt. 0.01) then
 247         ilwi(1,1)=0
 248       ENDIF
 249
 250       ! Volumetric soil moisture can be tuned here with a dust_smtune
 251       ! set in the namelist.
 252
 253       volsm=max(smois(i,1,j)*smtune,0.)
 254
 255       ! Calculate gravimetric soil moisture.
 256
 257       gravsm(1,1)=100*volsm/((1.-porosity(isltyp(i,j)))*(2.65*(1-clay(i,j))+2.50*clay(i,j)))
 258
 259       ! Choose an LSM option and drylimit option.
 260       ! Drylimit calculations based on look-up table in
 261       ! Clapp and Hornberger (1978) for RUC and PXLSM and
 262       ! Cosby et al. (1984) for Noah LSM.
 263
 264       smois_opt = 0
 265       IF (config_flags%dust_smois == 1) then
 266          sfc_select: SELECT CASE(config_flags%sf_surface_physics)
 267             CASE ( RUCLSMSCHEME, PXLSMSCHEME )
 268                drylimit(1,1) =0.035*(13.52*clay(i,j)+3.53)**2.68
 269                smois_opt = 1
 270             CASE ( LSMSCHEME )
 271                drylimit(1,1) =0.0756*(15.127*clay(i,j)+3.09)**2.3211
 272                smois_opt = 1
 273             CASE DEFAULT
 274
 275                ! Don't currently support volumetric soil moisture
 276                ! for this scheme, use drylimit based on gravimetric
 277
 278                drylimit(1,1)=14.0*clay(i,j)*clay(i,j)+17.0*clay(i,j)
 279          END SELECT sfc_select
 280       ELSE
 281
 282          ! use drylimit based on gravimetric soil moisture
 283
 284          drylimit(1,1)=14.0*clay(i,j)*clay(i,j)+17.0*clay(i,j)
 285       END IF
 286
 287       ! Call dust emission routine.
 288
 289       call source_dust(imx, jmx, lmx, nmx, smx, dt, tc, ustar, massfrac, &
 290                        erodtot, ilwi, gravsm, volsm, airden,dz_lowest,   &
 291                        bems, ustart, g, drylimit, alpha, gamma, smois_opt)
 292
 293       IF(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
 294         dustin(i,j,1:5)=tc(1:5)*converi
 295       ELSE IF(config_flags%chem_opt == CB05_SORG_AQ_KPP .or.  &
 296               config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP) then
 297 !KW            chem(i,kts,j,p_p25i)=chem(i,kts,j,p_p25i) &
 298 !                         +.25*(tc(1)+.286*tc(2))*converi
 299 !            chem(i,kts,j,p_p25i)=max(chem(i,kts,j,p_p25i),1.e-16)
 300 !            chem(i,kts,j,p_p25j)=chem(i,kts,j,p_p25j) &
 301 !                         +.75*(tc(1)+.286*tc(2))*converi
 302 !            chem(i,kts,j,p_p25j)=max(chem(i,kts,j,p_p25j),1.e-16)
 303 !            chem(i,kts,j,p_soila)=chem(i,kts,j,p_soila) &
 304 !                         +(.714*tc(2)+tc(3))*converi
 305 !            chem(i,kts,j,p_soila)=max(chem(i,kts,j,p_soila),1.e-16)
 306          chem(i,kts,j,p_p25j)=chem(i,kts,j,p_p25j) + 0.03*sum(tc(1:5))*converi
 307          chem(i,kts,j,p_soila)=chem(i,kts,j,p_soila) + 0.97*1.02*sum(tc(1:5))*converi
 308
 309          chem(i,kts,j,p_ac0) =  chem(i,kts,j,p_ac0) +  0.03*sum(tc(1:5))*converi*factnuma*soilfac
 310          chem(i,kts,j,p_corn) =  chem(i,kts,j,p_corn) + 0.97*1.02*sum(tc(1:5))*converi*factnumc*soilfac
 311       ELSE
 312         chem(i,kts,j,p_dust_1)=tc(1)*converi ! tc(1...5) is (kg/kg), p_dust_1...5 (ug/kg)
 313         chem(i,kts,j,p_dust_2)=tc(2)*converi
 314         chem(i,kts,j,p_dust_3)=tc(3)*converi
 315         chem(i,kts,j,p_dust_4)=tc(4)*converi
 316         chem(i,kts,j,p_dust_5)=tc(5)*converi
 317       ENDIF
 318
 319       ! Diagnostic dust lofting potential diagnostic calculation
 320
 321       psi=0.
 322       w10=(u10(i,j)**2.+v10(i,j)**2.)**0.5
 323       IF (ustar(1,1) .ne. 0. .and. znt(i,j) .ne. 0.) THEN
 324         psi=0.4*w10/ustar(1,1)-LOG(10.0/znt(i,j))
 325       ENDIF
 326       IF (erodtot(1,1) .gt. 0.) then
 327         afwa_dustloft(i,j)=ustune*w10-ustart*(LOG(10.0/znt(i,j))+psi)/0.4
 328       ENDIF
 329
 330      ! A. Ukhov
 331      ! for output diagnostics
 332      ! bems (kg/m2) per dt
 333      ! p_edust1...5 is accumulated dust emission (kg/m2)
 334      emis_dust(i,1,j,p_edust1)=emis_dust(i,1,j,p_edust1)+bems(1)
 335      emis_dust(i,1,j,p_edust2)=emis_dust(i,1,j,p_edust2)+bems(2)
 336      emis_dust(i,1,j,p_edust3)=emis_dust(i,1,j,p_edust3)+bems(3)
 337      emis_dust(i,1,j,p_edust4)=emis_dust(i,1,j,p_edust4)+bems(4)
 338      emis_dust(i,1,j,p_edust5)=emis_dust(i,1,j,p_edust5)+bems(5)
 339
 340
 341      ! Diagnostic accumulated total emitted dust (kg/m2)
 342      tot_edust(i,j)=tot_edust(i,j)+(bems(1)+bems(2)+bems(3)+bems(4)+bems(5))
 343
 344     ENDIF
 345
 346     ! Cumulative dust concentration (ug m^-3) and visibility (m) diagnostics
 347     ! Note visibility is capped at 20 km.  Simple visibility algorithm based
 348     ! on mean particle diameter for each dust bin - perfect spheres - perfect
 349     ! blackbodies.
 350
 351     DO k=kts,kte
 352       tot_dust(i,k,j)=(chem(i,k,j,p_dust_1)+chem(i,k,j,p_dust_2)+       &
 353                        chem(i,k,j,p_dust_3)+chem(i,k,j,p_dust_4)+       &
 354                        chem(i,k,j,p_dust_5))*rho_phy(i,k,j)
 355       IF ( tot_dust(i,k,j) .gt. 0. ) THEN
 356         vis_dust(i,k,j)=MIN(3.912/                                      &
 357                           ((1.470E-6*chem(i,k,j,p_dust_1)+              &
 358                             7.877E-7*chem(i,k,j,p_dust_2)+              &
 359                             4.623E-7*chem(i,k,j,p_dust_3)+              &
 360                             2.429E-7*chem(i,k,j,p_dust_4)+              &
 361                             1.387E-7*chem(i,k,j,p_dust_5))*             &
 362                            rho_phy(i,k,j)),999999.)
 363       ELSE
 364         vis_dust(i,k,j)=999999.
 365       ENDIF
 366     ENDDO
 367
 368   ENDDO
 369   ENDDO
 370
 371   END SUBROUTINE gocart_dust_afwa_driver
 372
 373   SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,&
 374                          erod, ilwi, gravsm, volsm, airden, dz_lowest,     &
 375                          bems, ustart, g0, drylimit, alpha, gamma, smois_opt)
 376
 377   ! ****************************************************************************
 378   ! *  Evaluate the source of each dust particles size bin by soil emission
 379   ! *
 380   ! *  Input:
 381   ! *         EROD      Fraction of erodible grid cell                (-)
 382   ! *         ILWI      Land/water flag                               (-)
 383   ! *         GRAVSM    Gravimetric soil moisture                     (g/g)
 384   ! *         VOLSM     Volumetric soil moisture                      (g/g)
 385   ! *         SOILM_OPT Soil moisture option (1:Use GRAVSM 2:VOLSM)   (-)
 386   ! *         DRYLIMIT  Upper GRAVSM (VOLSM) limit for air-dry soil   (g/g)
 387   ! *         ALPHA     Constant to fudge the total emission of dust  (1/m)
 388   ! *         GAMMA     Exponential tuning constant for erodibility   (-)
 389   ! *         AIRDEN    Density of air for each grid box              (kg/m3)
 390   ! *         USTAR     Friction velocity                             (m/s)
 391   ! *         MASSFRAC  Fraction of mass in each of 3 soil classes    (-)
 392   ! *         DT1       Time step                                     (s)
 393   ! *         NMX       Number of dust bins                           (-)
 394   ! *         SMX       Number of saltation bins                      (-)
 395   ! *         IMX       Number of I points                            (-)
 396   ! *         JMX       Number of J points                            (-)
 397   ! *         LMX       Number of L points                            (-)
 398   ! *         dz_lowest heigth of the lowest layer                    (m)
 399   ! *
 400   ! *  Data (see module_data_gocart_dust):
 401   ! *         SPOINT    Pointer to 3 soil classes                     (-)
 402   ! *         DEN_DUST  Dust density                                  (kg/m3)
 403   ! *         DEN_SALT  Saltation particle density                    (kg/m3)
 404   ! *         REFF_SALT Reference saltation particle diameter         (m)
 405   ! *         REFF_DUST Reference dust particle diameter              (m)
 406   ! *         LO_DUST   Lower diameter limits for dust bins           (m)
 407   ! *         UP_DUST   Upper diameter limits for dust bins           (m)
 408   ! *         FRAC_SALT Soil class mass fraction for saltation bins   (-)
 409   ! *
 410   ! *  Parameters:
 411   ! *         BETAMAX   Maximum sandblasting mass efficiency          (-)
 412   ! *         CMB       Constant of proportionality                   (-)
 413   ! *         MMD_DUST  Mass median diameter of dust                  (m)
 414   ! *         GSD_DUST  Geometric standard deviation of dust          (-)
 415   ! *         LAMBDA    Side crack propogation length                 (m)
 416   ! *         CV        Normalization constant                        (-)
 417   ! *         G0        Gravitational acceleration                    (m/s2)
 418   ! *         G         Gravitational acceleration in cgs             (cm/s2)
 419   ! *
 420   ! *  Working:
 421   ! *         BETA      Sandblasting mass efficiency                  (-)
 422   ! *         U_TS0     "Dry" threshold friction velocity             (m/s)
 423   ! *         U_TS      Moisture-adjusted threshold friction velocity (m/s)
 424   ! *         RHOA      Density of air in cgs                         (g/cm3)
 425   ! *         DEN       Dust density in cgs                           (g/cm3)
 426   ! *         DIAM      Dust diameter in cgs                          (cm)
 427   ! *         DMASS     Saltation mass distribution                   (-)
 428   ! *         DSURFACE  Saltation surface area per unit mass          (m2/kg)
 429   ! *         DS_REL    Saltation surface area distribution           (-)
 430   ! *         SALT      Saltation flux                                (kg/m/s)
 431   ! *         DLNDP     Dust bin width                                (-)
 432   ! *         EMIT      Total vertical mass flux                      (kg/m2/s)
 433   ! *         EMIT_VOL  Total vertical volume flux                    (m/s)
 434   ! *         DSRC      Mass of emitted dust               (kg/timestep/m2)
 435   ! *
 436   ! *  Output:
 437   ! *         TC        Total concentration of dust                 (kg/kg)
 438   ! *         BEMS      Source of each dust type           (kg/timestep/m2)
 439   ! *         USTART    Threshold friction vel. (bin 7)               (m/s)
 440   ! *
 441   ! ****************************************************************************
 442
 443   INTEGER, INTENT(IN)   :: nmx,imx,jmx,lmx,smx
 444   INTEGER, INTENT(IN)   :: ilwi(imx,jmx)
 445   REAL*8, INTENT(IN)    :: erod(imx,jmx)
 446   REAL*8, INTENT(IN)    :: ustar(imx,jmx)
 447   REAL*8, INTENT(IN)    :: gravsm(imx,jmx)
 448   REAL*8, INTENT(IN)    :: drylimit(imx,jmx)
 449   REAL*8, INTENT(IN)    :: airden(imx,jmx,lmx)
 450   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 451   REAL*8, INTENT(OUT)   :: bems(imx,jmx,nmx)
 452   REAL, INTENT(IN)    :: g0,dt1
 453   REAL, INTENT(OUT)   :: ustart
 454   INTEGER, INTENT(IN) :: smois_opt
 455   REAL*8, INTENT(IN)    :: volsm
 456   REAL,   INTENT(IN)  :: dz_lowest
 457
 458   REAL*8    :: den(smx), diam(smx)
 459 ! REAL*8    :: dvol(nmx), dlndp(nmx)
 460 ! REAL*8    :: distr_dust(nmx)
 461   REAL*8    :: dsurface(smx), ds_rel(smx)
 462   REAL*8    :: massfrac(3)
 463   REAL*8    :: u_ts0, u_ts, dsrc, dmass!, dvol_tot
 464   REAL*8    :: emit!, emit_vol
 465   REAL      :: rhoa, g
 466   REAL*8    :: salt, stotal
 467   INTEGER   :: i, j, m, n, s
 468
 469 ! Global tuning constant, alpha.  Sandblasting mass efficiency, beta.
 470 ! Beta maxes out for clay fractions above 0.2 => betamax.
 471
 472   REAL, INTENT(IN)  :: alpha
 473   REAL, PARAMETER :: betamax=5.25E-4
 474   REAL*8 :: beta
 475
 476 ! Experimental optional exponential tuning constant for erodibility.
 477 ! 0 < gamma < 1 -> more relative impact by low erodibility regions.
 478 ! 1 < gamma -> accentuates high erodibility regions.  Recommend this
 479 ! be set to 1 (default) unless looking for a way to increase spread
 480 ! within an ensemble framework.
 481
 482   REAL, INTENT(IN) :: gamma
 483
 484 ! Constant of proportionality from Marticorena et al, 1997 (unitless)
 485 ! Arguably more ~consistent~ fudge than alpha, which has many walnuts
 486 ! sprinkled throughout the literature.
 487
 488   REAL, PARAMETER :: cmb=1.0    ! Marticorena et al,1997
 489 ! REAL, PARAMETER :: cmb=2.61   ! White,1979
 490
 491 ! Parameters used in Kok distribution function. Advise not to play with
 492 ! these without the expressed written consent of someone who knows what
 493 ! they're doing. (See Kok, 2010 PNAS for details on this scheme).
 494
 495 ! REAL, PARAMETER :: mmd_dust=3.4D-6  ! median mass diameter (m)
 496 ! REAL, PARAMETER :: gsd_dust=3.0     ! geometric standard deviation
 497 ! REAL, PARAMETER :: lambda=12.0D-6   ! crack propogation length (m)
 498 ! REAL, PARAMETER :: cv=12.62D-6      ! normalization constant
 499
 500 ! Calculate saltation surface area distribution from sand, silt, and clay
 501 ! mass fractions and saltation bin fraction. This will later become a
 502 ! modifier to the total saltation flux.  The reasoning here is that the
 503 ! size and availability of saltators affects saltation efficiency. Based
 504 ! on Eqn. (32) in Marticorena & Bergametti, 1995 (hereon, MB95).
 505
 506   DO n=1,smx
 507     dmass=massfrac(spoint(n))*frac_salt(n)
 508     dsurface(n)=0.75*dmass/(den_salt(n)*reff_salt(n))
 509   ENDDO
 510
 511 ! The following equation yields relative surface area fraction.  It will only
 512 ! work if you are representing the "full range" of all three soil classes.
 513 ! For this reason alone, we have incorporated particle sizes that encompass
 514 ! the clay class, to account for its relative area over the basal
 515 ! surface, even though these smaller bins would be unlikely to play any large
 516 ! role in the actual saltation process.
 517
 518   stotal=SUM(dsurface(:))
 519   DO n=1,smx
 520     ds_rel(n)=dsurface(n)/stotal
 521   ENDDO
 522
 523 ! Calculate total dust emission due to saltation of sand sized particles.
 524 ! Begin by calculating DRY threshold friction velocity (u_ts0).  Next adjust
 525 ! u_ts0 for moisture to get threshold friction velocity (u_ts). Then
 526 ! calculate saltation flux (salt) where ustar has exceeded u_ts.  Finally,
 527 ! calculate total dust emission (tot_emit), taking into account erodibility.
 528
 529   g = g0*1.0E2                          ! (cm s^-2)
 530   emit=0.0
 531
 532   DO n = 1, smx                         ! Loop over saltation bins
 533     den(n) = den_salt(n)*1.0D-3         ! (g cm^-3)
 534     diam(n) = 2.0*reff_salt(n)*1.0D2    ! (cm)
 535     DO i = 1,imx
 536       DO j = 1,jmx
 537         rhoa = airden(i,j,1)*1.0D-3     ! (g cm^-3)
 538
 539         ! Threshold friction velocity as a function of the dust density and
 540         ! diameter from Bagnold (1941) (m s^-1).
 541
 542         u_ts0 = 0.13*1.0D-2*SQRT(den(n)*g*diam(n)/rhoa)* &
 543                 SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ &
 544                 SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0)
 545
 546         ! Friction velocity threshold correction function based on physical
 547         ! properties related to moisture tension. Soil moisture greater than
 548         ! dry limit serves to increase threshold friction velocity (making
 549         ! it more difficult to loft dust). When soil moisture has not reached
 550         ! dry limit, treat as dry (no correction to threshold friction
 551         ! velocity). GAC
 552
 553         ! Calculate threshold friction velocity. If volumetric (gravimetric)
 554         ! water content is less than the drylimit, then the threshold friction
 555         ! velocity (u_ts) will be equal to the dry threshold friction velocity
 556         ! (u_ts0). EDH
 557
 558         IF (smois_opt .EQ. 1) THEN
 559           IF (100.*volsm > drylimit(i,j)) THEN
 560             u_ts = MAX(0.0D+0,u_ts0*(sqrt(1.0+1.21*((100.*volsm)-drylimit(i,j))**0.68)))
 561           ELSE
 562             u_ts = u_ts0
 563           ENDIF
 564         ELSE
 565           IF (gravsm(i,j) > drylimit(i,j)) THEN
 566             u_ts = MAX(0.0D+0,u_ts0*(sqrt(1.0+1.21*(gravsm(i,j)-drylimit(i,j))**0.68)))
 567           ELSE
 568             u_ts = u_ts0
 569           END IF
 570         END IF
 571
 572         ! Bin 7 threshold friction velocity for diagnostic dust lofting
 573         ! potential product
 574
 575         IF (n .eq. 7) THEN  ! Saltation bin 7 is small sand
 576           ustart = u_ts
 577         ENDIF
 578
 579         ! Saltation flux (kg m^-1 s^-1) from MB95
 580         ! ds_rel is the relative surface area distribution
 581
 582         IF (ustar(i,j) .gt. u_ts .and. erod(i,j) .gt. 0.0 .and. ilwi(i,j) == 1) THEN
 583           salt = cmb*ds_rel(n)*(airden(i,j,1)/g0)*(ustar(i,j)**3)* &
 584                  (1. + u_ts/ustar(i,j))*(1. - (u_ts**2)/(ustar(i,j)**2))
 585         ELSE
 586           salt = 0.D0
 587         ENDIF
 588
 589         ! Calculate total vertical mass flux (note beta has units of m^-1)
 590         ! Beta acts to tone down dust in areas with so few dust-sized particles
 591         ! that the lofting efficiency decreases.  Otherwise, super sandy zones
 592         ! would be huge dust producers, which is generally not the case.
 593         ! Equation derived from wind-tunnel experiments (see MB95).
 594
 595         beta=10**(13.6*massfrac(1)-6.0)  ! (unitless)
 596         IF (beta .gt. betamax) THEN
 597           beta=betamax
 598         ENDIF
 599         emit=emit+salt*(erod(i,j)**gamma)*alpha*beta    ! (kg m^-2 s^-1)
 600       END DO
 601     END DO
 602   END DO                                 ! End do over saltation bins
 603
 604   ! Now that we have the total dust emission, distribute into dust bins using
 605   ! lognormal distribution (Dr. Jasper Kok, 2010), and
 606   ! calculate total mass emitted over the grid box over the timestep.
 607   !
 608   ! In calculating the Kok distribution, we assume upper and lower limits to
 609   ! each bin. For reff_dust=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/) (default),
 610   ! lower limits were ASSUMED at lo_dust=(/0.1D-6,1.0D-6,1.8D-6,3.0D-6,6.0D-6/)
 611   ! upper limits were ASSUMED at up_dust=(/1.0D-6,1.8D-6,3.0D-6,6.0D-6,10.0D-6/)
 612   ! These may be changed within module_data_gocart_dust.F, but make sure it is
 613   ! consistent with reff_dust values.  These values were taken from the original
 614   ! GOCART bin configuration. We use them here to calculate dust bin width,
 615   ! dlndp. dvol is the volume distribution. GAC
 616   !
 617   ! UPDATE: We bypass the calculation below and instead hardcode distr_dust for
 618   ! the five dust bins we are using here since this distribution is static and
 619   ! unnecessary to calculate at every time step. Keeping everything here to
 620   ! document the steps in obtaining distr_dust.  GAC 20140320
 621
 622 !  dvol_tot=0.
 623 !  DO n=1,nmx  ! Loop over all dust bins
 624 !    dlndp(n)=LOG(up_dust(n)/lo_dust(n))
 625 !    dvol(n)=(2.0*reff_dust(n)/cv)*(1.+ERF(LOG(2.0*reff_dust(n)/mmd_dust)/(SQRT(2.)*LOG(gsd_dust))))*&
 626 !          EXP(-(2.0*reff_dust(n)/lambda)**3.0)*dlndp(n)
 627 !    dvol_tot=dvol_tot+dvol(n)
 628 !
 629 !   ! Convert mass flux to volume flux
 630 !    emit_vol=emit/den_dust(n) ! (m s^-1)
 631 !  END DO
 632 !  DO n=1,nmx  ! Loop over all dust bins
 633 !    distr_dust(n)=dvol(n)/dvol_tot
 634 !  END DO
 635
 636   ! Now distribute total vertical emission into dust bins and update
 637   ! concentration.
 638
 639   DO n=1,nmx  ! Loop over all dust bins
 640     DO i=1,imx
 641       DO j=1,jmx
 642
 643         ! Calculate total mass emitted
 644
 645         dsrc = emit*distr_dust(n)*dt1  ! (kg m^-2) per dt1
 646         IF (dsrc < 0.0) dsrc = 0.0
 647
 648         ! Update dust mixing ratio at first model level.
 649         tc(i,j,1,n) = tc(i,j,1,n) + dsrc/dz_lowest/airden(i,j,1)  ! (kg/kg)
 650         bems(i,j,n) = dsrc                     ! diagnostic (kg/m2) per dt1
 651       END DO
 652     END DO
 653   END DO
 654
 655   END SUBROUTINE source_dust
 656
 657
 658 END MODULE GOCART_DUST_AFWA