wrfv2_fire/chem/module_hetn2o5.F

   1 MODULE module_hetn2o5
   2
   3
   4 !temporary:
   5 ! from dentener and crutzen for 5S, 160E on 100mb levels startin at 1000mb
   6     REAL, PARAMETER, DIMENSION(10) :: &
   7         kheti = (/4.2e-5, 3.8e-5, 1.e-5,3.e-6,4.e-6, 2.e-6,3.e-6,5.e-6,1.3e-5,1.8e-5/)
   8
   9     !kheta:  on aerosol
  10    REAL, SAVE, DIMENSION(120) :: kheta
  11
  12    LOGICAL :: do_aerosol=.TRUE.
  13
  14
  15 ! sticking coefficients for cloud water and cloud ice
  16
  17    REAL    , PARAMETER, PRIVATE :: gammacldw = 0.05,   &
  18                                    gammacldi = 0.03  !cms check !
  19
  20    REAL     , PARAMETER, PRIVATE ::     rhograul = 400. ,&
  21                                         rhoice   = 917.
  22
  23    REAL, PARAMETER, PRIVATE :: pi = 3.1415926
  24
  25
  26     ! D_g: diffusivity of species in gas phase in m^2/s
  27        REAL, PARAMETER, PRIVATE ::  D_g = 0.1E-6
  28                              !in Match this is 1.e-5 (Schwartz ?!)
  29
  30
  31
  32
  33     ! abar_c:  mass mean radius for cloud drops in m
  34        REAL, PARAMETER, PRIVATE :: abar_c = 10.E-6
  35
  36
  37
  38     ! RSTAR2 : universal gas constant in J/(kmol K)
  39        REAL, PARAMETER, PRIVATE  :: RSTAR2 = 8314.
  40
  41
  42     ! parameters from Lin scheme
  43        REAL    , PARAMETER, PRIVATE ::     xnor = 8.0e6
  44        REAL    , PARAMETER, PRIVATE ::     xnos = 3.0e6
  45        REAL    , PARAMETER, PRIVATE ::     xnog = 4.0e6
  46
  47
  48
  49  CONTAINS
  50   SUBROUTINE hetn2o5calc(hetn2o5, rho, T, QC, QR, QI, QS, QG, &
  51                        rhowater, rhosnow, M_n2o5,             &
  52                        P_QI,P_QS,P_QG,                        &
  53                        ids,ide, jds,jde, kds,kde,             &
  54                        ims,ime, jms,jme, kms,kme,             &
  55                        its,ite, jts,jte, kts,kte              )
  56
  57
  58
  59       IMPLICIT NONE
  60
  61
  62
  63   INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde , &
  64                                       ims,ime, jms,jme, kms,kme , &
  65                                       its,ite, jts,jte, kts,kte , &
  66                                       P_QI,P_QS,P_QG
  67
  68
  69    REAL, DIMENSION( its:ite , kts:kte , jts:jte),        &
  70          INTENT(INOUT )                               :: hetn2o5
  71
  72
  73
  74   REAL, INTENT(IN   ) ::                                rhosnow,  &
  75                                                         rhowater, &
  76                                                         M_n2o5
  77
  78
  79   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                 &
  80         INTENT(IN   ) ::                                          &
  81                                                               QC, &
  82                                                               QR, &
  83                                                               QI, &
  84                                                               QS, &
  85                                                               QG, &
  86                                                              rho, &
  87                                                                T
  88
  89
  90 !local vars
  91
  92     REAL :: orho, tmp1
  93
  94
  95   ! mass mean radii of the different hydrometeors
  96     REAL ::                &
  97                                  abar_r, abar_i, abar_s, abar_g
  98
  99
 100   ! nu_th: thermal velocity of species in m/s
 101     REAL :: nu_th
 102
 103     ! tau_Dg: gas diffusion timescale
 104     ! tau_i : timescale for mass transfer across interface of hydrometeor
 105        REAL :: tau_i, tau_Dg
 106
 107
 108     !  L : water contents in cloud drops, rain drops, hail stones,...
 109     !  in (cm^3 H_2O / cm^3 air)
 110        REAL :: L
 111
 112
 113     ! loss rates:
 114        REAL :: kc, kr, ki, ks, kg
 115
 116
 117       INTEGER :: i,j,k
 118
 119
 120 !----
 121
 122    j_loop:  DO j = jts, jte
 123    i_loop:  DO i = its, ite
 124
 125      DO k = kts, kte
 126            orho = 1./rho(i,k,j)
 127
 128
 129
 130          kc = 0.
 131          kr = 0.
 132          ki = 0.
 133          ks = 0.
 134          kg = 0.
 135
 136
 137
 138             !! abar_c
 139
 140
 141           IF(QR(i,k,j) .GT. 1.e-12) THEN
 142             tmp1=sqrt(pi*rhowater*xnor*orho/QR(i,k,j))
 143             !!xlambdar1Dr(k)=sqrt(tmp1)
 144             abar_r =  2./MAX(sqrt(tmp1), 1.E-8 ) !if abar is large, kt becomes small
 145           ENDIF
 146
 147
 148 !!$           IF(QI(i,k,j) .GT. 1.e-8) THEN
 149 !!$            tmp1=sqrt(pi*rhoice*xnoi*orho/QI(i,k,j))
 150 !!$            !!xlambdar1Di(k)=sqrt(tmp1)
 151 !!$            abar_i = 2./MAX(sqrt(tmp1), 1.E-8 )
 152 !!$           ENDIF
 153
 154
 155              abar_i =  abar_c
 156
 157            IF(QS(i,k,j) .GT. 1.e-12) THEN
 158             tmp1=sqrt(pi*rhosnow*xnos*orho/QS(i,k,j))
 159             !!xlambdar1Ds(k)=sqrt(tmp1)
 160              abar_s =  2./MAX(sqrt(tmp1), 1.E-8 )
 161            ENDIF
 162
 163
 164
 165            IF(QG(i,k,j) .GT. 1.e-12) THEN
 166             tmp1=sqrt(pi*rhograul*xnog*orho/QG(i,k,j))
 167             !!xlambdar1Dg(k)=sqrt(tmp1)
 168             abar_g = 2./MAX(sqrt(tmp1), 1.E-8 )
 169            ENDIF
 170
 171
 172
 173
 174 ! calculate thermal velocity of species
 175         nu_th = SQRT(8*RSTAR2*T(i,k,j)/(pi*M_n2o5))
 176
 177
 178
 179 ! calculate timescales
 180   !cloud droplets
 181     IF(QC(i,k,j) .GT. 1.e-12) THEN
 182       tau_i = (4.*abar_c) / (3.* nu_th * gammacldw)
 183       tau_Dg = (abar_c**2) / (3.* D_g)
 184
 185       L  = (rho(i,k,j) / rhowater) * QC(i,k,j)
 186       kc = L/(tau_i +  tau_Dg)
 187     ENDIF
 188
 189   !rain
 190     IF(QR(i,k,j) .GT. 1.e-12) THEN
 191       tau_i = (4.*abar_r) / (3.* nu_th * gammacldw)
 192       tau_Dg = (abar_r**2) / (3.* D_g)
 193
 194       L  = (rho(i,k,j) / rhowater) * QR(i,k,j)
 195       kr = L/(tau_i +  tau_Dg)
 196     ENDIF
 197
 198   !cloud ice
 199     IF(QI(i,k,j) .GT. 1.e-12) THEN
 200       tau_i = (4.*abar_i) / (3.* nu_th * gammacldi)
 201       tau_Dg = (abar_i**2) / (3.* D_g)
 202
 203       L  = (rho(i,k,j) / rhowater) * QI(i,k,j)
 204       ki = L/(tau_i +  tau_Dg)
 205     ENDIF
 206
 207   !snow
 208      IF(QS(i,k,j).GT. 1.e-12) THEN
 209       tau_i = (4.*abar_s) / (3.* nu_th * gammacldi)
 210       tau_Dg = (abar_s**2) / (3.* D_g)
 211
 212       L  = (rho(i,k,j) / rhowater) * QS(i,k,j)
 213       ks = L/(tau_i +  tau_Dg)
 214      ENDIF
 215
 216   !graupel
 217      IF(QG(i,k,j).GT. 1.e-12) THEN
 218       tau_i = (4.*abar_g) / (3.* nu_th * gammacldi)
 219       tau_Dg = (abar_g**2) / (3.* D_g)
 220
 221       L  = (rho(i,k,j) / rhowater) * QG(i,k,j)
 222       kg = L/(tau_i +  tau_Dg)
 223      ENDIF
 224
 225
 226
 227 !!HERE THE VENTILATION COEFF SHOULD BE INCLUDED AND IT PROBABLY DOES NOT HAPPEN ON ICE!!
 228
 229 !!    hetn2o5(i,k,j) = kc + kr + ki + ks + kg
 230
 231
 232
 233
 234
 235      hetn2o5(i,k,j) =  kc + kr + kheta(k)  !!+ kg
 236      !hetn2o5(i,k,j) = 0. !n2o5_test
 237
 238
 239      ENDDO
 240
 241    ENDDO i_loop
 242    ENDDO j_loop
 243
 244   END SUBROUTINE hetn2o5calc
 245
 246
 247
 248
 249
 250 !-----------------------------------------------------------------
 251
 252  SUBROUTINE hetn2o5_ini(      pb, pp, z,                     &
 253                                ids, ide, jds, jde, kds, kde, &
 254                                ims, ime, jms, jme, kms, kme, &
 255                                its, ite, jts, jte, kts, kte)
 256
 257 #ifdef DM_PARALLEL
 258    USE module_dm
 259 #endif
 260
 261     IMPLICIT NONE
 262
 263
 264     INTEGER, INTENT(IN    )  :: ids, ide, jds, jde, kds, kde, &
 265                                 ims, ime, jms, jme, kms, kme, &
 266                                 its, ite, jts, jte, kts, kte
 267
 268
 269     REAL, DIMENSION( ims:ime, kms:kme, jms:jme),  &
 270               INTENT(IN    )  ::                      z, pb, pp
 271
 272      REAL, DIMENSION(10) :: pin !moguntia
 273
 274
 275   LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
 276
 277
 278 !local
 279
 280      REAL, DIMENSION(kms:kme) :: zloc, ploc
 281
 282       INTEGER :: k, kk, l_low, l_up
 283       REAL :: m, b, dp
 284
 285
 286 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * RWORDSIZE )
 287   IF (.NOT. do_aerosol) RETURN
 288
 289 !currently not really necessary
 290 dm_on_monitor: IF ( wrf_dm_on_monitor() ) THEN
 291
 292   kheta(:) = 0.
 293
 294
 295
 296    pin(kms)=1000.
 297  DO k=2,10
 298    pin(k)=pin(k-1)-100.
 299  ENDDO
 300
 301
 302  DO k=kms,kme
 303
 304   zloc(k) = z(its,k,jts)
 305   ploc(k) = pb(its,k,jts) + pp(its,k,jts)
 306
 307    IF ( zloc(k) .GT. 1.e6 .OR. zloc(k) .LT. 0. ) THEN
 308       CALL wrf_error_fatal ("STOP: hetn2o5_ini")
 309    ENDIF
 310
 311  ENDDO
 312
 313  DO k=kms,kme
 314    IF ( ploc(k) .GT. 1.e5 ) THEN
 315       kheta(k) = kheti(1)
 316    ELSE
 317
 318     l_low = 11.-CEILING(ploc(k)/1.e4)
 319     l_up=l_low+1
 320
 321
 322      dp=ploc(k)/100. - pin(l_up)
 323      m=(kheti(l_low)-kheti(l_up))/100.
 324      kheta(k)=kheti(l_up) +  m*dp
 325
 326
 327
 328    !print *,l_low
 329    !print *, ploc(k), pin(l_low), pin(l_up)
 330    !print *,  "  ", kheti(l_low),kheti(l_up), kheta(k)
 331
 332
 333
 334    END IF
 335
 336
 337  ENDDO
 338
 339
 340
 341
 342 ENDIF dm_on_monitor
 343
 344   DM_BCAST_MACRO( kheta )
 345
 346
 347
 348  END SUBROUTINE hetn2o5_ini
 349
 350
 351 END MODULE module_hetn2o5
 352
 353
 354
 355
 356
 357