chem/module_vash_settling.F

   1 MODULE MODULE_VASH_SETTLING
   2
   3 CONTAINS
   4
   5 SUBROUTINE vash_settling_driver(dt,config_flags,t_phy,moist,               &
   6          chem,rho_phy,dz8w,p8w,p_phy,                                      &
   7          ash_fall,dx,g,                                                    &
   8          ids,ide, jds,jde, kds,kde,                                        &
   9          ims,ime, jms,jme, kms,kme,                                        &
  10          its,ite, jts,jte, kts,kte                                         )
  11   USE module_configure
  12   USE module_state_description
  13 ! USE module_data_gocart_dust
  14 ! USE module_data_gocart_seas
  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   ) ::                      &
  20                                   ids,ide, jds,jde, kds,kde,               &
  21                                   ims,ime, jms,jme, kms,kme,               &
  22                                   its,ite, jts,jte, kts,kte
  23     REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),                &
  24          INTENT(IN ) ::                                   moist
  25    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
  26          INTENT(INOUT ) ::                                   chem
  27    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  28           INTENT(IN   ) ::  t_phy,p_phy,dz8w,p8w,rho_phy
  29    REAL,  DIMENSION( ims:ime , jms:jme ),                        &
  30           INTENT(INOUT   ) ::  ash_fall
  31
  32   REAL, INTENT(IN   ) :: dt,dx,g
  33   integer :: nmx,i,j,k,kk,lmx,iseas,idust
  34   real*8, DIMENSION (1,1,kte-kts+1) :: tmp,airden,airmas,p_mid,delz,rh
  35   real*8, DIMENSION (1,1,kte-kts+1,4) :: sea_salt
  36 !srf
  37   real*8, DIMENSION (1,1,kte-kts+1,10) :: ash
  38   real*8, DIMENSION (10), PARAMETER :: den_ash(10)=(/2500.,2500.,2500.,2500.,2500., &
  39                                                      2500.,2500.,2500.,2500.,2500. /)
  40   real*8, DIMENSION (10), PARAMETER :: reff_ash(10)=(/0.5000D-3,&! 1.00 mm diameter
  41                                                       0.3750D-3,&! 0.75 mm
  42                                                       0.1875D-3,&!
  43                                                       93.750D-6,&!
  44                                                       46.875D-6,&!
  45                                                       23.437D-6,&!
  46                                                       11.719D-6,&!
  47                                                       05.859D-6,&!
  48                                                       02.930D-6,&!
  49                                                       00.975D-6 /)! 3.9 um
  50   real*8, DIMENSION (10) :: bstl_ash
  51   integer iash
  52 !srf
  53
  54 !
  55 ! bstl is for budgets
  56 !
  57
  58   real*8 conver,converi
  59        conver=1.e-9
  60        converi=1.e9
  61        lmx=kte-kts+1
  62        do j=jts,jte
  63        do i=its,ite
  64           kk=0
  65           bstl_ash(:)=0.
  66           ash(:,:,:,:)=0.
  67           do k=kts,kte
  68           kk=kk+1
  69           p_mid(1,1,kk)=.01*p_phy(i,kte-k+kts,j)
  70           delz(1,1,kk)=dz8w(i,kte-k+kts,j)
  71           airmas(1,1,kk)=-(p8w(i,k+1,j)-p8w(i,k,j))/g
  72           airden(1,1,kk)=rho_phy(i,k,j)
  73           tmp(1,1,kk)=t_phy(i,k,j)
  74           rh(1,1,kk) = .95
  75           rh(1,1,kk) = MIN( .95, moist(i,k,j,p_qv) / &
  76                (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
  77                (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))))
  78           rh(1,1,kk)=max(1.0D-1,rh(1,1,kk))
  79           enddo
  80
  81 !ash settling
  82
  83           iseas=0
  84           idust=0
  85           iash =1
  86
  87           kk=0
  88 !         write(0,*)'1',chem(i,1,j,p_dust_4)
  89           do k=kts,kte
  90           kk=kk+1
  91           ash(1,1,kk,7)=chem(i,k,j,p_vash_7)*conver
  92           ash(1,1,kk,8)=chem(i,k,j,p_vash_8)*conver
  93           ash(1,1,kk,9)=chem(i,k,j,p_vash_9)*conver
  94           ash(1,1,kk,10)=chem(i,k,j,p_vash_10)*conver
  95           enddo
  96           if(config_flags%chem_opt == 400  .or. config_flags%chem_opt == 402 ) then
  97           kk=0
  98           do k=kts,kte
  99           kk=kk+1
 100           ash(1,1,kk,1)=chem(i,k,j,p_vash_1)*conver
 101           ash(1,1,kk,2)=chem(i,k,j,p_vash_2)*conver
 102           ash(1,1,kk,3)=chem(i,k,j,p_vash_3)*conver
 103           ash(1,1,kk,4)=chem(i,k,j,p_vash_4)*conver
 104           ash(1,1,kk,5)=chem(i,k,j,p_vash_5)*conver
 105           ash(1,1,kk,6)=chem(i,k,j,p_vash_6)*conver
 106           enddo
 107           endif
 108
 109           call vsettling(1, 1, lmx, 10, g,&
 110                     ash, tmp, p_mid, delz, airmas, &
 111                     den_ash, reff_ash, dt, bstl_ash, rh, idust, iseas,iash)
 112           kk=0
 113           ash_fall(i,j)=ash_fall(i,j)+sum(bstl_ash(1:10))
 114           do k=kts,kte
 115           kk=kk+1
 116             chem(i,k,j,p_vash_7)=ash(1,1,kk,7)*converi
 117             chem(i,k,j,p_vash_8)=ash(1,1,kk,8)*converi
 118             chem(i,k,j,p_vash_9)=ash(1,1,kk,9)*converi
 119             chem(i,k,j,p_vash_10)=ash(1,1,kk,10)*converi
 120           enddo
 121           if(config_flags%chem_opt == 400  .or. config_flags%chem_opt == 402 ) then
 122           kk=0
 123           do k=kts,kte
 124           kk=kk+1
 125             chem(i,k,j,p_vash_1)=ash(1,1,kk,1)*converi
 126             chem(i,k,j,p_vash_2)=ash(1,1,kk,2)*converi
 127             chem(i,k,j,p_vash_3)=ash(1,1,kk,3)*converi
 128             chem(i,k,j,p_vash_4)=ash(1,1,kk,4)*converi
 129             chem(i,k,j,p_vash_5)=ash(1,1,kk,5)*converi
 130             chem(i,k,j,p_vash_6)=ash(1,1,kk,6)*converi
 131           enddo
 132           endif
 133
 134 !ash settling end
 135
 136        enddo
 137        enddo
 138 END SUBROUTINE vash_settling_driver
 139
 140
 141           subroutine vsettling(imx,jmx, lmx, nmx,g0, &
 142                     tc, tmp, p_mid, delz, airmas, &
 143                     den, reff, dt, bstl, rh, idust, iseas,iash)
 144 ! ****************************************************************************
 145 ! *                                                                          *
 146 ! *  Calculate the loss by settling, using an implicit method                *
 147 ! *                                                                          *
 148 ! *  Input variables:                                                        *
 149 ! *    SIGE(k)         - sigma coordinate of the vertical edges              *
 150 ! *    PS(i,j)         - Surface pressure (mb)                               *
 151 ! *    TMP(i,j,k)      - Air temperature  (K)                                *
 152 ! *    CT(i,j)         - Surface exchange coeff for moisture
 153 ! *                                                                          *
 154 ! ****************************************************************************
 155
 156
 157   IMPLICIT  NONE
 158
 159   INTEGER, INTENT(IN) :: imx, jmx, lmx, nmx,iseas,idust,iash
 160   INTEGER :: ntdt
 161   REAL, INTENT(IN) :: dt,g0 ! ,dyn_visc
 162   REAL*8,    INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx),  &
 163                          airmas(imx,jmx,lmx), rh(imx,jmx,lmx), &
 164                          den(nmx), reff(nmx), p_mid(imx,jmx,lmx)
 165   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 166   REAL*8, INTENT(OUT)   :: bstl(imx,jmx,nmx)
 167
 168   REAL*8    :: tc1(imx,jmx,lmx,nmx), dt_settl(nmx), rcm(nmx), rho(nmx)
 169   INTEGER :: ndt_settl(nmx)
 170   REAL*8    :: dzmin, vsettl, dtmax, pres, rhb, rwet(nmx), ratio_r(nmx)
 171   REAL*8    :: addmass,c_stokes, free_path, c_cun, viscosity, vd_cor, growth_fac
 172   REAL,    PARAMETER :: dyn_visc = 1.5E-5
 173   INTEGER :: k, n, i, j, l, l2
 174   ! for sea-salt:
 175   REAL*8, PARAMETER :: c1=0.7674, c2=3.079, c3=2.573E-11, c4=-1.424
 176
 177   ! for OMP:
 178   REAL*8 :: rwet_priv(nmx), rho_priv(nmx)
 179
 180   ! executable statements
 181
 182 ! IF (type) /= 'dust' .AND. TRIM(aero_type) /= 'sea_salt') RETURN
 183   if(idust.ne.1.and.iseas.ne.1.and.iash.ne.1)return
 184
 185   WHERE (tc(:,:,:,:) < 0.0) tc(:,:,:,:) = 1.0d-32
 186
 187   dzmin = MINVAL(delz(:,:,:))
 188   IF (idust == 1)     growth_fac = 1.0
 189   IF (iseas == 1)     growth_fac = 3.0
 190   IF (iash  == 1)     growth_fac = 1.0
 191
 192   DO k = 1,nmx
 193
 194      ! Settling velocity (m/s) for each tracer (Stokes Law)
 195      ! DEN         density                        (kg/m3)
 196      ! REFF        effective radius               (m)
 197      ! dyn_visc    dynamic viscosity              (kg/m/s)
 198      ! g0          gravity                        (m/s2)
 199      ! 3.0         corresponds to a growth of a factor 3 of radius with 100% RH
 200      ! 0.5         upper limit with temp correction
 201
 202      tc1(:,:,:,k) = tc(:,:,:,k)
 203      vsettl = 2.0/9.0 * g0 * den(k) * (growth_fac*reff(k))**2 / &
 204               (0.5*dyn_visc)
 205
 206      ! Determine the maximum time-step satisying the CFL condition:
 207      ! dt <= (dz)_min / v_settl
 208      ntdt=INT(dt)
 209      dtmax = dzmin / vsettl
 210      ndt_settl(k) = MAX( 1, INT( ntdt /dtmax) )
 211      ! limit maximum number of iterations
 212      IF (ndt_settl(k) > 12) ndt_settl(k) = 12
 213      dt_settl(k) = REAL(ntdt) / REAL(ndt_settl(k))
 214
 215      ! Particles radius in centimeters
 216      IF (iseas.eq.1)rcm(k) = reff(k)*100.0
 217 !srf     IF (idust.eq.1)then
 218      IF (idust.eq.1   .or. iash==1)then
 219           rwet(k) = reff(k)
 220           ratio_r(k) = 1.0
 221           rho(k) = den(k)
 222       endif
 223   END DO
 224
 225   ! Solve the bidiagonal matrix (l,l)
 226
 227 !$OMP PARALLEL DO &
 228 !$OMP DEFAULT( SHARED ) &
 229 !$OMP PRIVATE( i,   j,   l,   l2, n,   k,   rhb, rwet_priv, ratio_r, c_stokes)&
 230 !$OMP PRIVATE( free_path, c_cun, viscosity, rho_priv, vd_cor )
 231
 232   ! Loop over latitudes
 233   DO j = 1,jmx
 234
 235      DO k = 1,nmx
 236         IF (idust.eq.1 .or. iash==1) THEN
 237            rwet_priv(k) = rwet(k)
 238            rho_priv(k)  = rho(k)
 239         END IF
 240
 241         DO n = 1,ndt_settl(k)
 242
 243            ! Solve each vertical layer successively (layer l)
 244
 245            DO l = lmx,1,-1
 246               l2 = lmx - l + 1
 247
 248 !           DO j = 1,jmx
 249               DO i = 1,imx
 250
 251                  ! Dynamic viscosity
 252                  c_stokes = 1.458E-6 * tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4)
 253
 254                  ! Mean free path as a function of pressure (mb) and
 255                  ! temperature (K)
 256                  ! order of p_mid is top->sfc
 257                  free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l))
 258 !!!                 free_path = 1.1E-3/p_edge(i,j,l2)/SQRT(tmp(i,j,l))
 259
 260                  ! Slip Correction Factor
 261                  c_cun = 1.0+ free_path/rwet_priv(k)* &
 262                       (1.257 + 0.4*EXP(-1.1*rwet_priv(k)/free_path))
 263
 264                  ! Corrected dynamic viscosity (kg/m/s)
 265                  viscosity = c_stokes / c_cun
 266
 267                  ! Settling velocity
 268 !                IF (iseas.eq.1) THEN
 269 !                   rho_priv(k) = ratio_r(k)*den(k) + (1.0 - ratio_r(k))*1000.0
 270 !                END IF
 271
 272                  vd_cor = 2.0/9.0*g0*rho_priv(k)*rwet_priv(k)**2/viscosity
 273
 274                  ! Update mixing ratio
 275                  ! Order of delz is top->sfc
 276                  IF (l == lmx) THEN
 277                     tc(i,j,l,k) = tc(i,j,l,k) / &
 278                          (1.0 + dt_settl(k)*vd_cor/delz(i,j,l2))
 279                  ELSE
 280                     tc(i,j,l,k) = 1.0/(1.0+dt_settl(k)*vd_cor/delz(i,j,l2))&
 281                          *(tc(i,j,l,k) + dt_settl(k)*vd_cor /delz(i,j,l2-1) &
 282                          * tc(i,j,l+1,k))
 283                  END IF
 284               END DO   !i
 285 !           END DO   !j
 286         END DO  !l
 287
 288      END DO  !n
 289   END DO  !k
 290
 291   END DO   !j
 292 !$OMP END PARALLEL DO
 293
 294   DO n = 1,nmx
 295      DO i = 1,imx
 296         DO j = 1,jmx
 297            bstl(i,j,n) = 0.0
 298            addmass=0.
 299            DO l = 1,lmx
 300               addmass=addmass+(tc(i,j,l,n) - tc1(i,j,l,n)) * airmas(i,j,l)
 301               IF (tc(i,j,l,n) < 0.0) tc(i,j,l,n) = 1.0D-32
 302            END DO
 303            if(addmass.gt.0.)addmass=0
 304            bstl(i,j,n) = bstl(i,j,n) - addmass
 305         END DO
 306      END DO
 307   END DO
 308
 309 END SUBROUTINE vsettling
 310
 311 END MODULE MODULE_VASH_SETTLING