-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathadiabatic.f90
executable file
·117 lines (100 loc) · 4.54 KB
/
adiabatic.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
!>@author
!>Paul Connolly, The University of Manchester
!>@brief
!>routines for simple adiabatic parcel model
module adiabatic_vars
use numerics_type
implicit none
public
!private ::
real(wp), parameter :: r_gas=8.314_wp, ma=29.e-3_wp,mw=18.e-3_wp, &
ra=r_gas/ma,rv=r_gas/mw,cp=1005._wp, &
eps=ra/rv,lv=2.5e6_wp, ttr=273.15_wp, grav=9.81_wp
real(wp) :: p111,w111,theta111,theta_q_sat111, theta_surf,t1old,plcl1
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! saturation vapour pressure over liquid !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>calculates the saturation vapour pressure over liquid water according to buck fit
!>@param[in] t: temperature
!>@return svp_liq: saturation vapour pressure over liquid water
function svp_liq(t)
use numerics_type
implicit none
real(wp), intent(in) :: t
real(wp) :: svp_liq
svp_liq = 100._wp*6.1121_wp* &
exp((18.678_wp - (t-ttr)/ 234.5_wp)* &
(t-ttr)/(257.14_wp + (t-ttr)))
end function svp_liq
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! moist potential temperature !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>calculates the moist potential temperature
!>@param[in] t: temperature
!>@return calc_theta_q: theta_q -
!> can also be used to root-find if theta_q_sat111 set
function calc_theta_q(t111)
use numerics_type
implicit none
real(wp), intent(in) :: t111
real(wp) :: calc_theta_q
real(wp) :: ws
ws=eps*svp_liq(t111)/(p111-svp_liq(t111))
calc_theta_q=t111*(100000._wp/p111)**(ra/cp)*exp(lv*ws/cp/t111)-theta_q_sat111
end function calc_theta_q
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! lcl !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>root-finder helper for finding lcl
!>@param[in] p: pressure
!>@return find_lcl: zero if at lcl
function find_lcl(plcl)
use numerics_type
implicit none
real(wp), intent(in) :: plcl
real(wp) :: find_lcl
real(wp) :: ws,tlcl
tlcl=theta111*(plcl/100000._wp)**(ra/cp)
! calculate the mixing ratio
ws=eps*svp_liq(tlcl)/(plcl-svp_liq(tlcl))
find_lcl=ws-w111 ! w111 is the mixing ratio at the start
end function find_lcl
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! find dpdz !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!>@author
!>Paul J. Connolly, The University of Manchester
!>@brief
!>calculates the moist potential temperature
!>@param[in] z: height
!>@param[in] p: pressure
!>@param[out] dpdz: derivative of p wrt z
!> used to find pressure along adiabat
subroutine hydrostatic2a(z,p,dpdz)
use numerics_type
use numerics, only : zeroin
implicit none
real(wp), intent(in) :: z
real(wp), dimension(:), intent(in) :: p
real(wp), dimension(:), intent(out) :: dpdz
real(wp) :: t
p111=p(1)
! estimate the temperature using dry adiabat
t=theta_surf*(p111/1.e5_wp)**(ra/cp)
if(p(1)<plcl1) then
t=zeroin(t,t1old*1.01_wp,calc_theta_q,1.e-5_wp)
endif
! find the temperature by iteration
dpdz(1)=-(grav*p(1))/(ra*t)
end subroutine hydrostatic2a
end module adiabatic_vars