wrfv2_fire/chem/module_gocart_dust.F

   1 MODULE GOCART_DUST
   2
   3
   4   USE module_data_gocart_dust
   5
   6 CONTAINS
   7   subroutine gocart_dust_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy,  &
   8          v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,                  &
   9          ivgtyp,isltyp,vegfra,xland,xlat,xlong,gsw,dx,g,emis_dust,        &
  10          ids,ide, jds,jde, kds,kde,                                        &
  11          ims,ime, jms,jme, kms,kme,                                        &
  12          its,ite, jts,jte, kts,kte                                         )
  13   USE module_configure
  14   USE module_state_description
  15   USE module_model_constants, ONLY: mwdry
  16   IMPLICIT NONE
  17    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  18
  19    INTEGER,      INTENT(IN   ) :: julday, ktau,                     &
  20                                   ids,ide, jds,jde, kds,kde,               &
  21                                   ims,ime, jms,jme, kms,kme,               &
  22                                   its,ite, jts,jte, kts,kte
  23    INTEGER,DIMENSION( ims:ime , jms:jme )                  ,               &
  24           INTENT(IN   ) ::                                                 &
  25                                                      ivgtyp,               &
  26                                                      isltyp
  27    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),                &
  28          INTENT(IN ) ::                                   moist
  29    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
  30          INTENT(INOUT ) ::                                   chem
  31    REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL,&
  32          INTENT(INOUT ) ::                                                 &
  33          emis_dust
  34   REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) ,      &
  35       INTENT(INOUT) ::                               smois
  36    REAL,  DIMENSION( ims:ime , jms:jme, 3 )                   ,               &
  37           INTENT(IN   ) ::    erod
  38    REAL,  DIMENSION( ims:ime , jms:jme )                   ,               &
  39           INTENT(IN   ) ::                                                 &
  40                                                      u10,                  &
  41                                                      v10,                  &
  42                                                      gsw,                  &
  43                                                   vegfra,                  &
  44                                                      xland,                &
  45                                                      xlat,                 &
  46                                                      xlong
  47    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  48           INTENT(IN   ) ::                                                 &
  49                                                         alt,               &
  50                                                       t_phy,               &
  51                                                      dz8w,p8w,             &
  52                                               u_phy,v_phy,rho_phy
  53
  54   REAL, INTENT(IN   ) :: dt,dx,g
  55 !
  56 ! local variables
  57 !
  58   integer :: nmx,i,j,k,ndt,imx,jmx,lmx
  59   integer,dimension (1,1) :: ilwi
  60   real*8, DIMENSION (1,1,3,1) :: erodin
  61   real*8, DIMENSION (5) :: tc,bems
  62   real*8, dimension (1,1) :: w10m,gwet,airden,airmas
  63   real*8, dimension (1) :: dxy
  64   real*8  conver,converi
  65   real dttt
  66 ! conver=1.e-9*mwdry
  67 ! converi=1.e9/mwdry
  68   conver=1.e-9
  69   converi=1.e9
  70 !
  71 ! number of dust bins
  72 !
  73   imx=1
  74   jmx=1
  75   lmx=1
  76   nmx=5
  77   k=kts
  78   do j=jts,jte
  79   do i=its,ite
  80 !
  81 ! donṫ do dust over water!!!
  82 !
  83      if(xland(i,j).lt.1.5)then
  84      ilwi(1,1)=1
  85      tc(1)=chem(i,kts,j,p_dust_1)*conver
  86      tc(2)=chem(i,kts,j,p_dust_2)*conver
  87      tc(3)=chem(i,kts,j,p_dust_3)*conver
  88      tc(4)=chem(i,kts,j,p_dust_4)*conver
  89      tc(5)=chem(i,kts,j,p_dust_5)*conver
  90      w10m(1,1)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
  91      airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g
  92 !
  93 ! we donṫ trust the u10,v10 values, is model layers are very thin near surface
  94 !
  95      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))
  96      erodin(1,1,1,1)=erod(i,j,1)!/dx/dx
  97      erodin(1,1,2,1)=erod(i,j,2)!/dx/dx
  98      erodin(1,1,3,1)=erod(i,j,3)!/dx/dx
  99 !
 100 !  volumetric soil moisture over porosity
 101 !
 102      gwet(1,1)=smois(i,1,j)/porosity(isltyp(i,j))
 103 !    gwet(1,1)=.1
 104      ndt=ifix(dt)
 105      airden(1,1)=rho_phy(i,kts,j)
 106      dxy(1)=dx*dx
 107 !    erodin(1,1,1,1)=  0.149748762553862
 108 !    erodin(1,1,2,1)=  7.487438878070708E-002
 109 !    erodin(1,1,3,1)=  7.487438878070708E-002
 110 !    ilwi(1,1)=           1
 111 !    dxy(1)=   54585850453.7552
 112 !    w10m(1,1)=   10.6305338763678
 113 !    gwet(1,1)=  9.136307984590530E-002
 114 !    airden(1,1)=   1.16423276395132
 115 !    airmas(1,1)=   8114017750938.79
 116 !    tc (1) =  1.000000000000000D-030
 117 !    tc (2) =  1.000000000000000d-030
 118 !    tc (3) =  1.000000000000000d-030
 119 !    tc(4)  =  1.000000000000000d-030
 120 !    tc(5)  =  1.000000000000000d-030
 121 !    dttt=3600.
 122     call source_du( imx,jmx,lmx,nmx, dt, tc, &
 123                      erodin, ilwi, dxy, w10m, gwet, airden, airmas, &
 124                      bems,config_flags%start_month,g)
 125 !    write(0,*)tc(1)
 126 !    write(0,*)tc(2)
 127 !    write(0,*)tc(3)
 128 !    write(0,*)tc(4)
 129 !    write(0,*)tc(5)
 130      chem(i,kts,j,p_dust_1)=tc(1)*converi
 131      chem(i,kts,j,p_dust_2)=tc(2)*converi
 132      chem(i,kts,j,p_dust_3)=tc(3)*converi
 133      chem(i,kts,j,p_dust_4)=tc(4)*converi
 134      chem(i,kts,j,p_dust_5)=tc(5)*converi
 135 ! for output diagnostics
 136      emis_dust(i,1,j,p_edust1)=bems(1)
 137      emis_dust(i,1,j,p_edust2)=bems(2)
 138      emis_dust(i,1,j,p_edust3)=bems(3)
 139      emis_dust(i,1,j,p_edust4)=bems(4)
 140      emis_dust(i,1,j,p_edust5)=bems(5)
 141      endif
 142   enddo
 143   enddo
 144 !
 145
 146 end subroutine gocart_dust_driver
 147
 148
 149   SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, &
 150                      erod, ilwi, dxy, w10m, gwet, airden, airmas, &
 151                      bems,month,g0)
 152
 153 ! ****************************************************************************
 154 ! *  Evaluate the source of each dust particles size classes  (kg/m3)
 155 ! *  by soil emission.
 156 ! *  Input:
 157 ! *         EROD      Fraction of erodible grid cell                (-)
 158 ! *                   for 1: Sand, 2: Silt, 3: Clay
 159 ! *         DUSTDEN   Dust density                                  (kg/m3)
 160 ! *         DXY       Surface of each grid cell                     (m2)
 161 ! *         AIRVOL    Volume occupy by each grid boxes              (m3)
 162 ! *         NDT1      Time step                                     (s)
 163 ! *         W10m      Velocity at the anemometer level (10meters)   (m/s)
 164 ! *         u_tresh   Threshold velocity for particule uplifting    (m/s)
 165 ! *         CH_dust   Constant to fudge the total emission of dust  (s2/m2)
 166 ! *
 167 ! *  Output:
 168 ! *         DSRC      Source of each dust type           (kg/timestep/cell)
 169 ! *
 170 ! *  Working:
 171 ! *         SRC       Potential source                   (kg/m/timestep/cell)
 172 ! *
 173 ! ****************************************************************************
 174
 175 ! USE module_data_gocart
 176 ! USE module_data_gocart_dust
 177
 178
 179
 180   INTEGER, INTENT(IN)    :: nmx,imx,jmx,lmx
 181   REAL*8,    INTENT(IN)    :: erod(imx,jmx,ndcls,ndsrc)
 182   INTEGER, INTENT(IN)    :: ilwi(imx,jmx),month
 183
 184   REAL*8,    INTENT(IN)    :: w10m(imx,jmx), gwet(imx,jmx)
 185   REAL*8,    INTENT(IN)    :: dxy(jmx)
 186   REAL*8,    INTENT(IN)    :: airden(imx,jmx,lmx), airmas(imx,jmx,lmx)
 187   REAL*8,    INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 188   REAL*8,    INTENT(OUT)   :: bems(imx,jmx,nmx)
 189
 190   REAL*8    :: den(nmx), diam(nmx)
 191   REAL*8    :: tsrc, u_ts0, cw, u_ts, dsrc, srce
 192   REAL, intent(in)    :: g0
 193   REAL    :: rhoa, g,dt1
 194   INTEGER :: i, j, n, m, k
 195
 196
 197   REAL*8                  :: tcmw(nmx), ar(nmx), tcvv(nmx)
 198   REAL*8                  :: ar_wetdep(nmx), kc(nmx)
 199   CHARACTER(LEN=20)     :: tcname(nmx), tcunits(nmx)
 200   LOGICAL               :: aerosol(nmx)
 201
 202
 203 ! REAL*8 :: tc1(imx,jmx,lmx,nmx)
 204 ! REAL*8, TARGET :: tcms(imx,jmx,lmx,nmx) ! tracer mass (kg; kgS for sulfur case)
 205 ! REAL*8, TARGET :: tcgm(imx,jmx,lmx,nmx) ! g/m3
 206
 207   !-----------------------------------------------------------------------
 208   ! sea salt specific
 209   !-----------------------------------------------------------------------
 210 ! REAL*8, DIMENSION(nmx) :: ssaltden, ssaltreff, ra, rb
 211 ! REAL*8 :: ch_ss(nmx,12)
 212
 213   !-----------------------------------------------------------------------
 214   ! emissions (input)
 215   !-----------------------------------------------------------------------
 216 ! REAL*8 :: e_an(imx,jmx,2,nmx), e_bb(imx,jmx,nmx), &
 217 !         e_ac(imx,jmx,lmx,nmx)
 218
 219   !-----------------------------------------------------------------------
 220   ! diagnostics (budget)
 221   !-----------------------------------------------------------------------
 222 !  ! tendencies per time step and process
 223 !  REAL, TARGET :: bems(imx,jmx,nmx), bdry(imx,jmx,nmx), bstl(imx,jmx,nmx)
 224 !  REAL, TARGET :: bwet(imx,jmx,nmx), bcnv(imx,jmx,nmx)
 225 !
 226 !  ! integrated tendencies per process
 227 !  REAL, TARGET :: tems(imx,jmx,nmx), tstl(imx,jmx,nmx)
 228 !  REAL, TARGET :: tdry(imx,jmx,nmx), twet(imx,jmx,nmx), tcnv(imx,jmx,nmx)
 229
 230   ! global mass balance per time step
 231   REAL*8 :: tmas0(nmx), tmas1(nmx)
 232   REAL*8 :: dtems(nmx), dttrp(nmx), dtdif(nmx), dtcnv(nmx)
 233   REAL*8 :: dtwet(nmx), dtdry(nmx), dtstl(nmx)
 234   REAL*8 :: dtems2(nmx), dttrp2(nmx), dtdif2(nmx), dtcnv2(nmx)
 235   REAL*8 :: dtwet2(nmx), dtdry2(nmx), dtstl2(nmx)
 236
 237
 238
 239
 240   ! executable statemenst
 241
 242   DO n = 1, nmx
 243      ! Threshold velocity as a function of the dust density and the diameter
 244      ! from Bagnold (1941)
 245      den(n) = den_dust(n)*1.0D-3
 246      diam(n) = 2.0*reff_dust(n)*1.0D2
 247      g = g0*1.0E2
 248      ! Pointer to the 3 classes considered in the source data files
 249      m = ipoint(n)
 250      tsrc = 0.0
 251      DO k = 1, ndsrc
 252         ! No flux if wet soil
 253         DO i = 1,imx
 254            DO j = 1,jmx
 255               rhoa = airden(i,j,1)*1.0D-3
 256               u_ts0 = 0.13*1.0D-2*SQRT(den(n)*g*diam(n)/rhoa)* &
 257                    SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ &
 258                    SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0)
 259 !             write(0,*)u_ts0,den(n),diam(n),rhoa,g
 260               ! Fraction of emerged surfaces (subtract lakes, coastal ocean,..)
 261 !              cw = 1.0 - water(i,j)
 262
 263               ! Case of surface dry enough to erode
 264               IF (gwet(i,j) < 0.2) THEN
 265 !              IF (gwet(i,j) < 0.5) THEN  !  Pete's modified value
 266                  u_ts = MAX(0.0D+0,u_ts0*(1.2D+0+2.0D-1*LOG10(MAX(1.0D-3, gwet(i,j)))))
 267               ELSE
 268                  ! Case of wet surface, no erosion
 269                  u_ts = 100.0
 270               END IF
 271               srce = frac_s(n)*erod(i,j,m,k)*dxy(j)  ! (m2)
 272               IF (ilwi(i,j) == 1 ) THEN
 273                  dsrc = ch_dust(n,month)*srce*w10m(i,j)**2 &
 274                       * (w10m(i,j) - u_ts)*dt1  ! (kg)
 275 !                write(0,*)ch_dust(n,month),srce,w10m(i,j),u_ts,gwet(i,j)
 276               ELSE
 277                  dsrc = 0.0
 278               END IF
 279 !              dsrc = cw*ch_dust(k)*srce*w10m(i,j)**2 &
 280 !                   * (w10m(i,j) - u_ts)*dt1  ! (kg)
 281 !              dsrc = cw*ch_dust(n,dt(1)%mn)*srce*w10m(i,j)**2 &
 282 !                   * (w10m(i,j) - u_ts)*dt1  ! (kg)
 283               IF (dsrc < 0.0) dsrc = 0.0
 284
 285               ! Update dust mixing ratio at first model level.
 286               tc(i,j,1,n) = tc(i,j,1,n) + dsrc / airmas(i,j,1)
 287               bems(i,j,n) = dsrc
 288            END DO
 289         END DO
 290      END DO
 291   END DO
 292
 293 END SUBROUTINE source_du
 294
 295
 296 END MODULE GOCART_DUST