chem/module_gocart_settling.F

   1 MODULE MODULE_GOCART_SETTLING
   2 !
   3 ! Original GOCART
   4 ! KYKang, MKlose, CLWu various changes
   5 ! YShao [2011/09/30] complete update
   6 !
   7 ! 11/29/18 - A. Ukhov, bug fix: dust (and sea salt) mass balance was
   8 ! violated in the "settling" routine
   9
  10 CONTAINS
  11
  12 SUBROUTINE gocart_settling_driver(dt,config_flags,t_phy,moist,      &
  13            chem,rho_phy,dz8w,p8w,p_phy,                             &
  14            dustin,seasin,dx,g,                                      &
  15            dustgraset_1,dustgraset_2,dustgraset_3,                  &
  16            dustgraset_4,dustgraset_5,                               &
  17            setvel_1,setvel_2,setvel_3,setvel_4,setvel_5,imod,       &
  18            ids,ide,jds,jde,kds,kde,                                 &
  19            ims,ime,jms,jme,kms,kme,                                 &
  20            its,ite,jts,jte,kts,kte)
  21
  22   USE module_configure
  23   USE module_state_description
  24   USE module_data_gocart_dust
  25   USE module_data_gocart_seas
  26   USE module_model_constants, ONLY: mwdry
  27   IMPLICIT NONE
  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
  35   REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_moist), INTENT(IN) :: moist
  36   REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_chem), INTENT(INOUT) :: chem
  37   REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: t_phy,p_phy,dz8w,p8w,rho_phy
  38   REAL,  DIMENSION( ims:ime , jms:jme, 5 ),                          &
  39           INTENT(IN   ) ::  dustin,seasin
  40
  41   INTEGER, INTENT(IN) :: imod
  42   REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT)::                  &
  43         dustgraset_1,dustgraset_2,dustgraset_3,                      &
  44         dustgraset_4,dustgraset_5,                                   &
  45         setvel_1,setvel_2,setvel_3,setvel_4,setvel_5
  46
  47   REAL*8,  DIMENSION (1,1,5)          :: graset_dust,grasetvel_dust
  48   REAL*8,  DIMENSION (1,1,4)          :: graset_ss,grasetvel_ss
  49
  50   REAL, INTENT(IN) :: dt,dx,g
  51   INTEGER          :: kkk,nmx,i,j,k,kk,lmx,iseas,idust
  52   REAL*8, DIMENSION (1,1,kte-kts+1)   :: tmp,airden,p_mid,delz,rh
  53   REAL*8, DIMENSION (1,1,kte-kts+1,5) :: ddust
  54   REAL*8, DIMENSION (1,1,kte-kts+1,4) :: sea_salt
  55
  56   INTEGER :: uoc_flag  ! flag for UoC dust schemes
  57 !
  58 !   REAL*8, DIMENSION (5), PARAMETER :: den_dust(5)=(/2500.,2650.,2650.,2650.,2650./)
  59 !   REAL*8, DIMENSION (5), PARAMETER :: reff_dust(5)=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/)
  60 !   REAL*8, DIMENSION (4), PARAMETER :: den_seas(4)=(/2200.,2200.,2200.,2290./)
  61 !   REAL*8, DIMENSION (4), PARAMETER :: reff_seas(4)=(/0.30D-6,1.00D-6,3.25D-6,7.50D-6/)
  62 !
  63   REAL*8 conver, converi
  64   conver=1.e-9
  65   converi=1.e9
  66
  67   uoc_flag = 0
  68   if (config_flags%dust_opt .eq. 4) uoc_flag = 1
  69
  70   lmx=kte-kts+1
  71
  72   do j=jts,jte
  73     do i=its,ite
  74
  75        IF(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11) THEN
  76          do kkk=1,5
  77             ddust(1,1,kts,kkk)=dustin(i,j,kkk)*conver
  78          enddo
  79          kk=0
  80          do k=kts,kte
  81             kk=kk+1
  82
  83             p_mid(1,1,kk) =.01*p_phy(i,kte-k+kts,j)
  84             delz(1,1,kk)  =dz8w(i,kte-k+kts,j)
  85             airden(1,1,kk)=rho_phy(i,k,j)
  86             tmp(1,1,kk)= t_phy(i,k,j)
  87             rh(1,1,kk) = .95
  88             rh(1,1,kk) = MIN( .95, moist(i,k,j,p_qv) / &
  89                  (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
  90                  (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))))
  91             rh(1,1,kk) = max(1.0D-1,rh(1,1,kk))
  92          enddo
  93        ELSE
  94          kk=0
  95          DO k=kts,kte
  96             kk=kk+1
  97 !
  98             ddust(1,1,kk,1)=chem(i,k,j,p_dust_1)                         ! chem and dust in [ug/kg]
  99             ddust(1,1,kk,2)=chem(i,k,j,p_dust_2)
 100             ddust(1,1,kk,3)=chem(i,k,j,p_dust_3)
 101             ddust(1,1,kk,4)=chem(i,k,j,p_dust_4)
 102             ddust(1,1,kk,5)=chem(i,k,j,p_dust_5)
 103
 104             p_mid(1,1,kk)=.01*p_phy(i,kte-k+kts,j)
 105             delz(1,1,kk)=dz8w(i,kte-k+kts,j)                            ! delz(1) = dz8w(kte), delz(lmx)=dz8w(kts)
 106             airden(1,1,kk)=rho_phy(i,k,j)
 107             tmp(1,1,kk)   =t_phy(i,k,j)
 108             rh(1,1,kk)    = .95
 109             rh(1,1,kk)    = MIN( .95, moist(i,k,j,p_qv) /               &
 110                  (3.80*exp(17.27*(t_phy(i,k,j)-273.)/                   &
 111                  (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))))
 112             rh(1,1,kk)=max(1.0D-1,rh(1,1,kk))
 113          ENDDO
 114        ENDIF
 115        graset_dust(1,1,:)=0.
 116        graset_ss(1,1,:)=0.
 117
 118        iseas=0
 119        idust=1
 120        CALL settling(1,1,lmx,5,g,dyn_visc,ddust,tmp,p_mid,delz,       &
 121                      imod,graset_dust,grasetvel_dust, uoc_flag,       &
 122                      den_dust,reff_dust,dt,rh,idust,iseas,airden)
 123
 124        IF (config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) THEN
 125           kk=1
 126           do kkk=1,5
 127              if (kkk .le. 4) sea_salt(1,1,kts,kkk)=seasin(i,j,kkk)*conver
 128              if(ddust(1,1,kk,kkk) .ge. dustin(i,j,kkk)) ddust(1,1,kk,kkk)=dustin(i,j,kkk)
 129           enddo
 130           chem(i,kts,j,p_p25i)=chem(i,kts,j,p_p25i) &
 131                        +.25*(ddust(1,1,kk,1)+.286*ddust(1,1,kk,2))*converi
 132           chem(i,kts,j,p_p25i)=max(chem(i,kts,j,p_p25i),1.e-16)
 133           chem(i,kts,j,p_p25j)=chem(i,kts,j,p_p25j) &
 134                        +.75*(ddust(1,1,kk,1)+.286*ddust(1,1,kk,2))*converi
 135           chem(i,kts,j,p_p25j)=max(chem(i,kts,j,p_p25j),1.e-16)
 136           chem(i,kts,j,p_soila)=chem(i,kts,j,p_soila) &
 137                        +(.714*ddust(1,1,kk,2)+ddust(1,1,kk,3))*converi
 138           chem(i,kts,j,p_soila)=max(chem(i,kts,j,p_soila),1.e-16)
 139        ELSE
 140           kk = 0
 141           DO k = kts,kte
 142              kk = kk+1
 143              chem(i,k,j,p_dust_1)=ddust(1,1,kk,1)                         ! dust for size bin 1 [ug/kg]
 144              chem(i,k,j,p_dust_2)=ddust(1,1,kk,2)                         ! ...
 145              chem(i,k,j,p_dust_3)=ddust(1,1,kk,3)                         ! ...
 146              chem(i,k,j,p_dust_4)=ddust(1,1,kk,4)                         ! ...
 147              chem(i,k,j,p_dust_5)=ddust(1,1,kk,5)                         ! dust for size bin 5 (dust_opt 3: for all size bins) [ug/kg]
 148
 149              sea_salt(1,1,kk,1)=chem(i,k,j,p_seas_1)                    ! salt [ug/kg]
 150              sea_salt(1,1,kk,2)=chem(i,k,j,p_seas_2)
 151              sea_salt(1,1,kk,3)=chem(i,k,j,p_seas_3)
 152              sea_salt(1,1,kk,4)=chem(i,k,j,p_seas_4)
 153           ENDDO
 154        ENDIF
 155
 156         !graset_dust in [ug/kg][m/s]
 157        dustgraset_1(i,j)=dustgraset_1(i,j)+conver*graset_dust(1,1,1)*airden(1,1,1)*dt   !dustgraset [kg/m2]
 158        dustgraset_2(i,j)=dustgraset_2(i,j)+conver*graset_dust(1,1,2)*airden(1,1,1)*dt
 159        dustgraset_3(i,j)=dustgraset_3(i,j)+conver*graset_dust(1,1,3)*airden(1,1,1)*dt
 160        dustgraset_4(i,j)=dustgraset_4(i,j)+conver*graset_dust(1,1,4)*airden(1,1,1)*dt
 161        dustgraset_5(i,j)=dustgraset_5(i,j)+conver*graset_dust(1,1,5)*airden(1,1,1)*dt
 162
 163        setvel_1(i,j)=grasetvel_dust(1,1,1)                            ! settling velocity [m/s]
 164        setvel_2(i,j)=grasetvel_dust(1,1,2)
 165        setvel_3(i,j)=grasetvel_dust(1,1,3)
 166        setvel_4(i,j)=grasetvel_dust(1,1,4)
 167        setvel_5(i,j)=grasetvel_dust(1,1,5)
 168
 169        iseas=1
 170        idust=0
 171
 172        CALL settling(1,1,lmx,4,g,dyn_visc,sea_salt,tmp,p_mid,delz,        &
 173                      imod,graset_ss,grasetvel_ss, uoc_flag,               &
 174                      den_seas,reff_seas,dt,rh,idust,iseas,airden)
 175        IF (config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) THEN
 176           kk=1
 177           do kkk=1,4
 178              if(sea_salt(1,1,kk,kkk) .ge. seasin(i,j,kkk))sea_salt(1,1,kk,kkk)=seasin(i,j,kkk)
 179           enddo
 180           chem(i,kts,j,p_naai)=chem(i,kts,j,p_naai) &
 181                        +.25*(sea_salt(1,1,kk,1)+.942*sea_salt(1,1,kk,2))*converi
 182           chem(i,kts,j,p_naai)=max(1.e-16,chem(i,kts,j,p_naai))
 183           chem(i,kts,j,p_naaj)=chem(i,kts,j,p_naaj) &
 184                        +.75*(sea_salt(1,1,kk,1)+.942*sea_salt(1,1,kk,2))*converi
 185           chem(i,kts,j,p_naaj)=max(1.e-16,chem(i,kts,j,p_naaj))
 186           chem(i,kts,j,p_seas)=chem(i,kts,j,p_seas) &
 187                        +(.058*sea_salt(1,1,kk,2)+sea_salt(1,1,kk,3))*converi
 188           chem(i,kts,j,p_seas)=max(1.e-16,chem(i,kts,j,p_seas))
 189        ELSE
 190           kk=0
 191           DO k=kts,kte
 192              kk=kk+1
 193              chem(i,k,j,p_seas_1)=sea_salt(1,1,kk,1)
 194              chem(i,k,j,p_seas_2)=sea_salt(1,1,kk,2)
 195              chem(i,k,j,p_seas_3)=sea_salt(1,1,kk,3)
 196              chem(i,k,j,p_seas_4)=sea_salt(1,1,kk,4)
 197           ENDDO
 198        ENDIF
 199
 200     enddo  ! i
 201   enddo  ! j
 202
 203  END SUBROUTINE gocart_settling_driver
 204
 205
 206  subroutine settling(imx,jmx,lmx,nmx,g0,dyn_visc,tc,tmp,p_mid,delz, &
 207                     imod,graset, grasetvel, uoc,                    &
 208                     den_in,reff_in,dt,rh,idust,iseas,airden)
 209
 210 ! ****************************************************************************
 211 ! *                                                                          *
 212 ! *  Calculate the loss by settling, using an implicit method                *
 213 ! *                                                                          *
 214 ! *  Input variables:                                                        *
 215 ! *    SIGE(k)         - sigma coordinate of the vertical edges              *
 216 ! *    PS(i,j)         - Surface pressure (mb)                               *
 217 ! *    TMP(i,j,k)      - Air temperature  (K)                                *
 218 ! *    CT(i,j)         - Surface exchange coeff for moisture
 219 ! *                                                                          *
 220 ! ****************************************************************************
 221
 222   IMPLICIT  NONE
 223
 224   INTEGER, INTENT(IN) :: imx, jmx, lmx, nmx,iseas,idust
 225   INTEGER             :: ntdt
 226   REAL,    INTENT(IN) :: dt,g0,dyn_visc
 227   REAL*8,  INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx),  &
 228                          rh(imx,jmx,lmx), p_mid(imx,jmx,lmx),airden(imx,jmx,lmx)
 229 !
 230   REAL*8, INTENT(IN)  :: den_in(nmx), reff_in(nmx)
 231
 232   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 233
 234   INTEGER, INTENT(IN)   :: imod, uoc
 235   REAL*8, INTENT(INOUT) :: graset(imx,jmx,nmx)
 236   REAL*8, INTENT(OUT)   :: grasetvel(imx,jmx,nmx)
 237
 238   REAL*8    :: den(nmx), reff(nmx)                    ! local variables here
 239   REAL*8    :: dt_settl(nmx), rcm(nmx), rho(nmx)
 240   INTEGER   :: ndt_settl(nmx)
 241   REAL*8    :: dzmin, vsettl, dtmax, pres, rhb, rwet(nmx), ratio_r(nmx)
 242   REAL*8    :: c_stokes, free_path, c_cun, viscosity, growth_fac
 243   REAL*8    :: vd_cor(lmx), vd_wk1
 244   INTEGER   :: k, n, i, j, l, l2
 245
 246   REAL*8    :: transfer_to_below_level,temp_tc
 247
 248 ! for sea-salt:
 249   REAL*8, PARAMETER :: c1=0.7674, c2=3.079, c3=2.573E-11, c4=-1.424
 250
 251 ! for OMP:
 252   REAL*8 :: rwet_priv(nmx), rho_priv(nmx)
 253
 254 ! IF (type) /= 'dust' .AND. TRIM(aero_type) /= 'sea_salt') RETURN
 255   IF ( idust.ne.1 .and. iseas.ne.1 ) RETURN
 256   WHERE ( tc(:,:,:,:) < 0.0 ) tc(:,:,:,:) = 1.0D-32
 257
 258   den = den_in
 259   reff = reff_in
 260   dzmin = MINVAL(delz(:,:,:))
 261   IF (idust == 1) growth_fac = 1.0
 262   IF (iseas == 1) growth_fac = 3.0
 263 !
 264   IF (idust == 1 .and. uoc == 1) then
 265       den(1) = 2650.            ! constant density is use in UoC dust emission schemes;
 266 !                               dust radii are now consistent with GOCART dust schemes.  [mklose, 03082015]
 267   ENDIF
 268
 269   DO k = 1, nmx                                                   ! k for different size bins
 270 !
 271 ! Settling velocity (m/s) for each tracer (Stokes Law)
 272 ! DEN         density                        (kg/m3)
 273 ! REFF        effective radius               (m)
 274 ! dyn_visc    dynamic viscosity              (kg/m/s)
 275 ! g0          gravity                        (m/s2)
 276 ! 3.0         corresponds to a growth of a factor 3 of radius with 100% RH
 277 ! 0.5         upper limit with temp correction
 278 !
 279      vsettl = 4.0/9.0 * g0 * den(k) * (growth_fac*reff(k))**2 / dyn_visc
 280 !
 281 ! Determine the maximum time-step satisying CFL: dt <= (dz)_min / v_settl
 282 !
 283      ntdt = INT(dt)
 284      dtmax = dzmin / vsettl
 285      ndt_settl(k) = MAX( 1,INT(ntdt/dtmax) )
 286
 287 ! Limit maximum number of iterations
 288      IF (ndt_settl(k) > 12) ndt_settl(k) = 12
 289      dt_settl(k) = REAL(ntdt) / REAL(ndt_settl(k))
 290
 291 ! Particles radius in centimeters
 292      IF (iseas.eq.1)rcm(k) = reff(k)*100.0
 293      IF (idust.eq.1) THEN
 294       rwet(k) = reff(k)
 295       ratio_r(k) = 1.0
 296       rho(k) = den(k)
 297      ENDIF
 298   ENDDO
 299
 300 ! Solve the bidiagonal matrix (l,l)
 301
 302 !$OMP PARALLEL DO &
 303 !$OMP DEFAULT( SHARED ) &
 304 !$OMP PRIVATE( i,   j,   l,   l2, n,   k,   rhb, rwet_priv, ratio_r, c_stokes)&
 305 !$OMP PRIVATE( free_path, c_cun, viscosity, rho_priv, vd_cor )
 306
 307 ! Loop over latitudes
 308
 309   DO j = 1,jmx                      ! lat loop
 310     DO i = 1,imx                    ! lon loop
 311       DO k = 1,nmx                  ! bin loop
 312         graset(i,j,k)=0.
 313         grasetvel(i,j,k)=0.
 314
 315         IF (idust.eq.1) THEN
 316            rwet_priv(k) = rwet(k)
 317            rho_priv(k)  = rho(k)
 318         END IF
 319
 320         DO n = 1,ndt_settl(k)        ! time loop
 321
 322           transfer_to_below_level=0
 323
 324           DO l = lmx,1,-1            ! height loop, from top
 325              l2 = lmx - l + 1
 326
 327              IF (iseas.eq.1) THEN
 328                 rhb = MIN(9.9D-1, rh(i,j,l))                             ! Aerosol growth with relative humidity (Gerber, 1985)
 329                 rwet_priv(k) = 0.01*(c1*rcm(k)**c2/(c3*rcm(k)**c4 -  &
 330                 LOG10(rhb)) + rcm(k)**3)**0.33                           ! td changed to LOG10
 331                 ratio_r(k) = (reff(k)/rwet_priv(k))**3.0
 332              END IF
 333
 334              c_stokes = 1.458E-6*tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4)    ! Dynamic viscosity
 335              free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l))           ! Free path as func of pres(mb) and temp(K); order of p_mid: top->sfc
 336
 337              c_cun = 1.0+free_path/rwet_priv(k)*                     &
 338              (1.257 + 0.4*EXP(-1.1*rwet_priv(k)/free_path))              ! Slip correction
 339              viscosity = c_stokes / c_cun                                ! Corrected dynamic viscosity (kg/m/s)
 340
 341              IF (iseas.eq.1) THEN
 342                 rho_priv(k) = ratio_r(k)*den(k) + (1.-ratio_r(k))*1000.
 343              END IF
 344
 345              vd_cor(l) = 2./9.*g0*rho_priv(k)*rwet_priv(k)**2/viscosity  ! Settling velocity, depends on temp
 346
 347              ! Update mixing ratio; order of delz: top->sfc
 348              temp_tc=tc(i,j,l,k)                                          ! temp_tc - for temporal storage [ug/kg]
 349              vd_wk1 = dt_settl(k)*vd_cor(l)/delz(i,j,l2)                  ! fraction to leave level
 350
 351              tc(i,j,l,k)   =  tc(i,j,l,k)*(1.- vd_wk1)+transfer_to_below_level          ! [ug/kg]
 352
 353              IF (l==1) THEN
 354                 graset(i,j,k)=graset(i,j,k)+vd_cor(l)*temp_tc/ndt_settl(k) ! [ug/kg][m/s]
 355                 grasetvel(i,j,k)=vd_cor(l)                                 ! [m/s]
 356              ELSE
 357                 transfer_to_below_level = (temp_tc*vd_wk1)*((delz(i,j,l2)*airden(i,j,l))/(delz(i,j,l2+1)*airden(i,j,l-1))) ! [ug/kg]
 358              ENDIF
 359           ENDDO                 !l, height
 360         ENDDO                   !n, time
 361       ENDDO                     !k, bin
 362 !
 363     ENDDO                     !i
 364   END DO                      !j
 365 !$OMP END PARALLEL DO
 366
 367 END SUBROUTINE settling
 368
 369 END MODULE MODULE_GOCART_SETTLING