da_get_trh.inc

1 subroutine da_get_trh( input_file, dim1, dim2, dim3, k, temp, rh )
2 
3    !---------------------------------------------------------------------------
4    ! Purpose: Calculates T, RH from input WRF file.
5    !---------------------------------------------------------------------------
6 
7    implicit none
8 
9    character(len=200), intent(in) :: input_file       ! NETCDF file nane.
10    integer, intent(in)   :: dim1, dim2, dim3          ! Dimensions.
11    integer, intent(in)   :: k                         ! Model level.
12    real, intent(out)     :: temp(1:dim1,1:dim2)       ! Temperature.
13    real, intent(out)     :: rh(1:dim1,1:dim2)         ! Relative humidity.
14 
15    character(len=10)    :: var                       ! Variable to search for. var = "T"
16    integer               :: i, j                      ! Loop counters.
17 
18    real                  :: thetap(1:dim1,1:dim2)     ! Perturbation potential temperature.
19    real                  :: pb(1:dim1,1:dim2)         ! Base state pressure.
20    real                  :: pp(1:dim1,1:dim2)         ! Pressure perturbation.
21    real                  :: x(1:dim1,1:dim2)          ! Vapor mixing ratio.
22 
23    real                  :: theta                     ! Potential temperature.
24    real                  :: p                         ! Pressure.
25    real                  :: q                         ! Specific humidity.
26    real                  :: t_c                       ! Temp(Celsius).
27    real                  :: es                        ! Saturation vapor pressure.
28    real                  :: qs                        ! Saturation specific humidity.
29 
30    var = "T" ! Perturbation potential temperature in WRF.
31    call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, thetap)
32 
33    var = "PB"  ! Base state pressure in WRF.
34    call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, pb)
35 
36    var = "P" ! Perturbation pressure in WRF.
37    call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, pp)
38 
39    var = "QVAPOR"  ! Water vapor mixing ratio.
40    call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, x)
41 
42    do j = 1, dim2
43       do i = 1, dim1
44 
45          ! Convert p', theta' to T:
46          theta = ts0 + thetap(i,j)                 ! Theta = Theta0 + Thetap
47          p = pb(i,j) + pp(i,j)                     ! p = base p + pert p.
48          temp(i,j) = theta *( p * ps0_inv)**kappa ! Theta to T.
49 
50          ! Convert to specific humidity.
51          q = x(i,j) /( 1.0 + x(i,j))
52 
53          ! Calculate relative humidity:
54          t_c = temp(i,j) - t_kelvin
55          es = es_alpha * exp( es_beta * t_c /( t_c + es_gamma))
56          qs = rd_over_rv * es /( p - rd_over_rv1 * es)
57          rh(i,j) = 100.0 * q / qs
58       end do
59    end do
60 
61 end subroutine da_get_trh
62 
63