da_tp_to_qs1.inc

1 subroutine da_tp_to_qs1 (grid, es, qs)
2 
3    !---------------------------------------------------------------------------
4    !  Purpose: Convert T/p to saturation specific humidity.
5    !
6    !  Method: qs = es_alpha * es / ( p - ( 1 - rd_over_rv ) * es ).
7    !          use Rogers & Yau (1989) formula: es = a exp( bTc / (T_c + c) ).
8    !
9    !  This da_tp_to_qs1 was added and called by the corrected subroutine
10    !       da_tpq_to_rh_lin.
11    !---------------------------------------------------------------------------
12 
13    implicit none
14 
15    type (domain), intent(in)  :: grid
16    real,          intent(out) :: es(its:ite,jts:jte,kts:kte)
17    real,          intent(out) :: qs(its:ite,jts:jte,kts:kte)
18 
19    integer :: i, j, k      ! Loop counters.
20    real    :: t_c          ! Working variable.
21 
22    if (trace_use_dull) call da_trace_entry("da_tp_to_qs1")
23 
24    !---------------------------------------------------------------------------
25    ! [1.0] initialise:
26    !---------------------------------------------------------------------------
27     
28 
29    do k = kts, kte
30       do j = jts, jte
31          do i = its, ite
32 
33             !------------------------------------------------------------------
34             ! [1.0] initialise:
35             !------------------------------------------------------------------
36 
37             t_c = grid%xb % t(i,j,k) - t_kelvin
38    
39             !------------------------------------------------------------------
40             ! [2.0] Calculate saturation vapour pressure:
41             !------------------------------------------------------------------
42 
43             es(i,j,k) = es_alpha * exp( es_beta * t_c / ( t_c + es_gamma ) )
44    
45             !------------------------------------------------------------------
46             ! [3.0] Calculate saturation specific humidity:
47             !------------------------------------------------------------------
48 
49             qs(i,j,k) = rd_over_rv * es(i,j,k) / &
50                      (grid%xb % p(i,j,k) - rd_over_rv1 * es(i,j,k))
51 
52          end do
53       end do
54    end do
55 
56    if (trace_use_dull) call da_trace_exit("da_tp_to_qs1")
57 
58 end subroutine da_tp_to_qs1
59 
60