chem/module_vertmx_wrf.F

   1 MODULE module_vertmx_wrf
   2
   3 CONTAINS
   4
   5 !-----------------------------------------------------------------------
   6     SUBROUTINE vertmx( dt, phi, kt_turb, dryrho,  &
   7        zsigma, zsigma_half, vd, kts, ktem1 )
   8 !  !! purpose - calculate change in time of phi due to vertical mixing
   9 !  !!     and dry deposition (for 1 species, 1 vertical column, 1 time step)
  10 !  !! Mariusz Pagowski, March 2001
  11 !  !! conventions used:
  12 !  !! input is lower case
  13 !  !! output is upper case
  14 !
  15 !  !! modifications by R Easter, May 2006
  16 !  !!    added dryrho so this routine conserves column mass burde
  17 !  !!       when dry deposition velocity is zero
  18 !  !!    changed "kte" to "ktem1" for consistency with the kte in WRF
  19 !
  20 ! ARGUMENTS
  21 !
  22 ! dt = time step (s)
  23 ! phi = initial/final (at input/output) species mixing ratios at "T" points
  24 ! kt_turb = turbulent exchange coefficients (m^2/s) at "W" points
  25 ! dryrho = dry air density (kg/m^3) at "T" points
  26 ! zsigma = heights (m) at "W" points
  27 ! zsigma_half = heights (m) at "T" points
  28 ! vd = dry deposition velocity (m/s)
  29 ! kts, ktem1 = vertical indices of bottom and top "T" points
  30 !
  31       IMPLICIT NONE
  32
  33 ! .. Scalar Arguments ..
  34       INTEGER, INTENT(IN) :: kts,ktem1
  35       REAL, INTENT(IN) :: dt, vd
  36 ! ..
  37 ! .. Array Arguments ..
  38       REAL, INTENT(IN), DIMENSION (kts:ktem1+1) :: kt_turb, zsigma
  39       REAL, INTENT(IN), DIMENSION (kts:ktem1) :: dryrho, zsigma_half
  40       REAL, INTENT(INOUT), DIMENSION (kts:ktem1) :: phi
  41 ! ..
  42 ! .. Local Scalars ..
  43       INTEGER :: k
  44 ! ..
  45 ! .. Local Arrays ..
  46       REAL, DIMENSION (kts+1:ktem1) :: a_coeff
  47       REAL, DIMENSION (kts:ktem1) :: b_coeff, lhs1, lhs2, lhs3, rhs
  48 ! ..
  49 ! .. External Subroutines ..
  50 !     EXTERNAL coeffs, rlhside, tridiag
  51 ! ..
  52       CALL coeffs( kts, ktem1, dryrho, zsigma, zsigma_half, a_coeff, b_coeff )
  53
  54       CALL rlhside( kts, ktem1, kt_turb, dryrho, a_coeff, b_coeff, &
  55         phi, dt, vd, rhs, lhs1, lhs2, lhs3 )
  56
  57       CALL tridiag( kts, ktem1, lhs1, lhs2, lhs3, rhs )
  58
  59       DO k = kts,ktem1
  60         phi(k) = rhs(k)
  61       END DO
  62
  63     END SUBROUTINE vertmx
  64
  65
  66 !-----------------------------------------------------------------------
  67     SUBROUTINE rlhside( kts, ktem1, k_turb, dryrho, a_coeff, b_coeff, &
  68       phi, dt, vd, rhs, lhs1, lhs2, lhs3 )
  69   !! to calculate right and left hand sides in diffusion equation
  70   !! for the tridiagonal solver
  71   !! Mariusz Pagowski, March 2001
  72   !! conventions used:
  73   !! input is lower case
  74   !! output is upper case
  75       IMPLICIT NONE
  76
  77 ! .. Scalar Arguments ..
  78       INTEGER, INTENT(IN) :: kts,ktem1
  79       REAL, INTENT(IN) :: dt, vd
  80 ! ..
  81 ! .. Array Arguments ..
  82       REAL, INTENT(IN), DIMENSION (kts:ktem1+1) :: k_turb
  83       REAL, INTENT(IN), DIMENSION (kts+1:ktem1) :: a_coeff
  84       REAL, INTENT(IN), DIMENSION (kts:ktem1) :: b_coeff, dryrho, phi
  85       REAL, INTENT(OUT), DIMENSION (kts:ktem1) :: lhs1, lhs2, lhs3, rhs
  86 ! ..
  87 ! .. Local Scalars ..
  88       REAL :: a1, a2, alfa_explicit = .25, beta_implicit = .75
  89       INTEGER :: i
  90
  91 ! ..
  92       i = kts
  93       a2 = a_coeff(i+1)*k_turb(i+1)
  94       rhs(i) = (1./(dt*b_coeff(i)) - alfa_explicit*(vd*dryrho(i)+a2))*phi(i) + &
  95         alfa_explicit*(a2*phi(i+1))
  96       lhs1(i) = 0.
  97       lhs2(i) = 1./(dt*b_coeff(i)) + beta_implicit*(vd*dryrho(i)+a2)
  98       lhs3(i) = -beta_implicit*a2
  99
 100       DO i = kts+1, ktem1-1
 101         a1 = a_coeff(i)*k_turb(i)
 102         a2 = a_coeff(i+1)*k_turb(i+1)
 103
 104         rhs(i) = (1./(dt*b_coeff(i)) - alfa_explicit*(a1+a2))*phi(i) + &
 105           alfa_explicit*(a1*phi(i-1) + a2*phi(i+1))
 106
 107         lhs1(i) = -beta_implicit*a1
 108         lhs2(i) = 1./(dt*b_coeff(i)) + beta_implicit*(a1+a2)
 109         lhs3(i) = -beta_implicit*a2
 110       END DO
 111
 112       i = ktem1
 113       a1 = a_coeff(i)*k_turb(i)
 114       rhs(i) = (1./(dt*b_coeff(i)) - alfa_explicit*(a1   ))*phi(i) + &
 115         alfa_explicit*(a1*phi(i-1))
 116       lhs1(i) = -beta_implicit*a1
 117       lhs2(i) = 1./(dt*b_coeff(i)) + beta_implicit*(a1   )
 118       lhs3(i) = 0.
 119
 120     END SUBROUTINE rlhside
 121
 122
 123 !-----------------------------------------------------------------------
 124     SUBROUTINE tridiag( kts, ktem1, a, b, c, f )
 125   !! to solve system of linear eqs on tridiagonal matrix n times n
 126   !! after Peaceman and Rachford, 1955
 127   !! a,b,c,F - are vectors of order n
 128   !! a,b,c - are coefficients on the LHS
 129   !! F - is initially RHS on the output becomes a solution vector
 130   !! Mariusz Pagowski, March 2001
 131   !! conventions used:
 132   !! input is lower case
 133   !! output is upper case
 134       IMPLICIT NONE
 135
 136 ! .. Scalar Arguments ..
 137       INTEGER, INTENT(IN) :: kts,ktem1
 138 ! ..
 139 ! .. Array Arguments ..
 140       REAL, INTENT(IN), DIMENSION (kts:ktem1) :: a, b, c
 141       REAL, INTENT(INOUT), DIMENSION (kts:ktem1) :: f
 142 ! ..
 143 ! .. Local Scalars ..
 144       REAL :: p
 145       INTEGER :: i
 146 ! ..
 147 ! .. Local Arrays ..
 148       REAL, DIMENSION (kts:ktem1) :: q
 149 ! ..
 150       q(kts) = -c(kts)/b(kts)
 151       f(kts) = f(kts)/b(kts)
 152
 153       DO i = kts+1, ktem1
 154         p = 1./(b(i)+a(i)*q(i-1))
 155         q(i) = -c(i)*p
 156         f(i) = (f(i)-a(i)*f(i-1))*p
 157       END DO
 158
 159       DO i = ktem1 - 1, kts, -1
 160         f(i) = f(i) + q(i)*f(i+1)
 161       END DO
 162
 163     END SUBROUTINE tridiag
 164
 165
 166 !-----------------------------------------------------------------------
 167     SUBROUTINE coeffs( kts, ktem1, dryrho, &
 168       z_sigma, z_sigma_half, a_coeff, b_coeff )
 169 ! !! to calculate coefficients in diffusion equation
 170 ! !! Mariusz Pagowski, March 2001
 171 ! !! conventions used:
 172 ! !! input is lower case
 173 ! !! output is upper case
 174 ! .. Scalar Arguments ..
 175       IMPLICIT NONE
 176
 177       INTEGER, INTENT(IN) :: kts,ktem1
 178 ! ..
 179 ! .. Array Arguments ..
 180       REAL, INTENT(IN), DIMENSION (kts:ktem1+1) :: z_sigma
 181       REAL, INTENT(IN), DIMENSION (kts:ktem1) :: z_sigma_half, dryrho
 182       REAL, INTENT(OUT), DIMENSION (kts+1:ktem1) :: a_coeff
 183       REAL, INTENT(OUT), DIMENSION (kts:ktem1) :: b_coeff
 184 ! ..
 185 ! .. Local Scalars ..
 186       INTEGER :: i
 187       REAL :: dryrho_at_w
 188 ! ..
 189       DO i = kts, ktem1
 190         b_coeff(i) = 1./(dryrho(i)*(z_sigma(i+1)-z_sigma(i)))
 191       END DO
 192
 193       DO i = kts+1, ktem1
 194         dryrho_at_w = 0.5*(dryrho(i)+dryrho(i-1))
 195         a_coeff(i) = dryrho_at_w/(z_sigma_half(i)-z_sigma_half(i-1))
 196       END DO
 197
 198     END SUBROUTINE coeffs
 199
 200 !-----------------------------------------------------------------------
 201 END MODULE module_vertmx_wrf