da_tpq_to_slp.inc

1 subroutine da_tpq_to_slp ( t, q, p, terr, psfc, slp)
2 
3    !-----------------------------------------------------------------------
4    ! purpose:  computes sea level pressure from the rule                
5    !              t1/t2=(p1/p2)**(gamma*r/g).                              
6    !                                                                       
7    !     input       t        temperature
8    !                 q        mixing ratio
9    !                 p        pressure
10    !                 terr     terrain
11    !                 psfc     surface pressure
12    !                                                                       
13    !     output      slp      sea level pressure    
14    !-----------------------------------------------------------------------        
15 
16    implicit none
17 
18    real, intent(in)    :: terr, psfc
19    real, intent(in)    :: t(kms:kme)
20    real, intent(in)    :: q(kms:kme)
21    real, intent(in)    :: p(kms:kme)
22    real, intent(inout) :: slp
23 
24    integer         :: k, klo, khi
25    real            :: pl, t0, ts, xterm,tlo, thi, tl
26                                           
27    real, parameter :: gamma  = 6.5e-3
28    real, parameter :: tc     = t_kelvin+17.5
29    real, parameter :: pconst = 10000.0
30    real, parameter :: eps    = 0.622
31 
32    if (trace_use) call da_trace_entry("da_tpq_to_slp")
33                                                                        
34    ! sea level pressure                                            
35                                                                          
36    xterm=gamma* gas_constant / gravity                                                   
37                                                                        
38    ! compute pressure at pconst mb above surface (pl)              
39                                                                         
40    if (terr <= 0.0) then
41       slp = psfc
42       return
43    end if
44 
45    pl  = psfc - pconst                                        
46    klo = 0
47 
48    ! find 2 levels on sigma surfaces surrounding pl at each i,j    
49 
50    do k=kts, kte-1
51       if ((p(k) >= pl) .and. (p(k+1) < pl)) then
52          khi = k+1
53          klo = k
54          exit
55       end if
56    end do
57 
58    if (klo < 1) then                                      
59       write(unit=message(1),fmt='(a,f11.3,a)') &
60          'error finding pressure level ',pconst,' mb above the surface'
61       write(unit=message(2),fmt='(a,f11.3,2x,a,f11.3)') 'pl=',pl,'  psfc=',psfc
62       call da_error(__FILE__,__LINE__,message(1:2))                                
63    end if                                                         
64 
65    ! get temperature at pl (tl), extrapolate t at surface (ts)     
66    ! and t at sea level (t0) with 6.5 k/km lapse rate              
67 
68    tlo=t(klo) * (eps+q(klo))/(eps*(1.0+q(klo)))
69    thi=t(khi) * (eps+q(khi))/(eps*(1.0+q(khi)))
70    tl=thi-(thi-tlo)*log(pl/p(khi)) &
71                       /log(p(klo)/p(khi))               
72    ts=tl*(psfc/pl)**xterm                           
73    t0=ts +gamma*terr
74 
75    ! correct sea level temperature if too hot                      
76 
77    if ( t0 >= tc ) then
78       if ( ts <= tc ) then
79         t0 = tc
80       else
81         t0 = tc-0.005*(ts-tc)**2
82       end if
83    end if
84 
85    ! compute sea level pressure                                    
86 
87    slp=psfc*exp(2.0*gravity*terr/(gas_constant*(ts+t0)))
88 
89    if (trace_use) call da_trace_exit("da_tpq_to_slp")
90 
91 end subroutine da_tpq_to_slp
92 
93