chem/module_gocart_dust.F

   1 MODULE GOCART_DUST
   2
   3 ! A. Ukhov, 11 March 2021, Now "emis_dust" is accumulated dust
   4 ! emission (kg/m2). Before was instantenious flux (kg/cell).
   5 ! Bug fix in the loop over cells: Cells near domain boundaries
   6 ! were not processed. Code cleanup, remove unused variables.
   7
   8   USE module_data_gocart_dust
   9
  10 CONTAINS
  11   subroutine gocart_dust_driver(dt,config_flags,alt,t_phy,u_phy,          &
  12          v_phy,chem,rho_phy,dz8w,smois,u10,v10,erod,dustin,               &
  13          isltyp,xland,g,emis_dust,                                        &
  14          ids,ide, jds,jde, kds,kde,                                       &
  15          ims,ime, jms,jme, kms,kme,                                       &
  16          its,ite, jts,jte, kts,kte                                        )
  17   USE module_configure
  18   USE module_state_description
  19   IMPLICIT NONE
  20    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  21
  22    INTEGER,      INTENT(IN   ) :: ids,ide, jds,jde, kds,kde,               &
  23                                   ims,ime, jms,jme, kms,kme,               &
  24                                   its,ite, jts,jte, kts,kte
  25
  26    INTEGER,DIMENSION( ims:ime , jms:jme ),                                 &
  27           INTENT(IN   ) ::                           isltyp
  28    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
  29          INTENT(INOUT ) ::                           chem
  30    REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL,&
  31          INTENT(INOUT ) ::                           emis_dust
  32   REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) ,      &
  33       INTENT(INOUT) ::                               smois
  34    REAL,  DIMENSION( ims:ime , jms:jme, 3 )                   ,            &
  35           INTENT(IN   ) ::    erod
  36    REAL,  DIMENSION( ims:ime , jms:jme, 5 )                   ,            &
  37           INTENT(INout   ) ::    dustin
  38    REAL,  DIMENSION( ims:ime , jms:jme )                   ,               &
  39           INTENT(IN ) ::    u10,v10,xland
  40    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  41           INTENT(IN ) ::                                alt,               &
  42                                                       t_phy,               &
  43                                                      dz8w,                 &
  44                                               u_phy,v_phy,rho_phy
  45
  46   REAL, INTENT(IN   ) :: dt,g
  47 !
  48 ! local variables
  49 !
  50   integer :: nmx,i,j,k,imx,jmx,lmx
  51   integer,dimension (1,1) :: ilwi
  52   real*8, DIMENSION (1,1,3,1) :: erodin
  53   real*8, DIMENSION (5) :: tc,bems
  54   real*8, dimension (1,1) :: w10m,gwet,airden
  55   real                    :: dz_lowest
  56   real*8  conver,converi
  57   conver=1.e-9
  58   converi=1.e9
  59 !
  60 ! number of dust bins
  61 !
  62   imx=1
  63   jmx=1
  64   lmx=1
  65   nmx=5
  66   k=kts
  67   DO j=jts,jte
  68   DO i=its,ite
  69
  70 ! no dust over water!!!
  71 !
  72      if(xland(i,j).lt.1.5)then
  73      ilwi(1,1)=1
  74      if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
  75         tc(:)=1.e-16*conver
  76      else
  77         tc(1)=chem(i,kts,j,p_dust_1)*conver
  78         tc(2)=chem(i,kts,j,p_dust_2)*conver
  79         tc(3)=chem(i,kts,j,p_dust_3)*conver
  80         tc(4)=chem(i,kts,j,p_dust_4)*conver
  81         tc(5)=chem(i,kts,j,p_dust_5)*conver
  82      endif
  83      w10m(1,1)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
  84 !
  85 ! don't trust the u10,v10 values, is model layers are very thin near surface
  86 !
  87      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))
  88      erodin(1,1,1,1)=erod(i,j,1)
  89      erodin(1,1,2,1)=erod(i,j,2)
  90      erodin(1,1,3,1)=erod(i,j,3)
  91 !
  92 !  volumetric soil moisture over porosity
  93 !
  94      gwet(1,1)=smois(i,1,j)/porosity(isltyp(i,j))
  95      airden(1,1)=rho_phy(i,kts,j)
  96      dz_lowest = dz8w(i,1,j)
  97
  98     call source_du( imx,jmx,lmx,nmx, dt, tc, &
  99                      erodin, ilwi, w10m, gwet, airden, &
 100                      dz_lowest,bems,config_flags%start_month,g)
 101
 102     if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
 103         dustin(i,j,1:5)=tc(1:5)*converi
 104     else
 105         chem(i,kts,j,p_dust_1)=tc(1)*converi  ! tc(1...5) is (kg/kg), p_dust_1...5 (ug/kg)
 106         chem(i,kts,j,p_dust_2)=tc(2)*converi
 107         chem(i,kts,j,p_dust_3)=tc(3)*converi
 108         chem(i,kts,j,p_dust_4)=tc(4)*converi
 109         chem(i,kts,j,p_dust_5)=tc(5)*converi
 110     endif
 111
 112      ! A. Ukhov
 113      ! for output diagnostics
 114      ! bems (kg/m2) per dt
 115      ! p_edust1...5 is accumulated dust emission (kg/m2)
 116      emis_dust(i,1,j,p_edust1)=emis_dust(i,1,j,p_edust1)+bems(1)
 117      emis_dust(i,1,j,p_edust2)=emis_dust(i,1,j,p_edust2)+bems(2)
 118      emis_dust(i,1,j,p_edust3)=emis_dust(i,1,j,p_edust3)+bems(3)
 119      emis_dust(i,1,j,p_edust4)=emis_dust(i,1,j,p_edust4)+bems(4)
 120      emis_dust(i,1,j,p_edust5)=emis_dust(i,1,j,p_edust5)+bems(5)
 121      endif
 122   enddo
 123   enddo
 124 !
 125
 126 end subroutine gocart_dust_driver
 127
 128
 129   SUBROUTINE source_du( imx,jmx,lmx,nmx, dt1, tc, &
 130                      erod, ilwi, w10m, gwet, airden, &
 131                      dz_lowest,bems,month,g0)
 132
 133 ! ****************************************************************************
 134 ! *  Evaluate the source of each dust particles size classes  (kg/m3)
 135 ! *  by soil emission.
 136 ! *  Input:
 137 ! *         EROD      Fraction of erodible grid cell                (-)
 138 ! *                   for 1: Sand, 2: Silt, 3: Clay
 139 ! *         AIRVOL    Volume occupy by each grid boxes              (m3)
 140 ! *         DT1       Time step                                      (s)
 141 ! *         W10m      Velocity at the anemometer level (10meters)   (m/s)
 142 ! *         u_tresh   Threshold velocity for particule uplifting    (m/s)
 143 ! *         CH_dust   Constant to fudge the total emission of dust  (s2/m2)
 144 ! *         dz_lowest heigth of the lowest layer                     (m)
 145 ! *
 146 ! *  Output:
 147 ! *         DSRC      Source of each dust type           (kg/timestep/m2)
 148 ! *         BEMS      Source of each dust type           (kg/timestep/m2)
 149 ! *
 150 ! *  Working:
 151 ! *         SRC       Potential source                   (kg/m/timestep/cell)
 152 ! *
 153 ! ****************************************************************************
 154
 155   INTEGER,   INTENT(IN)    :: nmx,imx,jmx,lmx
 156   REAL*8,    INTENT(IN)    :: erod(imx,jmx,ndcls,ndsrc)
 157   INTEGER,   INTENT(IN)    :: ilwi(imx,jmx),month
 158
 159   REAL*8,    INTENT(IN)    :: w10m(imx,jmx), gwet(imx,jmx)
 160   REAL*8,    INTENT(IN)    :: airden(imx,jmx,lmx)
 161   REAL*8,    INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 162   REAL*8,    INTENT(OUT)   :: bems(imx,jmx,nmx)
 163   REAL,      INTENT(IN  )  :: dz_lowest
 164
 165   REAL*8    :: den(nmx), diam(nmx)
 166   REAL*8    :: u_ts0, u_ts, dsrc, srce
 167   REAL, intent(in)    :: g0
 168   REAL    :: rhoa, g,dt1
 169   INTEGER :: i, j, n, m, k
 170
 171   ! executable statemenst
 172
 173   DO n = 1, nmx
 174      ! Threshold velocity as a function of the dust density and the diameter
 175      ! from Bagnold (1941)
 176      den(n) = den_dust(n)*1.0D-3
 177      diam(n) = 2.0*reff_dust(n)*1.0D2
 178      g = g0*1.0E2
 179      ! Pointer to the 3 classes considered in the source data files
 180      m = ipoint(n)
 181      DO k = 1, ndsrc
 182         ! No flux if wet soil
 183         DO i = 1,imx
 184            DO j = 1,jmx
 185               rhoa = airden(i,j,1)*1.0D-3
 186               u_ts0 = 0.13*1.0D-2*SQRT(den(n)*g*diam(n)/rhoa)* &
 187                    SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ &
 188                    SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0)
 189 !             write(0,*)u_ts0,den(n),diam(n),rhoa,g
 190
 191               ! Case of surface dry enough to erode
 192               IF (gwet(i,j) < 0.5) THEN
 193                  u_ts = MAX(0.0D+0,u_ts0*(1.2D+0+2.0D-1*LOG10(MAX(1.0D-3, gwet(i,j)))))
 194               ELSE
 195                  ! Case of wet surface, no erosion
 196                  u_ts = 100.0
 197               END IF
 198               srce = frac_s(n)*erod(i,j,m,k)   ! (kg s^2 m^-5)
 199               IF (ilwi(i,j) == 1 ) THEN
 200                                     !(kg s^2 m^-5)*(m^3 s^-3)*s = (kg/m2) per dt1
 201                  dsrc = ch_dust(n,month)*srce*w10m(i,j)**2 * (w10m(i,j) - u_ts)*dt1
 202               ELSE
 203                  dsrc = 0.0
 204               END IF
 205               IF (dsrc < 0.0) dsrc = 0.0
 206
 207               ! Update dust mixing ratio at first model level.
 208               tc(i,j,1,n) = tc(i,j,1,n) + dsrc/dz_lowest/airden(i,j,1)  ! (kg/kg)
 209               bems(i,j,n) = dsrc                     ! diagnostic (kg/m2) per dt1
 210            END DO
 211         END DO
 212      END DO
 213   END DO
 214
 215 END SUBROUTINE source_du
 216
 217
 218 END MODULE GOCART_DUST